On the fundamental resonant mode of inhomogeneous soil deposits

The problem of estimating seismic ground deformation is central to state-of-practice procedures of designing and maintaining infrastructure in earthquake-prone areas. Particularly, the problem of estimating the displacement ﬁeld in a soft shallow layer overlying rigid bedrock induced by simple SH wave excitation has been favored by engineers due to its simplicity combined with inherent relevance for practical scenarios. We here derive analytical, accurate estimates for both the fundamental frequency and the amplitude of the ﬁrst resonant mode of such systems by applying an intuitive argument based on resonance of single-degree-of-freedom systems. Our estimates do not presuppose a continuous velocity distribution, and can be used for fast assessment of site response in seismic hazard assessment and engineering design. On the basis of the said estimates of fundamental frequency and amplitude, we next propose a novel deﬁnition of “equivalent homogeneous shear modulus” of the inhomogeneous deposit; and we show that the response of the fundamental mode of these systems is governed by the mechanical properties of the layers closer to the bedrock. We ﬁnally discuss the validity of our argument, and evaluate the accuracy of our results by comparison with analytical and numerical solutions.


Introduction
It has been long established that site conditions and soil layering have meaningful effects on ground motion characteristics during seismic events. The variation of soil mechanical properties with depth can induce a concentration of seismic energy in the softer, shallow layers [1], which in turn can aggravate the seismic risk for man-made above-ground and underground structures.
In the past 50 years, numerous researchers have studied the modification of seismic shaking by inhomogeneous stratified media, usually by means of simplified one dimensional (1D) wave propagation models under the assumption of linear-elastic medium [2,3,4,5,6,7,8,9]. These assumptions can be shown to be appropriate when 1) the soil is subjected to small strains ergo the material response remains linear-elastic (absence non-linear material phenomena), 2) the soil is horizontally stratified (thus the problem reduces to one-dimensional, vertical, wave propagation), 3) the excitation comprises of vertically-propagating SH waves, and 4) the material (intrinsic) damping is frequency-independent.
The customary approach focuses on finding a steady-state solution in the frequency domain (frequently referred to a site response transfer function), which happens to be the dominant component of the response, except under some specific circumstances (low-frequency fundamental mode, high-frequency excitation [10]).

