We revisit previous results in small-strain One-dimensional Site Response Analysis of heterogeneous soil deposits. Speciﬁcally, we focus on sites whose shear modulus distribution is described by means of a power law that yields zero stiﬀness at the free surface. First, we show that in some cases (which we characterize in detail) considerations of energy ﬁnitude should prevail over considerations of vanishing tractions at the free-surface, as these may pose acuter constrains. We re-evaluate previous contributions in light of this result. Second, we analyze the previously-reported occurrence of “energy accumulation in upper layers”, providing a physical explanation for it. In passing, we supply estimates of the natural frequencies, and compare these with our previous results. Keywords— Geotechnical Engineering Site Response Analysis Vanishing stiﬀness 1 Introduction One-dimensional site response analysis is a popular engineering approach to estimation of the ground deformation caused by seismic events. A particular conﬁguration that has received much attention since the seventies [1] is characterized as: Geometrically, an inﬁnite soft soil layer overlying rigid bedrock. There can be inhomogeneity along depth, but all vertical cross-sections are considered equal (1D problem, plane strains). Impinging S waves (at zero incidence angle) that can be modeled as a forced-displacement boundary condition at the bedrock level. The strains in the soil are taken to remain small enough as to approximate material behavior as linear-elastic. The depth of the deposit (width of the layer) is small enough so as to assume that steady-state conditions are rapidly achieved, so that, in conclusion, the wave reverberation disembogues into the appearance of of a standing wave so the problem can shift into the frequency domain. An alternative justiﬁcation to working in frequency domain is provided by invoking Random Vibration Theory as the framework for the analysis [2]. Corresponding author: ajgarcia@caltech.edu 1 Such a conﬁguration combines relative simplicity, what make it suitable to be tackled with analytical techniques, with relevance for practical engineering purposes. It is a known fact that the stiﬀness displayed by the soil is mediated by the overburden pressure [3], thus, in general, the deeper soil layers are more diﬃcult to deform than the upper ones, a fact that is translated into the models as a increasing stiﬀness with depths. Logically, the stiﬀness of the ﬁrst inﬁnitesimal slice of soil at the ground would present no stiﬀness if one assumes that the soil behaves as a granular ﬂow in absence of the conﬁning pressure (soil column weight). This feature (“stiﬀnessless” of ground surface) leads to some striking consequences. Towhata [4] studied the response of a deposit when the inhomogeneity is given by a power law. It was concluded that if the stiﬀness vanishes at the surface, the displacements could become inﬁnite. Conversely, Gezetas and Trevasarou [5], building on Dobry’s results [1], observed that: “The shear strain at the surface would either tend to zero or to inﬁnity, depending only on the rate of increase of the shear wave velocity with depth”. More recently, Rovithis et al. [6] found results that seemed to contradict Towhata’s conclusions, at least for a certain range of exponents of the power law. In the same article, they also reported a certain “accumulation phenomenon”, leading to large displacement ampliﬁcations at the free surface. Our aim in this note is two-fold: To characterize whether displacement-ﬁeld singularities are compatible with ﬁnite energy levels, in similar spirit of Williams’ work in Fracture Mechanics [7]. These criteria would rule out the physical feasibility of some results in [4]. To revisit the problem addressed in [6], in order to provide a physical explanation of the occurrence of “surface energy-accumulation” therein detailed. 2 The physical model cs,top cs,basecsy u = 0𝑦τxz = 0 HzX ሷg cs,top cs,basecsy 𝑦τxz = 0 Hz ut=Xgeiϖt Figure 1: Scheme of the system considered in this study (consider cs;top= 0). The harmonic wave-propagation problem to be considered, in terms of the total displacement amplitude ^ut, is [6] d dz (z)d^ut dz +$2^ut= 0s.t. ^ut(z=H) =Xg,lim z!0^ xz= lim z!0 (z)d~u dz = 0;(1) 2 wherezis the coordinate stretching from the free surface to the rigid base, represents the density of the soil (which is assumed to be constant), the seismic S-wave excitation is modeled as a uniform displacement at the soil base Xgei$t, and(z)represents the variation of shear modulus with depth. As per usual, the density is taken as constant throughout the deposit. Despite damping not being not shown explicitly in eq. (1), the shear modulus could readily be interpreted as an storage modulus (real number) plus a loss modulus (imaginary number) by means of a substitution (z)!(z)(1 + id), wheredis referred as coeﬃcient of material damping. For the sake of clarity, let us restate that the time variation comes given by ei$t, hence both Xgand^utare theamplitudes that accompany the phase: Xgis a positive real number, whereas ^utrepresents a complex number whose modulus is the total displacement magnitude and whose argument represents the time shift between load (input) and displacement (response, output). 3 Energy-ﬁnitude as boundary condition at the free surface We have chosen to expresses the shear-wave velocity as a limit in order to stress that, technically, the shear strain could tend to inﬁnity as z!0, as long as the decay of (z)compensates the singularity, satisfying in the end the traction-free condition [6]. We would like to propose a complementary approach to considering the eﬀect of the vanishing stiﬀness of the ground, one focused on energetic conditions instead of mechanical ones. Let us consider per-cycle-averaged strain energy density and the kinetic energy density [8] for this purpose: W(z) =1 2(z) d^ut dz 2 =1 2(z)d^ut dz conjd^ut dz ; (2) K(z) =1 2$2j^utj2=1 2$2^utconj(^ut); (3) where conj(x)represents the complex conjugate of x. 3.1 Scaling of the deformation and the kinetic energy at the surface Let us assume scaling ^u=O(zk)asz!0for somek2R, and similarly for (z) =O(zn), forn>0. We shall ponder the sign of klater. Next, given these assumptions, consider both the total strain energy per unit of area as well its kinetic counterpart in a slice of soil running from z= 0to somez>0: Zz 01 2(z) d^ut dz 2 dzznz2(k1)z=zn1+2k; (4a) Zz 01 2$2j^utj2dzz2kz=z1+2k: (4b) Shouldkbe positive, the displacement ﬁeld is constrained by eq. (4a), it necessarily must be the case that k>(1n)=2, otherwise the deformation energy would be singular. In practice we see this to be the case, for example, in the homogeneous case ( n= 0), wherein the displacement ﬁeld goes 1and the strain as z1aszapproaches the ground surface at 0, satisfying this condition. Consider we had k=jkj<0: the energy density in the inﬁnitesimal soil slide must remain ﬁnite and vanish as the slice tends to zero. This does not preclude singularities in the displacement ﬁeld necessarily, but it limits the range of singular behaviors that are physically admissible: given a stiﬀness decay znat the free surface, both the kinetic and the strain energy remain bounded, eq. (4a) and eq. (4b), as long as k>max 1 2;n1 2 : (5) 3 -1.0 -0.5 0.0 0.5 1.001234 knTractions Deformation Kinetic Tractions∧Deformation∧KineticFigure 2: Graphical representation of the conditions that energy ﬁnitude and tractionless surface impose over scalingO(zk)of displacement ﬁeld as a function of stiﬀness scaling around surface O(zn). Admissible region does not include boundary. In other words, kcannot be more negative than either 1=2or(n1)=2. The least negative of the two factors sets the limit in singularity. In summary, the singularity is constrained diﬀerent depending on n: Forn>2, the displacement ﬁeld cannot scale as z1=2or faster (ﬁnitude of kinetic energy). Forn= 2, idem, but due to ﬁnitude of both kinetic and deformation energy. For1<n< 2, the displacement ﬁeld cannot scale as z(n1)=2or faster (ﬁnitude of deformation density). Forn= 1, the displacement ﬁeld has to scale faster than logarithmically (zero shear at surface and ﬁnitude of deformation energy). For0< n < 1, the displacement ﬁeld can not be singular as z!0, while the singularity of the strains has to be less acute than nto verify the traction-free condition, the controlling constraint in this case. Towhata [4] proposed that for n>0(includingn>2) that the displacement ﬁeld is singular as z1=2. Such outcome conﬂicts with the condition of ﬁnitude of kinetic energy eq. (4b). ThisenergyargumentationisinspiredbyWilliams’[7]. Forreferencepurposes, recallthat, inahomogeneous body, the displacement ﬁeld at a crack tip scales with the square root of the distance to the tip, while the strain and the stress does it proportionally to the inverse of the square root. 4 4 Response of a system with power inhomogeneity 0< n < 2 Consider the speciﬁc case (z) =base(z=H)n, wheren>0. Any exponent will make (z= 0) = 0 , yet the rates of change around the free-surface vary greatly. If 0<n< 1, the rate tends to inﬁnity, hence, in a plot, the stiﬀness tends to zero asymptotically. When n>1, the rate also tend to zero. Letusrewritethelinear-momentumbalanceineq.(1)indimensionlessform, expandingtheﬁrstterm: d2~u d2+n d~u d+1 n~u= 0; (6) where=z=H2[0;1],r=$H=cs;base,^ut=~u=Xgandf() =(z)=base. See that f(= 1) = 1 , f(= 0) = 0 andf0(= 0)>0. See that, for n2, this is a conventional Bessel equation, wherein there is a regular singular point at = 0and an irregular one at =1. See also that in our case we are deﬁning the equation over the interval2[0;1]. The general solution, with no boundary conditions yet, is [9] ~u=1n 2 c1J1n 2n2r 2n1n 2 +c2J1n 2n2r 2n1n 2 (7) whereJmrepresents the Bessel function of the ﬁrst kind and degree m. The constants c1,c2depend onr andn. The local behavior of the solution around the free surface is described by the asymptotics as the argument of the Bessel function goes to zero. Recall: J1n 2n2r 2n1n 2 1n 2; (8a) J1n 2n2r 2n1n 2 n1 2: (8b) Calling ~u1and~u2to the part of the general solution eq. (7) proportional to c1andc2respectively, the following scaling around the free-surface ( !0) is concluded: ~u11n;~u21; (9) hence one of the potential solutions goes as a constant ( ~u2) whereas the other one ( ~u1) scales proportionally to, its scaling depending on the homogeneity distribution by means of 1n. The fact that ~u2is constant entails that this part of the displacement ﬁeld attains a local maximum at the ground surface, thus the there is no straining of the upper soil slice, thus this part of the solution satisﬁes the energy conditions discussed in Section 3, as well as the tractionless condition. On the other hand, ~u11nmeans that there is potentially a singular component of the displacement ﬁeld if 1<n< 2. See howd~u1=dz1=zn, thus the deformation energy, eq. (4a), scales as z1n, thus it is singular, thus this solution is unphysical, hence c1= 0. The stress-free condition would not be satisﬁed either, and neither would the kinetic energy condition if n>3=2. In conclusion, ~u1must be discarded, thence: ~u= ~u2=c21n 2J1n 2n2r 2n1n 2 ; (10) 5 the constant c2is identiﬁed from the boundary condition at the base, ~u(= 1) = 1 . Lastly, ^ut Xg=1n 2J1n 2n 2r 2n1n 2 J1n 2n 2r 2n: (11) This solution satisﬁes both the mechanical and energetic conditions, yet it also holds the following: ^ut^ut(z= 0) Xg2n!d~u dz1n; (12) thus the strain at the top will be zero for n < 1whereas it will be singular for n > 1, despite the displacement ﬁeld remaining ﬁnite for as long as 0<n< 2. 4.1 Accumulation phenomena The prior scaling is the same as the one presented by Rovithis, Parashakis and Mylonakis in [6]. These researchers also reported a certain phenomenon, which was also noted by Gazetas and Travasarou in [5]: for1<n< 2, asn!2, the amplitude of the soil deformation appears ever-increasing as the frequency of the excitation increases. A milder version was shown to occur also when the ratio between top and base stiﬀness is very small, yet not zero. Even more, even for cases relatively removed from the limit value n= 2, in spite of not appreciating this ever-increasing amplitude, the resonance spikes seems to stand over a baseline level well above 1. We aim to explain this feature next. First, let us recall the asymptotic behavior of Bessel functions of the ﬁrst kind when its argument is large: Jm(x)r 2 xsin xm 2 4 ; (13) whenx jm21=4j[9]. Now, take eq. (11), and particularize at = 0in order to get the base-to-top ampliﬁcation (note that we are not taking the modulus of this expression yet): A(r) =^ut(= 0) Xg= r 2n1n 2n  1 2n J1n 2n 2r 2n; (14) and combine it with eq. (13), A(r)Ahf(r) =p r 2n n 2(2n)  1 2nsec2r 2nn 2n 4 ; (15) where “hf” denotes a high-frequency (short-wavelength approximation). Recalling the condition for eq. (13), we note that the exact transfer function shall converge to Ahfif r=$H cs;base 2n 2 1n 2n2 1 4 ; (16) this condition is easily satisﬁed for most values of n<2, see Figure 3. According to eq. (15), the exponent of rwill always be positive, and the larger the closer ngets to 2. But the eﬀect of damping remains unaccounted for in the expressions. Recall the substitution !(1 + id), 6 eq.(16)rhs0 0.5 1 1.5 20.00.20.40.60.81.0 nFigure 3: Right-hand side of the inequality (16) which amounts to r!r=p1 + idr=(1 + id=2)r(1id=2). When carrying out the substitution into eq. (15), the substitution into the power of rin the denominator changes the magnitude almost imperceptibly, while having an impact on the phase lag between ground deformation and bedrock input. Conversely, the substitution into the secant function does yield important consequences over the magnitude of the response: see that sec2r 2nn 2n 4idr 2n = (17a) =1 cos 2r 2n n 2n  4 cosh dr 2n + i sin 2r 2n n 2n  4 sinh dr 2n; (17b) which decays exponentially fast as r!1due to the hyperbolic functions in the denominator brought about by the inclusion of damping. Let us consolidate these results before moving forward: in the presence of hysteretic damping, eq. (15) becomes A(r)Ahf(r) =p r 2n n 2(2n)= 1 2n cos 2r 2n n 2n  4 cosh dr 2n + i sin 2r 2n n 2n  4 sinh dr 2n:(18) By means of eq. (18) one can realize that there is an interplay between a tendency by the numerator to increase the amplitude, according to a certain power dependent on n, and the denominator which tends to decrease the amplitude exponentially yet it is abated by the small value of the damping coeﬃcient, d. Eventually, the hyperbolic functions overcomes the initial disadvantage due to the small value of d, and these will start to dominate the numerator, as the hyperbolic functions scale faster than any power of r. This means that the envelope of the magnitude of A(r)asrincreases is the one of a bell curve, initially increasing as the behavior is dominated by the numerator to fall to zero as the denominator takes over. Such baseline envelope can be calculated explicitly: ﬁrst, note that the absolute value of the denominator 7 eq. (18) can be written as cos2r 2nn 2n 4 coshdr 2n + i sin2r 2nn 2n 4 sinhdr 2n = =s 1 2 cos4r 2nn 2n 2 + cosh2dr 2n :(19) The cosine addend of the expression gives the resonance peaks as the hyperbolic cosine yields the baseline envelope. Retaining the latter plus ignoring the contribution of the damping to the magnitude of the numerator to obtainAenv, i.e., the envelope of the magnitude of eq. (18): Aenv(r) = r 2n n 2(2n)  1 2nvuut2 cosh 2dr 2n: (20) The reason behind this behavior going unnoticed and unreported hitherto, we surmise, may be due to the transfer function visualizations being either limited to a range of frequencies before the pick was attained, or the the combinations of nanddused were such that the resonant peaks “camouﬂaged” the bell shape. 8 0 1 2 3 4 5051015 a) d=0.10, n=0.33 A()Ref. [4] Exact eq.(14) Envelope eq.(20) 0 0.5 101 y/H 0 1 2 3 4 505101520 b) d=0.10, n=1.25 A() 0 0.5 101 y/H 0 1 2 3 4 5 /s,base0102030405060 c) d=0.25, n=1.67 A() 0 0.5 1 (y) / base01 y/HFigure 4: Comparison of eq. (7) to results in [6], when there is small but ﬁnite stiﬀness at the top (cs;top=cs;base = 103in (a) and (b), 108in (c)) 9 0 2 4 6 8012345 a) d=0.30, n=0.40 A()Exact eq.(14) Asymptotics eq.(15) Envelope eq.(20) 0 0.5 101 y/H 0 2 4 6 8020406080100120140 b) d=0.05, n=1.50 A() 0 0.5 101 y/H 0 2 4 6 8 /s,base0100200300400500 c) d=0.20, n=1.75 A() 0 0.5 1 (y) / base01 y/HFigure 5: Depiction of the accumulation phenomenon for various combinations of inhomogeneity exponent nand damping d. 10 No similar behavior occurs in homogeneous sites, as there one has similar scaling of the denominator while the numerator is just a constant. 4.1.1 Physical interpretation The phenomenon can be described, fair and square, in mathematical terms. However, the physics of the underlying mechanisms, not so much. We propose a tentative qualitative explanation based on resonance of layered deposits. In such a con- ﬁguration, one appreciates a distinctive bump in base-to-top ampliﬁcation around certain frequencies in addition to resonance peaks that can take place simultaneously. In those cases, there is a simultaneous “quasi-resonance” in some layers that prompts overall large amplitudes. Quasi-resonance is meant in the sense that the lower, and stiﬀer, one almost resonates as if it was on its own, prompted by the bedrock displacement imposed at its base; at the same time, a similar circumstance happens to some upper layer, it just so happens that the ﬁrst layer plays the role of the bedrock now, exciting the upper layer at a frequency at which inertial and elastic forces are out of phase [10].The ﬁnal outcome of the chained quasi-resonant deformation in each layer is a global resonance mode in the stratum that seems to raise from a baseline level above one. In the author’s recent article, [10], we identify the phenomenon described above with the bump present around$=!s;base = 3:1in the corresponding Figure 8(c). We posit that a similar phenomenon happens here, but in this case is continuous instead of localized. The continuity of the behavior obeys to the fact that there is a continuous gradation of potential upper (softer) and lower (stiﬀer) layers, with any impedance ratio, as the top value of stiﬀness reaches zero value. Hence, this could be an artifact arising from allowing for exactly-zero stiﬀness at the top layer. In reality, where there is always some stiﬀness or some extra physics entering the picture, there can be localized instances of the discussed per-layer quasi-resonance giving rise to global extra deformation, but the presence of limit stiﬀness dismisses the possibility of the continuous gradation of baseline resonances. 4.2 Estimates for the natural frequencies of the stratum See that eq. (15) allows a simple closed-form estimate of the natural frequencies, given by the values that make the secant function tend to inﬁnity, that is: cos2r 2nn 2n 4 = 0!r= 2h m 1n 2 +n 4i ; (21) wherem= 1;3;:::denotes the ﬁrst, second and subsequent resonant peaks. Given that r=(=2) = $=(cs;base=2H) =$=!s;base, where!s;base =cs;base=2Hrefers to the fundamental frequency of the homogeneous stratum with constant base stiﬀness, the last result can also be expressed as !m !s;base=m 1n 2 +n 4; (22) where!mrepresentsthe m-thresonancefrequencyoftheinhomogeneousstratumwhosestiﬀnessdistribution is given by (z) =base(z=H)n. In particular, for the fundamental mode, m= 1, the estimate yields !1 !s;base= 1n 4; (23) which converges to the homogeneous value in the limit as n!0and posits that, in the range 0<n< 2the fundamental frequency of the layer will be smaller than the homogeneous-based-on-bottom-properties, while 11 not smaller than one half of the value. This results concurs with previous work by the authors [10] in a more general setting. It was found that (see eq.(21b) in the reference, adapted to the current notation): !1 !s;base3 16Z1 0 f()1=2 ; (24) whence, given f() =n, !1 !s;base3 162 2n1=2 r 1n 2= 1n 4+O(n2); (25) the last approximation assuming nnot tending to 2. See how the smaller n, the more localized the change around the free surface, and the better the approximation, as spelled out in [10]. Likewise, for the higher modes ( m> 1) the geometrical optics approximation, presented in [11], also holds well ifndoes not tend to 2. In this case, it was found !m !s;basemZ1 01 f()1=21=2 =m 1n 2 ; (26) what agrees perfectly if one assumes mlarge enough as to being able to neglect the second addend in in favor of the ﬁrst. 4.2.1 Equivalent homogeneous properties One could estimate the natural frequencies trough the deﬁnition of equivalent-homogeneous stiﬀness or shear-wave velocity. As discussed in previous work by the authors [10 –12], coming up with equivalent homogeneous properties valid over the whole frequency range appears an irreducible task. Instead, it seems more feasible to assess equivalent homogeneous properties for either the fundamental mode [10], or relatively high frequencies in general [11]. The current scenario characterized by (z)=base= (z=H)n,0<n< 2, is no exception to this appraisal. Thus, from section 4.2, we recommend as equivalent-homogeneous shear-wave velocity cs;eq= 1n 4 cs;base; (27) whereas for higher modes cs;eq= 1n 2 cs;base; (28) 5 Final remarks In the context of 1DSR of soft-soil deposits displaying a power-law stiﬀness evolution with null stiﬀness at the free surface, this note has been concerned with providing: 1)Complementary boundary conditions in terms of energy ﬁniteness at the free surface, and using these to discuss the physical viability of some previously-reported results [4]. 2)A mathematical explanation of the accumulation phenomenon described in [6], plus a preliminary qualitative interpretation of its onset. 3)Providing estimates for the natural frequencies of these sites, and, by extension, also of its equivalent- homogeneous shear-wave velocities. 12 The future work should be concerned with extending these results to the case n>2, still in the small-strain regime, and a general re-assessment the previous conclusions in the presence of non-linear material behavior (large-amplitude excitation regime). The argumentation developed in Section 4.1.1 can and should be substantiated more profusely. 6 Acknowledgments The authors would like to thank Prof. Elnaz Seylabi (UNR) for her invaluable input. 7 Supplementary Material AMathematica notebook [13] containing the computations leading to all the plots shown in the text can be found in the repository named 1DSR_stiffnessless_surface in the ﬁrst author GitHub page github.com/jgarciasuarez . An animation of diﬀerent types of frequency responses depending on stiﬀness proﬁle is also included in the same repository.