Seismic Response of Inhomogeneous Soil Deposits with Exponentially Varying Stiffness

Available analytical solutions for the seismic response of continuously inhomogeneous media reveal the possibility of higher amounts of seismic energy reaching the ground surface over media with discontinuous variations in elastic properties with depth (Towhata 1996). Problems involving continuous inhomogeneity have been treated for different velocity models (Ewing et al. 1957; Ambraseys 1959; Toki and Cherry 1972; Schreyer 1977; Gazetas 1982; Dakoulas and Gazetas 1985; Aki 1993; Towhata 1996; Semblat and Pecker 2009; Paraschakis et al. 2010; Rovithis et al. 2011; Mylonakis et al. 2013; Vrettos 2013; Garcia-Suarez 2020) without necessarily acknowledging the importance of continuous versus discontinuous variations in soil material properties with depth. For the common case of a parabolic (power law) velocity model, the special condition of vanishing stiffness at the soil surface has been examined by Dobry et al. (1971), Towhata (1996), Travasarou and Gazetas (2004), Rovithis et al. (2011), and Kausel (2013). More recently, Mylonakis et al. (2013) demonstrated analytically that under viscoelastic conditions and regardless of excitation frequency and material damping, the near-surface shear strain may be either zero or infinite (but never finite) depending solely on rate of inhomogeneity. This result is in contrast to surface displacements that, naturally, are always finite.

The effect of soil inhomogeneity has been studied in multiple dimensions for different types of stress waves. Vrettos (1990) tackled the problem of dispersive horizontally polarized shear surface waves propagating in an inhomogeneous half-space for both bounded and unbounded variations in the shear modulus with depth. A collection of available solutions (not restricted to elastic waves) for wave propagation in inhomogeneous media is available in Manolis and Shaw (1997, 1999, 2000) and Brekhovskikh and Bayer (1976).

The seismic response of inhomogeneous media modeled using “equivalent” homogeneous counterparts, that is, defined by means of an average wave propagation velocity, has been explored against exact solutions (Idriss and Seed 1967; Dobry et al. 1971; Zhao 1997; Rovithis et al. 2011; Garcia-Suarez 2020). Additionally, simple analytical formulae have been proposed for the fundamental dynamic characteristics of inhomogeneous soils at resonance conditions using Rayleigh approximations (Mylonakis et al. 2013) and asymptotic analysis (Garcia-Suarez 2020; Garcia-Suarez et al. 2020).

Despite research on the subject for more than half a century, nearly all analytical solutions for inhomogeneous media are limited to cases in which stiffness and shear wave propagation velocity vary as a power of depth. Although some dependable solutions exist for other variations in stiffness (e.g., the constant minus the exponential function considered in the seminal study by Vrettos 2013), these solutions are semianalytical in nature because they are expressed in the form of an infinite power series and not special functions (e.g., of the Bessel type) associated with a wealth of mathematical properties that facilitate their analytical treatment. Closed-form analytical solutions, including relevant asymptotic expressions and approximations pertaining to low and high frequencies, can shed light on the physics of the problem, which is often not possible with other methods. The gap in knowledge regarding the dynamic behavior of more general velocity profiles provided the motivation for the herein reported work.

In this paper, the seismic response of an inhomogeneous soil layer resting on a rigid base is investigated in the realm of one-dimensional viscoelastic wave propagation theory. An exponential function is adopted, which can describe both positive and negative velocity gradients. The latter may be relevant for cases of large confining stresses under building foundations, stiff overconsolidated surface crusts, or ground improvement near the soil surface, which may result in soil stiffness decreasing with depth (Dobry et al. 1976; Vrettos 2013). For simplicity and as a first approximation, soil mass density and material damping are assumed to be constant with depth. The problem is treated analytically, leading to an exact solution expressed in terms of Bessel functions for the vibrational characteristics of the system, the base-to-surface transfer function, and the variation in ground displacements, shear strains, and curvatures with depth. The solution is compared against Rayleigh approximations for the fundamental frequency of the layer by implementing different shape functions for the corresponding mode shapes. Asymptotic and approximate solutions are proposed for the response in low and high frequencies, which are combined through a smooth transition function to provide a full-domain approximation. The model is validated against a numerical layer transfer-matrix (Haskell-Thomson) formulation for two actual soil profiles. Apart from its intrinsic theoretical interest, the proposed exact solution can be used to assess other relevant solutions for wave propagation and site response analysis.