The governing equations
The harmonic wave-propagation problem that we consider, in terms of the total displacement amplitudeû t , is [11] where z is the coordinate stretching from the free surface to the rigid base, ρ represents the density of the soil (which is assumed to remain 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.
Hysteretic damping is not shown explicitly in eq. (1) because the shear modulus is interpreted as a storage modulus (real number) plus a loss modulus (imaginary number) as µ(z) → µ(z)(1 + iξ d ), where ξ d is referred to as coefficient of material damping or as loss factor.
For the sake of clarity, let us restate that the time variation is 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).
The prior equation is classified as an homogeneous ODE of the Helmhotztype with inhomogeneous boundary conditions. By introducingû =û t − X g , the relative-to-base horizontal displacement amplitude, and y = H − z, the complementary coordinate to z, the problem (1) becomes (see thatẌ which represents the same phenomenon, although in a subtly, yet consequential, different way: instead of a homogeneous equation, now a forced Helmhotz equation with homogeneous boundary conditions is obtained. This fact will come in handy later.

Modeling the distribution of mechanical properties across depth
The main difficulty associated with the solution of either one of the equations above is that real soil profiles have complex variation of properties with depth.
Consider µ(y) = ρc s (y) 2 = ρc 2 s,base f 2 (η), where c s (y) is a function that supplies the value of the S-wave velocity at any point in the stratum; c s,base corresponds to the value at the base; η is the dimensionless vertical coordinate (normalized with respect to the profile thickness, H); and f (η) is the dimensionless function that describes the distribution of material properties along the deposit (see that f (0) = 1 and f (1) = β = c s,top /c s,base ). This difficulty hinders the procurement of simple closed-form solutions for the natural frequencies and the displacement field for arbitrary profiles.
A widely-used expression for f (η) for which the problem has an analytical solution is the so-called "generalized-parabola evolution", given by where b = β 1/n = (c s,top /c s,base ) 1/n and n is dubbed "inhomogeneity factor" (also referred to as "inhomogeneity parameter" [12]). The parameter b controls the maximum variation of c s , which in this case always corresponds to top-to-bottom or vice versa, since the function only allows for monotonic evolution across the layer.

Previous results by Mylonakis, Rovithis and Parashakis (2011, 2013)
Mylonakis and collaborators [13] proposed applying the Rayleigh quotient of the governing ODE eq. (1), to estimate the natural frequencies: where ω s represents the fundamental frequency of the system, η = y/H is a dimensionless vertical coordinate, ψ 1 is a shape function that must be pointwisesimilar to the fundamental shape in order for the estimation to work, and f (η) = c s (y)/c s,base is a dimensionless S-wave velocity evolution along the layer.
The method relies on "guessing" a good approximation of the first mode shape, but simple approximations such as parabolic or sinusoidal, can yield good results, at least in those models conforming to the eq. (3) (check Figure 6 in order to gauge these assertions).
The same authors also provided an exact solution of eq. (2) when f (η) is given by eq. (3) [11]. They found that the total displacement evolution was given to a homogeneous layer [14]: Still focusing specifically on profiles given by eq. (3), the exact amplitude of the base-to-top transfer function was also provided by simply particularizing the eq. (5) at η = 1: a function which depends on the parameters introduced with eq. (5). Therefore, the exact natural frequencies of the inhomogeneous system where f (η) statisfies eq. (3) are the roots of the denominator of eq. (7) (namely, the bracketed factor).

Previous results by Garcia-Suarez, Seylabi and Asimaki (2019)
Likewise, an estimate particularly well-suited for the high-frequency regime has been furnished previously by the authors [15]. Regarding the ground surface amplification relative to the profile base, it was found A asy (β, n, r) = c s,base /c s,top cos cs,eq/H , where c s,eq represents the harmonic mean of the distribution c s (z): These results were derived through analysis of the exact solution eq. (7). The relation between the high-frequency limit of the exact solution and the so-called nominal frequency of the site [6] (the fundamental frequency of the homogeneous system having constant wave velocity given as eq. (9)) was established mathematically, which allowed to specify the validity conditions for the approximation to hold. The natural frequencies of the site (particularly the higher modes) were shown to be approximated as for k = 1, 2, 3, . . .

Derivation of estimates for the fundamental mode response
Given the system depicted in Section 1, governed by eq. (1), we would like to find estimates for the displacement field and the natural frequency concerning not the whole range of possible dynamic behavior but only the fundamental (viz. the first) resonant mode. To fulfil this endeavor, the following approach is proposed.

An "easing" argument based on single-degree-of-freedom resonance
See that eq. (2) represents but a classic equation of a one-degree-of-freedom system, wherein elastic, damping, and inertial forces balance an external load [16].
It may not be so easy to appraise this as the damping is not of viscous type but hysteretic, so let us show it clearly by explicitly substituting µ(y) by µ(y)(1 + iξ d ) in eq. (2) (at this point let us also allow for evolving density, ρ = ρ(y), too): Now it is easy to recognize each term in the left-hand side: from left to right, the first addend corresponds to elastic forces, the second one to hysteretic damping forces and the third to inertial forces, while the right-hand side is the "external body force". The diagrams in Figure 2 represent two possible admissible configurations in eq. (11): • To the left, a configuration wherein elastic forces are dominant and hence, for the most part, these take care of balancing the external load (lowfrequency regime). See that elastic force magnitude is proportional to the gradient of displacements (strains), and thus large displacements are not necessarily required to balance the excitation, only enough change from point to point. Damping and inertial forces remain relatively small. See also that the phase difference between inertial and elastic forces is π rad, that is, half a period.
• To the right, a configuration we may refer as resonance: the elastic forces and the inertial forces reach equal magnitude, and, since they are out of phase, cancel each other, leaving the damping force solely in charge of compensating the external load. Given that ξ d tends to be small (ξ d 1 usually), it follows that the gradient magnitude must surge to compensate the detrimental influence of ξ d , but in this situation the magnitudes of both displacement and displacement gradient are directly linked through the condition of inertia balancing elastic forces, and thus the displacement amplitude surges as well. In conclusion, let us, from this point onward, consider resonance as the dynamic setting wherein elastic and inertial forces cancel each other (as they attain the same magnitude but their phases are shifted by π radians), and thus damping forces alone balance the "external" loading. In such a scenario we may particularize eq. (11) as This equation still represents a second-order ODE, hence it can be integrated to yield the horizontal displacement field when = ω s while enforcing the actual boundary conditions in (2). Note how in these developments we do not require that the soil properties satisfy continuity or smoothness conditions, but only that the final function can be integrated. Proceeding with the integration and imposing the boundary conditions in (2) one attains the displacement at resonanceû (y) where F (y) represents where y represents a dummy integration variable. This expression can be evaluated for any arbitrary layered medium, regardless of the impedance contrast across layers. If one next assumes, as it is customary, that the density is constant, then µ(y)/ρ = c 2 s (y) = (c s,base f (η)) 2 , and so is the damping ξ d too, eq. (13) can be simplified as:û where, recall, f (η) = c s (y)/c s,base . It still remains to identify the value of the fundamental frequency, ω s , to complete the previous expression. The estimate provided by the Rayleigh quotient, eq. (4), could be used. Hence, assuming once again that the stiffness varies with depth, whereas the density and the damping are constant across the stratum, we propose the following estimate for the fundamental frequency: This estimate is completely general, but it still requires a good guess of the modal shape (what may be difficult to ascertain in some cases, e.g., stiff sites with interbedded soft layers).

A dichotomy arising from the resonance argument
There is a consequential presupposition that we should point out before going any further: eq. (11) represents equilibrium point-wise, that is, the forces therein are forces at a point on the cross-section, not the total force acting on the cross-section under consideration. The fact that at a given depth inertial and elastic forces balance each other, does not entail necessarily that the same happens at all depths simultaneously. In other words, it is not guaranteed that one can integrate eq. (12) to obtain eq. (13) and subsequent, and therefore this result must be considered a simplification.
There is a simple way of realizing this limitation: since elastic and inertial forces are equal, eq. (12) can be also written as From this, one can derive an expression for the displacement alternative to eq. (16):û which, again, may hold at the point level but it cannot be extended to the whole layer as it would not satisfy the boundary conditions at the base (unlike eq. (12), eq. (17) is not a second-order ODE). One can hence readily see that eq. (16) (after integrating across the layer) and eq. (18) cannot hold simultaneously at all points of the vertical section.
Let us pause briefly to appreciate the meaning of eq. (18): it indicates that a point that is undergoing resonance experiences a relative displacement equal to the displacement at the base divided by the damping factor, and a change of phase of half a period. See that both eq. (18) and eq. (16) assign the same phase lag to such "resonating points".

Fundamental frequency and equivalent shear modulus
The dichotomy discussed in section 2.1.1 yields two contradictory displacement expressions, eq. (15) and eq. (18). One could think of reconciliating them, and such equality would yield an expression for the fundamental frequency. Such approach does not provide a functioning expression, but it inspires an alternative strategy nonetheless: consider the resonance amplitude, in terms of relative displacement in the homogeneous stratum from eq. (6b). The fundamental resonant amplitude at the free-surface can be approximated as [17] u(η = 1) where it is assumed that ξ d 1. Then, equate this expression to the modulus of eq. (15) to find ω s and the equivalent homogeneous shear modulus so that the inhomogeneous deposit resonates as a homogeneous layer : whence Thus the inhomogeneous resonating system behaves as an homogeneous one of an equivalent shear modulus equal to a linearly-weighted harmonic mean of the inhomogeneous distribution. See that these linear weights (π 3 z/16 ≈ 1. This result differs fundamentally from the one that was obtained when considering the long-wavelength limit, eq. (9) , where the equivalent homogeneous property was the harmonic mean of the shear-wave velocity. Note, however, that both results agree for the case of a homogeneous layer.

Ground surface amplification relative to bedrock during resonance
Let us proceed to use eq. (16) to obtain the amplification during resonance, ignoring the phase difference: Note how this results points out that the whole heterogeneity distribution affects the amplification. In order to check this estimate, consider a homogeneous stratum, f (η) = 1, and the exact mode shape sin(πη/2): where 1 has been neglected as it will be much smaller than the second term, since ξ d 1 is presupposed. Comparing to the value suggested by Röesset [17], that is, 4/πξ d ≈ 1.27/ξ d , it seems that, as it also happens in the case of viscously-damped 1-dof systems [18], the maximum amplitude is slightly shifted towards a higher frequency due to the presence of damping. Nevertheless, the estimate, at least in this case, yields a value just about a 3% smaller.
Incidentally, the expression ω s in eq. (20b) could be used to furnish another estimate, but this would yield simply the already-mentioned homogeneous value for amplification proposed by Röesset, namely eq. (19).

Comparison between exact and estimated response at resonance
for profiles idealized as "generalized parabolas" In the forthcoming sections, we shall focus our attention to verifying the estimates against results that assume an S-wave distribution given by eq. (3).

Fundamental frequency
As a token of suitability, assume a profile given by eq. (3) and consider

Ground surface amplification relative to bedrock
A comparison of the performance of both eq. (22) with eq. (4) and eq. (19) with eq. (20b) against eq. (7) and eq. (8) is desirable. Before being able to deliver it, the integral in eq. (22) must be evaluated using f (η) as given in eq. (3): The expression for the Rayleigh Quotient (assuming parabolic shape ψ 1 = η 2 /2 − η) is available in the literature [13]. Figure 7 displays the comparison.
The estimates for the fundamental mode amplitude work fairly well in all three cases, as reported in the following table: Parameters Equation this amplitude misfit is associated with an overestimation of the fundamental frequency, hence improving the estimate for the latter (e.g. using a sinusoidal mode shape instead of a parabolic one) would yield a better estimate of the former. The data regarding the fundamental frequency estimation are already consigned in Figure 6.
One should highlight that using eqs. (19) and (22) in conjunction with eq. (4) for the fundamental mode, while using eq. (8) for the higher modes, may provide decent estimates of the whole range of frequency behavior for these systems in some cases.

Particularization for profiles with acute gradients of stiffness localized in the near-surface
In many practical cases, there is evidence of substantial stiffness changes being confined in a narrow region close to the surface, followed by a milder stiffness increase at depth. This behavior adjusts better to values of the inhomogeneity factor n ≤ 0.1.
Let us re-assess eq. (22) bringing this fact to bear. a) In the first place, as R depends on the integral of the mode shapes, if most inhomogeneity is intensely localized in a thin shallow layer, it seems logical to assume that the shape will resemble the one of the homogeneous profile except in that shallow section of the profile. However, as the quotient entails a global comparison to the exact shape (through the integral), as opposed to point-wise local comparison, we argue that the value of the integrals will not differ too much from the one of the homogeneous stratum with properties as at the base, and hence R ≈ π/2. As test to this argument, revisit Figure 4, the intermediate tile. Let us deploy this arguments within eq. (22): This is the same result that we obtained in eq. (23), that is, when we used the formula to calculate the amplification corresponding to the homogeneous layer. See that this expression does not require a prior estimate of the natural frequency, unlike eq. (22).
The conclusion is that, for the case when n is relatively small, the behavior of the inhomogeneous stratum shall resemble the one of the homogeneous one having as mechanical properties those of the base of the layer. Under these premises, the response of the fundamental mode can be considered as independent of both the contrast, i.e., β, and the fine details of the distribution of shear-wave velocity, i.e., exact value of n (provided that n < 0.1 to begin with). In order to better appraise this corollary, revisit the results in Figure 4 and Figure 5 in light of this premise, particularly the middle tile in Figure 4.

Verification for profiles idealized as "generalized parabolas"
For instance, compare the value that the estimate would deliver (δ d = 0.05) for the fundamental mode to those of eq. (7), which are shown in next  Table 2: Exact amplitude at the fundamental mode for a number of models, from eq. (7).
On the same time, the estimate from eq. (25) is All the prior results concern ξ d = 0.05. The maximum error, is 24.3% (amplitude is underestimated), but it corresponds to n = 0.9 and β = c s,top /c s,base = 0.1, a 90% decrease in stiffness happening almost linearly, hence the premises upon which the reasoning was developed (gradient localized at the top) do not hold in any way, shape or form for this case, which is nonetheless an non-physical variation of stiffness for natural deposits.
Is interesting to note that for β = c s,top /c s,base = 0.5 we still obtain less than 10% error (underestimation), a fact that we ascribe to the helpful effect of the weight in the integral factor.
Likewise, as long as n = 0.1, the error is just testimonial. Again, the agreement is credited to the weights of the integral favoring the bottom layers.
Clearly, the trends that were predicted from eq. (25) are the ones observed.

Verification for discontinuous velocity profiles
Hitherto, everything has been referred to profiles under the umbrella of the "generalized parabola", eq.     Let us point out that these three sites are classified as " LG" according to the taxonomy proposed by Thompson and collaborators [20], that is that onedimensional site response analyses can inform realistic estimates of amplification are these sites (provided realistic data regarding damping at specific sites was used). In addition, consider the last layer in each case as rigid bedrock, then calculate the transfer function numerically following the classic methodology (give, for instance, in [14]). The calculation result is contained in Figure 8, the damping values used have been ξ d = 0.147 for IBRH17, 0.068 for TKCH08 and 0.035 for IBRH10.
As it can be seen from inspecting Figure 8, the value of the natural frequency is remarkably well estimated by eq. (20b), not only in the second case (TKCH08) which is the one that would better adjust to n 1, but also the other two.
The performance of eq. (20b) estimating the fundamental frequency of the sites and is assessed qualitatively in the following  Table 6: Amplitudes at the fundamental mode for three KiK-net model in Figure 8 (the number in parenthesis represent percent deviation from numerical result).
Thus the amplitude error ranges from less than 1% (overestimation) to less than 8% at most (underestimation). In similar fashion, the error in (under)estimating the fundamental frequency is around 2% of the numerically-obtained value.

Conclusions and future work
The main objective of this work has been to deliver simple estimates for the response of an inhomogeneous soil layer overlying rigid bedrock, subjected to incoming SH waves which excite the fundamental resonance mode of the stratum.
Customary simplifying assumptions (linear-elastic response, vertically-layered medium, SH wave-front, frequency-independent damping) are taken to be in effect.
The expression eq. (22) is valid for soils having any stiffness variation with depth (although it requires a good frequency estimate R), whereas eq. (25) is adequate to estimate such amplitude for any distribution, so long as the dominant stiffness gradients happen relatively close to the surface.
Regarding estimation of the fundamental frequency of the site, eq. (20b) has proved to work satisfactorily in both scenarios concerning continuous profiles and sites presenting discontinuous profiles with sudden jumps.
The first path to improvement must be replacing the rigid bedrock assumption with a more realistic elastic bedrock [14], which would extend the applicability of the results from a limited number of strata with large impedance contrast (effective rigid bedrock) to any two horizontal sections in a layered half-space.
A two-layer system resting on bedrock [11], another classic configuration, may also be readily used to expand these results (intuition hints that the estimate eq. (25) will perform satisfactorily in these cases, based on the argumentation carried out in section 4).
Insofar mechanical properties distribution is concerned, one should re-stress that the expression eq. (22) is completely general, so actually any velocity profile could be used, including those that present stiffness reversal (soft intermediate layers), but it requires an estimate of the fundamental frequency (e.g. the Rayleigh quotient), what may not be always easy to obtain. On the other hand, it also remains to be verified that eq. (20b) works in scenarios presenting either increasing stiffness gradients (β > 1) or intermediate softer layers.
Once again, the previous results were derived under very specific suppositions, see section 1. One of those presuppositions is linear-elastic behavior. The possibility of extending the argument used in section 2 to the non-linear regime via equivalent linear assumptions is enticing and should be explored.

Supplementary material
A Jupyter notebook containing code to generate the result figures in this paper can be found in the first author GitHub (github.com/jgarciasuarez).
In the repository fundamental_mode_resonance_inhomogeneous_soil one can find the notebook after execution and a release containing the corresponding .ipynb file. Using this notebook in conjunction with Binder is advised.