Linear One-dimensional Site Response Analysis in the presence of stiffness-less free surface for certain power-law heterogeneities

We revisit previous results in small-strain One-dimensional Site Response Analysis of heterogeneous soil deposits. Specifically, we focus on sites whose shear modulus distribution is described by means of a power law that yields zero stiffness at the free surface. First, we show that in some cases (which we characterize in detail) considerations of energy finitude should prevail over considerations of vanishingtractions 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.


Introduction
One-dimensional site response analysis is a popular engineering approach to estimation of the ground deformation caused by seismic events.
A particular configuration that has received much attention since the seventies [1] is characterized as: • Geometrically, an infinite 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 justification to working in frequency domain is provided by invoking Random Vibration Theory as the framework for the analysis [2].
Such a configuration 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 stiffness displayed by the soil is mediated by the overburden pressure [3], thus, in general, the deeper soil layers are more difficult to deform than the upper ones, a fact that is translated into the models as a increasing stiffness with depths. Logically, the stiffness of the first infinitesimal slice of soil at the ground would present no stiffness if one assumes that the soil behaves as a granular flow in absence of the confining pressure (soil column weight). This feature ("stiffnessless" 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 stiffness vanishes at the surface, the displacements could become infinite. 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 infinity, 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 amplifications at the free surface.
Our aim in this note is two-fold: • To characterize whether displacement-field singularities are compatible with finite 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.  The harmonic wave-propagation problem to be considered, in terms of the total displacement amplitudê u t , is [6] d dz where z is 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 X g e i 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 δ d is referred as coefficient of material damping.
For the sake of clarity, let us restate that the time variation comes given by e i t , hence both X g andû t are the amplitudes that accompany the phase: X g is a positive real number, whereasû t represents 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-finitude 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 infinity 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 effect of the vanishing stiffness 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: 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û = O(z k ) as z → 0 for some k ∈ R, and similarly for µ(z) = O(z n ), for n > 0. We shall ponder the sign of k later. 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 = 0 to some z > 0: Should k be positive, the displacement field is constrained by eq. (4a), it necessarily must be the case that k > (1 − n)/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 field goes ∼ 1 and the strain as ∼ z 1 as z approaches the ground surface at 0, satisfying this condition.
Consider we had k = −|k| < 0: the energy density in the infinitesimal soil slide must remain finite and vanish as the slice tends to zero. This does not preclude singularities in the displacement field necessarily, but it limits the range of singular behaviors that are physically admissible: given a stiffness decay ∼ z n at the free surface, both the kinetic and the strain energy remain bounded, eq. (4a) and eq. (4b), as long as  In other words, k cannot be more negative than either −1/2 or −(n − 1)/2. The least negative of the two factors sets the limit in singularity.
In summary, the singularity is constrained different depending on n: • For n > 2, the displacement field cannot scale as z −1/2 or faster (finitude of kinetic energy).
• For n = 2, idem, but due to finitude of both kinetic and deformation energy.
• For n = 1, the displacement field has to scale faster than logarithmically (zero shear at surface and finitude of deformation energy).
• For 0 < n < 1, the displacement field can not be singular as z → 0, while the singularity of the strains has to be less acute than −n to verify the traction-free condition, the controlling constraint in this case.
Towhata [4] proposed that for n > 0 (including n > 2) that the displacement field is singular as z −1/2 . Such outcome conflicts with the condition of finitude of kinetic energy eq. (4b).
This energy argumentation is inspired by Williams' [7]. For reference purposes, recall that, in a homogeneous body, the displacement field 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 Response of a system with power inhomogeneity 0 < n < 2 Consider the specific case µ(z) = µ base (z/H) n , where n > 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 infinity, hence, in a plot, the stiffness tends to zero asymptotically. When n > 1, the rate also tend to zero.
Let us rewrite the linear-momentum balance in eq. (1) in dimensionless form, expanding the first term: where See that, for n ≤ 2, this is a conventional Bessel equation, wherein there is a regular singular point at ξ = 0 and an irregular one at ξ = ∞. See also that in our case we are defining the equation over the interval ξ ∈ [0, 1].
The general solution, with no boundary conditions yet, is [9] u = ξ 1−n 2 where J m represents the Bessel function of the first kind and degree m. The constants c 1 , c 2 depend on r and n.
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: Callingũ 1 andũ 2 to the part of the general solution eq. (7) proportional to c 1 and c 2 respectively, the following scaling around the free-surface (ξ → 0) is concluded: hence one of the potential solutions goes as a constant (ũ 2 ) whereas the other one (ũ 1 ) scales proportionally to ξ, its scaling depending on the homogeneity distribution by means of 1 − n.
The fact thatũ 2 is constant entails that this part of the displacement field 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 satisfies the energy conditions discussed in Section 3, as well as the tractionless condition.
On the other hand,ũ 1 ∼ ξ 1−n means that there is potentially a singular component of the displacement field if 1 < n < 2. See how dũ 1 /dz ∼ 1/z n , thus the deformation energy, eq. (4a), scales as z 1−n , thus it is singular, thus this solution is unphysical, hence c 1 = 0. The stress-free condition would not be satisfied either, and neither would the kinetic energy condition if n > 3/2. In conclusion,ũ 1 must be discarded, thence:ũ the constant c 2 is identified from the boundary condition at the base,ũ(ξ = 1) = 1. Lastly, This solution satisfies both the mechanical and energetic conditions, yet it also holds the following: thus the strain at the top will be zero for n < 1 whereas it will be singular for n > 1, despite the displacement field remaining finite for as long as 0 < n < 2.

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]: for 1 < n < 2, as n → 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 stiffness 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 first kind when its argument is large: when x |m 2 − 1/4| [9]. Now, take eq. (11), and particularize at ξ = 0 in order to get the base-to-top amplification (note that we are not taking the modulus of this expression yet): and combine it with eq. (13), 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 A hf if this condition is easily satisfied for most values of n < 2, see Figure 3.
). When carrying out the substitution into eq. (15), the substitution into the power of r in 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 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 coefficient, δ 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) as r increases 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: first, note that the absolute value of the denominator eq. (18) can be written as 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 obtain A env , i.e., the envelope of the magnitude of eq. (18): (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 n and δ d used were such that the resonant peaks "camouflaged" the bell shape. No similar behavior occurs in homogeneous sites, as there one has similar scaling of the denominator while the numerator is just a constant.

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 configuration, one appreciates a distinctive bump in base-to-top amplification 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 stiffer, 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 first 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 final 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.1 in 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 (stiffer) layers, with any impedance ratio, as the top value of stiffness reaches zero value. Hence, this could be an artifact arising from allowing for exactly-zero stiffness at the top layer. In reality, where there is always some stiffness 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 stiffness dismisses the possibility of the continuous gradation of baseline resonances.

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 infinity, that is: where m = 1, 3, . . . denotes the first, second and subsequent resonant peaks. Given that r/(π/2) = /(πc s,base /2H) = /ω s,base , where ω s,base = πc s,base /2H refers to the fundamental frequency of the homogeneous stratum with constant base stiffness, the last result can also be expressed as where ω m represents the m-th resonance frequency of the inhomogeneous stratum whose stiffness distribution is given by µ(z) = µ base (z/H) n .
In particular, for the fundamental mode, m = 1, the estimate yields which converges to the homogeneous value in the limit as n → 0 and posits that, in the range 0 < n < 2 the fundamental frequency of the layer will be smaller than the homogeneous-based-on-bottom-properties, while 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): whence, given f (ξ) = ξ n , the last approximation assuming n not 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 if n does not tend to 2. In this case, it was found what agrees perfectly if one assumes m large enough as to being able to neglect the second addend in in favor of the first.

Equivalent homogeneous properties
One could estimate the natural frequencies trough the definition of equivalent-homogeneous stiffness or shear-wave velocity.
As discussed in previous work by the authors [10][11][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 c s,eq = 1 − n 4 c s,base , whereas for higher modes c s,eq = 1 − n 2 c s,base ,

Final remarks
In the context of 1DSR of soft-soil deposits displaying a power-law stiffness evolution with null stiffness at the free surface, this note has been concerned with providing: 1) Complementary boundary conditions in terms of energy finiteness 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 equivalenthomogeneous shear-wave velocities.
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.

Acknowledgments
The authors would like to thank Prof. Elnaz Seylabi (UNR) for her invaluable input.

Supplementary Material
A Mathematica 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 first author GitHub page github.com/jgarciasuarez. An animation of different types of frequency responses depending on stiffness profile is also included in the same repository.