The problem considered in this work refers to a continuously inhomogeneous viscoelastic soil layer of thickness $H$, excited by vertically propagating seismic shear ($SH$) waves applied at its base (Fig. 1).

Soil mass density, $\rho$, and hysteretic damping ratio, $\xi$, are assumed constant with depth, whereas the shear wave propagation velocity follows the exponential variationwhere $V_0$ = shear wave propagation velocity at the ground surface; $z_r$ = reference depth; and $\alpha$ = dimensionless inhomogeneity parameter.

Rearranging Eq. (1), parameter $\alpha$ can be defined aswhere $V_H(=V_0{\mathrm{e}}^{\alpha })$ = shear wave propagation velocity at the base ($z_{\mathrm{r}}=H$). Evidently, positive and negative values of $\alpha$ correspond to stiffness increasing ($V_0<V_H$) and decreasing ($V_0>V_H$) with depth, respectively. For $\alpha =0$, the model degenerates to a uniform velocity distribution ($V_0=V_H$).

The input motion is specified at the base of the system in the form of a harmonic horizontal displacement, $u=u_0\exp (i\omega t)$, with $\omega =(2\pi f)$ as the cyclic excitation frequency, thereby generating shear waves propagating vertically up and down (on reflection) in the profile.

Evidently, the increase in the shear wave propagation velocity according to this model is superlinear (when $\alpha >0$), whereas the corresponding variation in uniform normally consolidated soils is known to be sublinear (Kramer 1996; Towhata 1996). Although this discrepancy poses difficulties in simulating deep uniform soil profiles, the model can successfully describe certain classes of nonuniform soil profiles with discontinuous variations in shear wave velocity with depth via curve fitting. Examples are provided later in this paper. In addition, the model can be used to simulate arbitrary layered velocity profiles (e.g., via a transfer matrix formulation as shown in the Supplemental Materials).

Under harmonic, steady-state oscillations, one-dimensional shear waves in an inhomogeneous soil layer are described using the familiar ordinary differential equation (Ewing et al. 1957; Towhata 1996)

After introducing the displacement potential $Y(z)$ through the transformation $dY/dz=u(z)$ and integrating with respect to $z$, the governing equation in Eq. (3) can be rewritten in canonical formin which the term associated with the first derivative of the new dependent variable, $dY/dz$, has disappeared, and the depth-varying shear modulus has escaped differentiation. In the latter equationstands for the depth-varying wavenumber of the associated stress waves, where $k_0=(\omega /V_0)$ is the value at ground surface. The solution to Eq. (4) is (Wylie and Barret 1986; Trachanas 2004—see Supplemental Materials)where $J_0()$ and $N_0()$ are the zero-order Bessel functions of the first and the second kind, respectively; and $C_1$, $C_2$ = integration constants to be determined from the boundary conditions.

When differentiating with depth, the general solution for displacements is obtained from Eq. (6) aswhere $J_1()$ and $N_1()$ denote the corresponding Bessel functions of the first order. The associated shear stresses can be derived by inspection because $\tau (z)=G(z)\gamma (z)=G(z)du/dz=G(z)d^{2}Y/d{z}^2=-\rho {\omega }^{2}Y$ to obtainwhich reveals that the displacement potential in Eq. (6) is merely a scaled form of shear stress $\tau (z)$ on division by the product ($-\rho {\omega }^2$). A relevant discussion is provided in Towhata (1996).

The boundary condition of a traction-free surface, $\tau (0)=0$, yields a proportionality relation between constants $C_1$ and $C_2$ (see Supplemental Materials); the solution for displacements becomes

From Eqs. (3) and (9), the corresponding solution for shear strain is

Finally, operating on Eq. (3), one can easily establish the second derivative $d^{2}u/d{z}^2$, which expresses the curvature of soil with depth

The last equation, which holds for any smooth velocity model, is particularly useful when shear systems, such as the one at hand, are in contact with embedded flexural structural elements, such as piles or diaphragm walls, to ensure that strains and curvatures are “transferred” from the soil to the structural elements. Relevant discussions on strain and curvature transmissibility pertaining to piles are provided in Mylonakis (2001), Anoyatis et al. (2013), Di Laora and Rovithis (2015), Di Laora et al. (2017) and in a recent review article by Di Laora and Rovithis (2021).

Eq. (11) holds for general variations in shear wave propagation velocity with depth and is not limited to a particular type of soil inhomogeneity. For the velocity model in Eq. (1), one obtains

In light of these analytical developments, the following observations are worth noting. First, the wavenumber [$k_o\exp (-\alpha z/z_r)$] appearing in the arguments of the Bessel functions coincides with the variable coefficient $k(z)$ in Eq. (4), which is not the case with other relevant solutions (e.g., Semblat and Pecker 2009; Rovithis et al. 2011; Mylonakis et al. 2013). Second, the multiplier ($z_r/\alpha$) in the arguments of the Bessel functions normalizes the wavenumber to produce a dimensionless argument that coincides with the derivative $d/dz[-k(z)]=d/dz[-k_o\exp (-\alpha z/z_r)]$. Again, this special property of the particular type of inhomogeneity is not encountered in other related solutions. Third, the multipliers $\exp (-2\alpha z/z_r)$ and $\exp (-\alpha z/z_r)$ appear, respectively, in the expressions for strains in Eq. (10) and displacements in Eq. (9), but not in the one for stresses [Eq. (8)]. For the common case of positive $\alpha$, this implies that strains and displacements attenuate with depth faster than stresses, whereas strains attenuate with depth faster than displacements. This, again, is different from other basic problems in geomechanics, such as the Boussinesq point-load solution, for which stresses attenuate with depth faster than displacements (Poulos and Davis 1974; Davis and Selvadurai 1996). The slower attenuation of stresses relative to strains and displacements in the problem at hand can be interpreted in light of the direct dependence of stresses on dynamic equilibrium with the imposed inertial loads, which forces stresses to exhibit a more “statically determined” behavior than the other response functions. Fourth, for positive values of $\alpha$, soil curvature ($1/R$) attenuates rapidly with depth, as evident from the exponential multipliers $\exp (-2\alpha z/z_r)$ and $\exp (-3\alpha z/z_r)$ in the first and second term in the right-hand side of Eq. (12). As the term associated with soil displacement attenuates faster than the one for strain, curvatures at depth tend to be governed by strains. Fifth, the soil curvature ($1/R$) in an inhomogeneous medium may be larger or smaller than the corresponding one in a homogeneous medium subjected to the same displacement and wavenumber depending on the sign of the first term in the right side of Eq. (11). Sixth, the rate of inhomogeneity, $\alpha$, does not affect the order of the Bessel functions in the solution. This, again, is different from other relevant problems involving wave propagation in inhomogeneous media (Towhata 1996). Seventh, shear strain at the soil surface, $\gamma (0)$, is always zero because shear stress is by definition zero, whereas the corresponding shear modulus, $G(0)$, is finite.

For $z=0$, the displacement at soil surface is obtained as (see Supplemental Materials)

The fundamental natural frequency of the system can be well approximated by the familiar Rayleigh quotient (Semblat and Pecker 2009; Mylonakis et al. 2013)

This expression is derived in the Supplemental Materials. By replacing $u(z)$ by a unitary dimensionless shape function, $\psi (z)$, representing the fundamental mode shape of the inhomogeneous layer, ${\mathrm{\Phi }}_1(z)$, considering constant mass density and employing the normalized shear wave propagation velocity ${\overline{V}}_s(\eta )=V_s(\eta )/V_{\mathrm{H}}=e^{+\alpha (\eta -1)}$, the Rayleigh quotient can be rewritten in the normalized form

Following Mylonakis et al. (2013), four shape functions are examined (Table 1): linear, parabolic, sinusoidal, and “compatible.” The last function is obtained by imposing a distributed horizontal static load proportional to the soil mass along the height of the column (see Supplemental Materials). The Rayleigh solutions are compared against the numerically evaluated exact solution in Fig. 4 for different values of the inhomogeneity factor $\alpha$. Specific values of the ($f/f_{1\mathrm{H}}$) ratio are provided in Table 2, whereas the exact fundamental mode shape in Eq. (16) is compared with these shape functions in Figs. 5(a–c) and (d–f) for $V_{\mathrm{o}}/V_{\mathrm{H}}$ ratios lower and higher than unity, respectively.

Note that in the special case in which $\alpha \to 0$, the compatible shape function in Table 1 reduces to the parabolic one. The normalized natural frequency for the linear shape is then obtained as $(12{)}^{1/2}/\pi \cong 1.1027$; the parabolic and stiffness compatible shapes yield $(10{)}^{1/2}/\pi \cong 1.0066$, whereas the sinusoidal shape yields the exact value of 1 (Kloukinas et al. 2012). These results are provided in Table 2.

Two actual profiles with discontinuous variation in shear wave propagation velocity with depth were employed to assess the proposed solutions against a numerical simulation of soil response using a series of a piecewise homogeneous layers, implemented through the classical Haskell-Thompson algorithm that involves a sequence of transfer-matrix multiplications (Thompson 1950; Roesset 1977). To this end, a real soil profile (hereafter referred to as “Profile A”) was selected from the Kiban-Kyoshin Network (KiK-net) database (Okada et al. 2004; NIED 2019) in the Tagahaki site in Ibarakiken Prefecture, where the IBRH13 station of the Kik-net accelerometric network is installed. This site is classified as “Level Ground” (LG), according to the taxonomy proposed in Thompson et al. (2012), which allows for the application of one-dimensional wave propagation theory in a site response analysis. The shear wave propagation velocity profile is plotted in Fig. 10(a) up to an elevation of $-34\mathrm{m}$, where a rigid base is assumed following the definition of seismic bedrock in EC8.

The second soil profile (hereafter referred to as “Profile B”) is a multilayer deposit typical of the San Francisco Bay Area, resting on bedrock at a depth of approximately 60 m (Gazetas and Dobry 1984), Fig. 10(b). The available data, as reported in these references, are provided in Tables 1 and 2 of the Supplemental Materials. In both case studies, the soil mass density and hysteretic damping ratio were assumed constant with depth, equal to $2\mathrm{\mu g}/{\mathrm{m}}^3$ and 5%, respectively. The exponential function velocity model was established by employing a least squares fit. In this manner, the vector ($V_0$, $V_0/V_{\mathrm{H}}$) was established at approximately ($197\mathrm{m}/\mathrm{s}$, 0.303) and ($134\mathrm{m}/\mathrm{s}$, 0.265) for Profiles A and B, respectively. Accordingly, the dimensionless inhomogeneity parameter $\alpha$ equals 1.194 and 1.326 in the two cases; the exponential velocity profile is given by (Fig. 10)

Worth mentioning is that the $V_s$ data reported in the Kik-net and other databases are typically crude approximations of the actual soil properties. Nevertheless, the intention of this study is to validate the proposed approximate/asymptotic solutions against “exact” solutions for a layered soil rather than evaluate the level of approximation associated with the databases. Studies encompassing uncertainty in site response analysis can be found in Stewart et al. (2017), Stewart and Afshari (2021), and de la Torre et al. (2022).

An analytical solution was derived for the dynamic response of a soil layer with stiffness varying exponentially with depth. The main conclusions drawn from the study are as follows:

As a final remark, it is fair to mention that the increase in wave velocity with depth in the model at hand is superlinear (when $\alpha >0$), whereas the corresponding variation in uniform, normally consolidated soils is sublinear. This discrepancy inevitably poses difficulties in simulating deep uniform soil profiles; however, the model can successfully describe certain classes of nonuniform soil profiles with discontinuous variation in shear wave velocity with depth via curve fitting (Fig. 10). In addition, the solution at hand can be readily extended to encompass an arbitrary number of inhomogeneous layers or a combination of stacked homogeneous and inhomogeneous layers by means of multiplications of transfer matrices. This is outlined in the Supplemental Materials.

All data and models generated or used during the study appear in the published article. All code that supports the findings of this study are available from the corresponding author on reasonable request.

The authors acknowledge the assistance of Ms. Vasiliki Miha and Ms. Kathleen Schulze, former students of the senior author at the University of Patras and the City University of New York, respectively, during the early stages of this research. Thanks are also due to Dr. Charalambos Paraschakis for fruitful discussions and to Dr. Elia Voyagaki for proofreading the manuscript. Partial funding was received by EU/H2020 under grant agreement number 730900 (SERA) with George Mylonakis as Principal Investigator for University of Bristol.