Minimal Conditioned Stiffness Matrices with Frequency-Dependent Path Following for Arbitrary Elastic Layers over Half-Spaces

: This paper introduces an efficient computational procedure for analyzing the propagation of harmonic waves in layered elastic media. This offers several advantages, including the ability to handle arbitrary frequencies, depths, and the number of layers above an elastic half-space, and efforts to follow dispersion curves and flag up possible singularities are investigated. While there are inherent limitations in terms of computational accuracy and capacity, this methodology is straightforward to implement for studying free or forced vibrations and obtaining relevant response data. We present computations of wavenumber dispersion diagrams, phase velocity plots, and response data in both the frequency and time domains. These computational results are provided for two example cases: plane strain and axisymmetry. Our methodology is grounded in a well-conditioned dynamic stiffness approach specifically tailored for deep-layered strata analysis. We introduce an innovative method for efficiently computing wavenumber dispersion curves. By tracking the slope of these curves, users can effectively manage continuation parameters. We illustrate this technique through numerical evidence of a layer resonance in a real-life case study characterized by a fold in the dispersion curves. Furthermore, this framework is particularly advantageous for engineers addressing problems related to ground-borne vibrations. It enables the analysis of phenomena such as zero group velocity (ZGV), where a singularity occurs, both in the frequency and time domains, shedding light on the unique characteristics of such cases. Given the reduced dimension of the problem, this formulation can considerably aid geophysicists and engineers in areas such as MASW or SASW techniques.


Introduction
The primary objective of this paper was to investigate the propagation of elastic waves in a waveguide from a rigid strip or a rigid disk placed on the surface of a two dimensional or axisymmetric three-dimensional layered semi-infinite space.Importantly, this study aimed to achieve this without the need for extensive layer subdivisions, even for deep layers.Our research builds upon the findings presented in a prior work (reference [1]), which provided a scaled dynamic stiffness matrix for bedrock layers.While extending this bedrock problem to a half-space formulation is relatively straightforward and requires minimal additional analysis, the novelty of our work lies in its practical applications.For exceptionally deep strata, we demonstrate that a path continuation methodology can identify an internal zero group velocity resonance, and the corresponding time domain solution vividly illustrates the resonance behavior.It is worth noting that larger and more complex models would yield the same results, but what sets our approach apart is its ability to achieve these outcomes with manageable computational demands, even on standard PC desktops.
The ground is modeled as a soil profile that consists of multiple elastic, homogeneous, isotropic layers.This problem has seen much activity in the past three decades but this formulation contains a novel approach to determine a dynamic stiffness matrix.What is new, but comparable to previous studies in this area, is the computation of natural modes of free vibration for very deep strata, which may be computed easily and used for further physical analysis.These are usually used to help explain a forced response calculation due to harmonic or seismic loads in the body of the ground, wherein an example related to a harmonic rigid-surface disk load is presented.Since the formulation is related primarily to P-SV wavetypes we will drop most references pertaining to either the plane strain or axisymmetric problems, where it is understood that the formulations can be applied to both cases.At the end of the paper, a forced-vibration result presented in terms of the time and frequency domains is presented for an axisymmetric problem.
Ground vibration literature contains a wealth of publications, with much work having been conducted on the natural propagating modes in a layer over a half-space.The case of a fluid layer over a higher-velocity fluid half-space is less complicated than the elastic wave problem; it was first investigated by [2], and it is well-known as one of the first attempts to study wave propagation through layers with overly infinite domains.Although the domains lacked shear wave coupling, Ref. [2]'s thorough analysis showed that unattenuated propagation of sound would occur in a slower fluid layer characterized by the minimization of something called 'the group velocity'.The group velocity basically indicates the speed at which the amplitude or envelope wave "groups" progress in a medium.Importantly, if the group velocity is constant then the wave propagation is commonly called non-dispersive, which means that all waves progress at the same speed.In layered media, where wave speeds vary between layers, group velocity is clearly non-constant, and wave propagation analysis, although linear in nature, is complicated to analyze.
Currently, elastic wave propagation analysis is usually performed via the [3] stiffness matrix method (SMM).This approach was devised based on the earlier transfer matrix method (TMM) originally proposed over 70 years ago by [4], later corrected and elaborated by [5], and iterated further by [6,7] relatively recently.Dunkin [8] also obtained modal solutions by using the TMM.As was pointed out by Roesset, and more recently by Kausel [9], the SMM has several advantages over the TMM approach: (i) the global stiffness matrix becomes symmetric in the SMM and, as a result, lesser storage is required and also fewer operations are required for executing the analysis; (ii) various loading profiles can be easily treated; (iii) sub-structuring techniques become readily applicable.In the SMM approach, two different methods exist, namely, (i) the exact approach based on the root search method (RSM), and (ii) the approximate thin layer method (TLM) based on a quadratic eigenvalue formulation.In the exact approach, while forming the global stiffness matrix, the contribution of each layer is duly incorporated without making any approximation.However, while obtaining the solution using the RSM, since the exact approach contains a number of hyperbolic and transcendental functions the solution can be obtained only by trial and iterative procedure.Tan has also described how to increase accuracy for multilayer problems [10,11], Ba et al. have described efficient methods applied to porous layers [12], and Huang et al. have formulated a two-step solution combining solutions in the Fourier and spatial domains [13].
The model considered here demonstrates the effect of a harmonic finite load over layered strata [14].The results derived by Fourier transform are valid for any frequency and, importantly, any depth of layer.In principle, following traditional methods [15], we could have used displacement and stress-continuity boundary conditions at the bottom of the layer with equations at the ground surface to generate equations for four subsequent unknowns of stress and displacement.However, this direct approach leads to formidable numerical problems, in part due to fundamental expressions for characteristic wave functions, such as cosh or sinh, which, when employed, can have a dramatic effect on the numerical evaluation of solutions.Moreover, problems can arise due to the cancellation or division of either very small or very large numbers [16].To overcome this, ref. [17] derived a well-conditioned propagator matrix (TMM) for radially symmetric problems.In this work, though, we constructed a single stiffness matrix method (SMM) for the physical layer for plane strain problems, which conveniently avoids these difficulties.We therefore proceeded in the establishment of an original global dynamic stiffness matrix using expressions that do not cause numerical problems.
A scaled dynamic stiffness matrix for an arbitrary thick layer was derived.To achieve this, the vibration components in the wavenumber domain for layer depths and a half-space were considered and arranged into a single matrix formulation.For the forced response, the load was modeled as an infinite strip, so that the problem was plane; for the axisymmetric case, the matrices were almost identical, so for a rigid uniform disk loading the vertical and radial responses could also be calculated by an inverse transform, noting that a Hankel transfom requires inversion in the axisymmetric case.Nevertheless, the overall purpose of the present study was to present a computational method that does not suffer numerical evaluation difficulties when predicting vibration transmission, in particular its attenuation on the surface of a deep layer.The usefulness of the method is illustrated by presenting numerical results from potentially computationally intensive application examples.Comparisons against numerical solutions using finite element methods (FEM) or boundary element methods (BEM) have been successfully performed previously but we wished to avoid details of truncation, convergence, or element type that might detract from the essence of this work.
This paper, which focuses on determining solutions for free-and forced-wave propagation problems, is structured as follows: In Section 2, we establish the scaled formulation of the problem for both plane strain and axisymmetric coordinates in the wavenumber domain.This section also provides a brief overview of the root-finding methodology.Subsequently, we present a numerical example that verifies the free-vibration formulation derived by previous researchers.Additionally, in a second example, we identify the presence of zero group velocity modes, which are associated with a well-known earthquake study.Section 3 includes frequency and time domain simulations, specifically addressing a stationary rigid disk load.Finally, in Section 4, we offer a summary of our findings and draw conclusions to wrap up the paper.

Computational Model
In this section, the well-known methodologies related to free and forced vibration are formulated for plane strain or axisymmetric cases.

Two Dimensional and Axisymmetric Models
A generic model considered is shown in Figure 1.A strip load has width 2b and is aligned with respect to the z-axis.It rests on a homogeneous, isotropic, elastic layer, with material properties E (Young's modulus), ρ (mass density), and ν (Poisson's ratio).A harmonic vertical rigidly supported load acts uniformly over a strip situated above elastic layers.The elastic layers of finite depth consist of homogeneous and isotropic material, overlying an infinite half-space of flexible material.The model is two dimensional, and the co-ordinate system and parameters are shown in Figure 1.Although two dimensional, the methods used here cover plane strain and axisymmetric cases.This figure shows a generalized example of a semi-infinite stratified soil profile with linear elastic layers.For computation of theoretical dispersion curves corresponding to the assumed layer structure the problem assumes plane strain in the (x, z) plane.The x-axis is taken parallel to the layers, with the x-axis in a horizontal direction to the surface wave propagation.The positive z-axis is directed downwards.To develop an axisymmetric model, we referred to [18] for inspiration.For plane strain conditions much of the analysis necessary for the derivation of the dynamic stiffness matrix has been presented before in references [14,15,19], so this will only be briefly summarized.For plane strains the behavior of the elastic material is described by Navier's elastodynamic equations [14].Without loss of generality, in the absence of body forces, these equations apply to any layer and are written in vector form as follows: where u = (u, w) represents the components of the displacement in the x and z directions and λ, µ are the Lamé constants.For the free-vibration problem the boundary conditions for the layered problem are as follows: [(a)] zero stresses at the surface, so on z = 0, [(b)] continuity of displacement and stress at internal interfaces, and in the far field, [(c)] displacements u, w → 0 as z → ∞.The quantities C p and C s are, respectively, the pressure (P) and shear (S) wave speeds given by Note that the boundary conditions are almost identical for the axisymmetric problem in (r, z), but on the surface σ zz = − P/πb 2 whereas σ zz = −(P/π2b) over a disk of radius b or half-strip length b for the forced vibration in axisymmetry and plane strain, respectively.The non-symmetric and complex-valued matrix [K] (e,i) , Appendix B, and the steps required in its derivation are given in [1].Specifically, [K] (e,i) is the dynamic stiffness matrix for a single elastic layer that is valid for any depth H > 0. To include the half-space in the formulation, we consider the addition of a 2 × 2 complex element stiffness matrix, [K] (e,n+1) : where D = 1/(k 2 − α 1 α 2 ) and at the lowest interface u (e,n+1) = [iw n+1 , u n+1 ] T and σ (e,n+1) = [−iσ n+1 , −τ n+1 ] T .Due to the way in which the ordering of the stresses σ (e,i) and displacements u (e,i) are conveniently organized, the stresses and tractions at any interface cancel so that the load vector on the right-hand side is zero except for stress that accounts for surface load.The element matrix equations obtained for each layer of the soil model are subsequently assembled at the common layer interfaces to form the global layered system.The stiffness matrices for each layer can now be combined to give a single matrix equation for an elastic layer over an elastic half-space, involving the scaled stiffness matrix for the elastic layers [K (e,i) ], i = 1..n and the matrix for the half-space [K (n+1) ]. [ Hence, it is straightforward to generalize this technique to n elastic layers supported by an elastic half-space where the size of the dynamic stiffness matrix becomes a single complex-valued 2n × 2n matrix, but remembering that the matrix entries are transcendental in wavenumber k.The right-hand side now only involves transformed stresses at the ground surface, since the neighboring interface values have canceled each other, as discussed earlier: Kausel and Roesset, in [3], previously introduced an alternative formulation using a transfer matrix method (TMM) approach with stiffness matrices (SMM) akin to those employed in structural dynamics.An element stiffness matrix, denoted as [K i ], is derived for each layer within the geodynamic model.This element stiffness matrix for a specific layer relates the stresses at the upper and lower interfaces of the layer to the corresponding displacements.In the case of a multi-layered model, the system stiffness matrices, combined with prescribed load vectors, can be employed to solve for displacements using techniques analogous to the standard finite element method.When a load is applied to the surface, the transformed stresses {τ} G become non-zero, allowing the matrix Equation ( 6) to be solved in spatial wavenumber.
Apart from a slight change in the transformed values for the displacements, identical expressions are obtained in the axisymmetric P-SV case.Hence, it is straightforward to compute plane strain propagating and axisymmetric P-SV propagating modes in the free-vibration case.

Free Vibration
To consider free vibration for the layered half-space problem the right-hand side of Equation ( 6) is zero, and, hence, we have the following expression: where the matrix has been defined above (see Equation ( 5)) and the unknowns are the displacements on interfaces.Non-trivial solutions of Equation ( 7) are found by equating the determinant of the matrix to zero.The non-zero solutions or eigenvectors of this equation give the actual (real-valued) u and w displacements representing the propagating modes (for real wavenumber solutions) through a vertical cross-section.
For calculations for the response of a layered half-space to a surface loading, the transformed loads in the plane strain and axisymmetric P-SV are not identical nor are the inverse transforms that need to be evaluated.

Forced Vibration
In the plane strain case, the transformed load on the surface is represented as a traction.In both cases of plane strain and axisymmetry, responses at points on the surface of the ground can be recovered by employing suitable quadrature to the products of the displacements in the wavenumber domain and the transform kernel.The solution of the system of Equations ( 6) yields the transformed displacements u(k, z, ω) at the surface and, hence, by suitable inverse-transformations solutions can be obtained.The response in the space-frequency domain is obtained by an inverse transformation from the wavenumber to the spatial domain, as outlined in [18], where elsewhere.( 8)

Numerical Evaluation of the Dynamic Stiffness Matrix
In general, when dealing with non-dimensional wavenumbers where k s h ≥ 36 (with k s representing the shear wavenumber), the conventional approach encounters difficulties.Specifically, when the depth exceeds approximately six shear wavelengths (h ≥ 6λ s ), a numerical bottleneck problem arises when attempting to solve the system of linear algebraic equations.Additionally, high-frequency computations can become ill-conditioned.It is important to note that results cannot be reliably presented in cases where a "bottleneck" occurs or is imminent since matrix elements may become unbounded.
To overcome the challenges associated with the global coordinate system, we introduce a set of local coordinate systems to characterize vibrations within layered media.This approach is designed to address bottleneck issues as outlined in [16].
In this work, the necessary manipulations to provide a robust scaled matrix, as a projection method, have already been performed to assist in coding the matrix.The utilization of the projected method ensures stable numerical evaluations for stiffness matrix entries, thereby mitigating concerns related to numerical round-off errors.This is particularly crucial when dealing with operations such as dividing large numbers by small ones or subtracting very large numbers.For instance, round-off errors can lead to exponential growth in relative errors when evaluating the difference between functions like |cosh(s 1 ) − cosh(s 2 )|, where s 1 < s < s 2 and the difference |s 1 − s 2 | is small.Additionally, the numerical assessment of hyperbolic functions such as cosh and sinh can pose challenges when handling large argument values.

Root Finding and Parameter Continuation
Considering free vibration, the expression, Equation (7), may be written succintly as where each entry in the stiffness matrix [K] represents a nonlinear function in terms of the spatial wavenumber k and frequency ω, giving rise to a nonlinear root-finding exercise when non-trivial solutions are sought.While we shall not delve into the specifics of the function itself, our focus is on locating its roots.To achieve this, the method described below is employed to identify the roots of the determinant det[K], denoted as DK and commonly called the dispersion equation.A notable consequence of this approach is that the phase velocity ω/k for any propagating mode Φ can be easily derived for points along the dispersion curves where A root-finding routine has been implemented, utilizing the complex plane winding number integral to verify the count of real roots for systems such as Equation (11).These routines for root finding in the complex plane rely partly on the application of the complexargument principle.
Wave propagation frequencies or wavenumbers, along with their associated modes carrying energy, can be readily evaluated by determining the value of the roots of Equation (11).Numerous methods are available for finding roots or zeroes of nonlinear functions in one or multiple variables.However, the root-finding algorithm employed in this study is based on the procedure developed by Ivansson and Karasalo [20].Further elaboration on the specific details of this algorithm will not be provided here except to mention that no sophisticated strategy for determining the initial guesses for the solvers was adopted.For example, at a specific frequency of interest the winding number integral determines the number of roots in a domain, and a suitable nonlinear solver finds all these roots.Supposing 10 distinct roots have been found, then 10 curves are followed, using the continuation process outlined below.
Instead of seeking all the roots of the dispersion Equation ( 11) in a step-by-step manner through stepping in one of the parameters, a technique different from more traditional methods has been adopted.This approach proves effective in systematically seeking roots along curves in (ω, k)-space.The method involves tracing dispersion curves through a parametric continuation approach by using previous roots as initial estimates.In this context, the parameter can be either the frequency ω or, alternatively, the wavenumber k.For simplicity, let us assume all the roots, k i , at step i have been found and that the roots are now sought at the next frequency, ω = ω i + ∆(ω), where, in the case of constant steps, ∆(ω) = ∆ω = ω i+1 − ω i .If |∆ω| is sufficiently small, and using simple linear extrapolation to determine an initial estimate k 0 i+1 , then we should expect to find solutions, k 0 i+1 , of Equation ( 13) efficiently: This method works well while the dispersion curves are slowly varying; moreover, continuation for all the roots simultaneously is possible, but when a steep slope is encountered an adaptive technique for individual curves is required.A basic alternative formulation yielding variable (usually smaller) step sizes, ∆(ω), which are inversely proportional to the rate of change of the slope at previous steps, appeared to yield convergent results, while also flagging up possible fold points.Further research into adaptive methods, polynomial extrapolation, or arc continuation is continuing.However, this simple approach introduces an additional layer of flexibility and proves highly effective, especially when each frequency corresponds to a predetermined fixed number of single roots.Currently, we focus on tracing a single set of curves in terms of frequency, as illustrated in the upper plot, Figure 2. Once a warning occurs, interchanging the parameters ω and k in the dispersion equation means practically that fixing values of the frequency ω can lead to solutions of the dispersion equation in k (11).Hence, if |∆(k)| > 0 is small we expect to find solution curves in ω.In other words, continuing the solution path in either ω or k, as a means to obtaining solutions ω i+1 of Equation ( 13) efficiently.
The methodology of interchanging is a powerful technique that requires automating in the present algorithm.In the next section, the idea of parameter continuation is proposed as a means to increase speed for computing relevant curves.Each curve is determined by the roots of a specific formulation.

Numerical Examples
Surface waves, except for Rayleigh waves on an isotropic half-space, exhibit a phenomenon called dispersion, where the apparent velocity along the surface varies with frequency.When seismic sources generate waves, they contain a range of frequencies, and each frequency component has its own velocity, known as the phase velocity, denoted as c(ω).If we were dealing with a single-frequency (monochromatic) wave, we would only need to consider the phase velocity for that specific frequency to fully characterize the disturbance.However, in real scenarios, waves of multiple frequencies interfere with each other, creating constructive and destructive patterns that influence the overall ground motion.These constructive interference patterns behave as wave packets, which themselves propagate along the surface with well-defined group velocities, represented as U(ω).Therefore, the phase velocity is primarily determined by medium properties (such as layer thickness, intrinsic P and S velocities, rigidity, etc.) and how well a particular harmonic component fits into the associated boundary conditions.This was discussed in the previous section.On the other hand, the group velocity depends on both medium properties (through their influence on phase velocity) and how phase velocity varies with frequency.This variation controls the interference between different harmonics.The group velocity is especially significant because energy mainly propagates within the constructively interfering wave packets, which move at the group velocity rather than individual phase velocities.
Additionally, there are backward waves, which are guided waves with opposite phase and group velocities, causing the phase to move toward the wave source.The study of backward waves dates back over a century, beginning with Lamb's pioneering work in 1904 [21].A characteristic feature of a backward-wave mode is a distinctive bend in the dispersion curve with a negative slope, resulting in the negative group velocity of the corresponding guided wave.At the lowest frequency within the backward-wave range, the group velocity becomes zero while the wavenumber remains non-zero, leading to what we call a zero group velocity (ZGV) mode.
The detailed physics of situations involving negative group velocities is not discussed here but is exemplified through example results in both the time and frequency domains, drawn from the existing literature.To validate the accuracy of the method, several scenarios are compared in the examples presented.It is important to note that these scenarios are not considered benchmark solutions as no such benchmarks exist, to the best of the authors' knowledge.Instead, they serve to illustrate the formulation's efficiency in terms of computational speed.The soil characteristic values used in these examples are provided for reference in two cases.

Free Vibration-Dispersion Relations
An application of the free vibration, Section 2.2, is to aid understanding of characteristics that show up for a forced-response case.In this respect, Figures 3 and 4 show the principal results for the ground profiles, Tables 1 and 2. Verification of the method in the first example and demonstration in the second example of the "unusual" ZGV behavior mentioned above is presented.To show the accuracy and effectiveness of the proposed technique, we proceed to evaluate the dispersion relation and phase velocity spectra for the two ground profiles.

Verification of Model, [22]: Example 1
A ground profile was chosen in the thesis by [22], which was further evaluated by [23], Table 1 using various popular approaches.The ground profile is made up from four different layers, each 1.0 km in depth, and of a semi-infinite space, with the second and fourth softer layers sandwiched by stiffer layers.Generally, the stiffness of the layer and mass density vary non-uniformly with depth and the shear wavespeeds vary from C s = 2.7 km/s to C s = 4.7km/s for the half-space.There is little physical information included in this referred example but the computations are quite extensive due to the depths of the layers involved, the present method requiring inversions of 10 × 10 matrices; hence, this example was chosen for computational reasons.Figure 3 was produced by the root solver described in Section 2.5, and by digitization of the curves presented in Kumar and Naskar it was possible to reproduce the phase velocity curves presented.
It was not possible to directly compare the accuracy or computation times between Kumar and Naskar and the Kausel method for this single example, but working on an i5 processor 32MB RAM machine running MATLAB 2020A calculation software the tasks were clocked.Here, wavespeed variation with frequencies at 140 equal frequency steps, between 0 and 9 Hz, for the first nine natural modes were computed.Adopting the rootfinding method beginning at 9 Hz and ending at 0.1 Hz, the method adopted here took approximately 12.0 seconds to complete the 140-frequency-step analysis.The half-space shear wavespeed is shown as a line across the graph, to show the limit of the solution.It is clear that all three methods agree well.2, used for synthetic seismograms observed as aftershocks of the Latur earthquake, 1993, was evaluated in paper [6].It consists, simply, of two layers overlying an elastic half-space.Figure 4 is a more conventional dispersion diagram for wavenumbers and phase velocities and is suitable for physical interpretation, such as cut-on frequencies and Leaky modes, etc. Figure 4a, however, does not feature characteristics that allow straightforward interpretation.Two modes have real wavespeeds or wavenumbers at all frequencies, but only one, which we shall call the second mode, tends toward the Rayleigh wavespeed for the upper layer at low frequencies (190 m/s).At low frequencies or large wavelengths, the layer depth should become negligible and wavespeeds should tend towards the half-space Rayleigh wavespeed (which would exist without the two layers).At higher frequencies the modes do tend toward the upper-layer wavespeed as expected.Employing the path continuation algorithm described in Section 2.5 was straightforward and required no adaptive methods to adjust step sizes.  2 in [6], where → indicates a ZGV (zero group velocity) location and the circles indicate frequencies (28.5 Hz and 38.8 Hz) and wavenumbers where modal shapes are displayed for each (see Figure 5).(a) Propagating wavenumber plotted against frequency for ground profile, Table 2. Zero group velocity detected at 28.4 Hz and 38.7 Hz; (b) Phase velocity plotted against frequency for ground profile, Table 2. Zero group velocity detected at 28.4 Hz and 38.7 Hz.
It is well known in the non-destructive inspection community that Lamb modes can possess backward-wave phenomena for certain values of Poisson's ratio or combinations of layer and substrate.A backward-wave mode is indicated by a typical bend of the dispersion curve, with a negative slope resulting in the negative group velocity of the corresponding guided wave.At the minimum frequency of the backward-wave range, the group velocity becomes zero while the wavenumber remains non-zero; the corresponding wave is referred to as a ZGV mode, which occurs twice in the frequency range used here (Figure 4a).Aside from the opposite phase movement, for which detection is a rather challenging task, the backward wave is considered to manifest itself by a surface-displacement resonance at the ZGV frequency.The red dots in Figure 4b depict the positions at which the propagation modes are computed in the next section.Note well that employing the path continuation algorithm described in Section 2.5 requires an ad hoc adaptive technique that can reduce the step size sufficiently while finding distinct wavenumber solutions.2. All modes were evaluated at (a) Six mode shapes plotted against depth on a logarithmic scale at frequency f = 28.5 Hz (see Figure 4), (b) Six mode shapes plotted against depth on a logarithmic scale at frequency, f = 38.8Hz (see Figure 4).

Modal Evaluation, Latur Earthquake: Example 2
Figure 5 shows the variation of amplitude (to within an arbitrary amplitude factor) of the vertical and horizontal motion components of the first few modes, plotted against depth up to 1000 m into the ground at frequencies in the vicinity of the appearance of the zero group velocity (ZGV) mode.The depth scale is logarithmic, which highlights the modal behavior nearer the surface.The legends in this figure show the phase velocity and the wavenumber related to individual mode shapes.In Figure 5a, the variation of the radial mode (U) with depth, adjacent to the first ZGV frequency f = 28.5 Hz, shifts from maximum to minimum within the top 5 m layer then attenuates quickly to zero through the lower strata in the ground.At the surface, it appears that each modal wavetype can equally contribute to a surface-vibration response, except at the second mode c W = 2302 m/s, k = 0.078 rad/m, where it appears, on the contrary, that this second mode dominates in carrying the vertical surface wave energy for vertical mode shapes (W).At this frequency, mode shapes in the vicinity of the ZGV turning point are not remarkably different to other modes corresponding to other waves.At the second ZGV turning point, the mode shapes were evaluated at an adjacent frequency, 38.8 Hz (Figure 5b).Contributors to potential surface waves are evenly distributed through the lowest phase velocity modes across all seven mode shapes displayed here.Interestingly, the mode that relates to the half-space shear wave dispersion curve shows a potential significant contribution in the lower depths of the strata, up to around 600 m below the surface.

Forced Vibration-Frequency Domain, Latur Earthquake, Example 2
For this section, the results displayed in Figure 6 were obtained with the proposed solution strategy determined by the axisymmetric model, Section 2.3, and are compared with those obtained with a methodology proposed by Wang, [6], a method which has been adopted widely by many researchers in the ground vibration and seismic scientific communities.This approach evolved from the Thomson-Haskell layer transition matrix methods, with an orthonormalization procedure to overcome numerical stability issues with the original propagation algorithm.In the following, responses to a harmonic surface load were calculated using the hysteretic damping values established in Wang [6].To conduct this, a circular uniform traction load with radius b = 1 m, Equation ( 8), was applied on the surface of a layered soil deposit.The deposit consisted of an h = 5 m soft soil layer on top of a 300 m deep stiffer layer which overlaid a homogeneous half-space.The soil layer properties are shown in Table 2. Figure 6 shows the vertical and horizontal displacement magnitude at five response points located on the soil surface at distances d = 5, 10, 20, 40, 60 m from the disk center, as obtained with the proposed solution strategy.Attenuation of radial and vertical displacement with distance is evident throughout the frequency range.At 10 m and beyond it is clear that the greatest response featured close to frequencies 18 and 38 Hz for radial components and 20 and 40 Hz for vertical components beyond 10 m.Although this is not strong evidence of resonance due to the existence of a ZGV mode, an increased amplitude was observed and verified the conclusions by Wang [6].

Forced Vibration-Time Domain, Latur Earthquake, Example 2
To study the behavior of a time-dependent force located on the surface of the layered half-space we used the numerical fast Fourier transform (FFT) to convert solutions from the frequency domain to the time domain.This technique gives the displacements in the space-time domain, taking a source whose temporal variation is provided by a Ricker wavelet, as defined below.To reduce the computational time we used quadrature in space to calculate responses at a few selected receiver points.Also, since the Ricker wavelet decays rapidly, in both time and frequency, this allowed a computed time domain solution to be calculated in a reasonable time for interpretation.
The Ricker wavelet function is given by where W 0 is the amplitude, τ = (t − t d )/t 0 ; t d is the time delay at which maximum occurs, and t 0 is related to the "natural frequency" of the Ricker wavelet, 1/t 0 = 2π f 0 (see Figure 7).In our calculations we chose the width of the input pulse as f 0 = 14 Hz, somewhat arbitrarily, but to allow at least the cut-on of a few propagating modes.By applying exponential windowing with damping terms equivalent to the frequency domain, time domain responses were calculated using an FFT with 512 frequency steps up to the Nyquist frequency around 250 Hz.The two sets of results highlight the difference between a model that includes a soft top layer (see Figure 8a) and one which does not (see Figure 8b).Both sets of solutions naturally exhibit a strong dispersive nature as the envelope of the solution wave packet propagates along the surface.This is especially evident for the soft 5 m layer result: there appears a second wave packet in the vertical response at distances that develop from 40 m, which is due to the soft-layer response.As reported in [6], the "slow" speeds of the two wave packets effected by ZGV dispersion are clearly evident.

Conclusions
A uniformly applicable model has been developed for studying the propagation of surface vibrations through elastic layers of arbitrary depth.This model comprises an elastic, isotropic, and homogeneous layer overlaying a half-space.A dynamic stiffness matrix, wellconditioned for this model, has been established through the projection of characteristic functions onto depth-dimension endpoints.To explore the behavior of natural propagating modes within these layers under zero-surface-stress conditions, singular values of the dynamic stiffness were computed.The characteristic equation was also employed to determine the solution shapes of these propagating modes.For instance, when examining a ground profile with a high Poisson's ratio in the top layer over a deep strata, a zero group velocity mode was identified.This mode corresponded to a layer resonance, which had been previously observed.Vibrational responses due to a fixed disk load in the frequency domain indicated resonances at cut-off frequencies.Additionally, time domain results were obtained, confirming the existence of non-dispersive wave packets, as suggested by earlier researchers.This formulation's broad applicability for dynamic stiffness matrices opens the door to modeling various new problems where plane strain conditions are applicable.Furthermore, it can be readily expanded to account for sub-layers of different materials within the strata and extended to three-dimensional problems, where an exact 6x6 dynamic stiffness matrix could be established.
In summary, the following conclusions have been drawn: • A well-conditioned dynamic stiffness matrix has been developed for layers of arbitrary depths.

•
The characteristic equation has been formulated in a way that requires unknowns to be expressed only at material layer nodal points, eliminating the need for subdividing the strata.

•
The model has been applied to known ground profiles, revealing the existence of zero group velocity modes.Examples involving time histories and frequency domain responses have been provided.

•
Anisotropic constitutive models, such as transversely isotropic or orthotropic models, may be used to describe the mechanical behavior of anisotropic soils in numerical analyses, and this is an area in which we are interested in extending our knowledge for future research.

•
The appearance of zero group velocity modes may be studied efficiently in a sensitivity analysis, since the size of the problem is small and the continuation methodology can aid in testing its existence.• Future advancements towards a 3D dynamic stiffness matrix, involving the analytical determination of 36 entries, could significantly reduce computational costs in full-sized models well within reach on small desktop PCs.

Figure 2 .
Figure 2. Sketch showing a turning point or fold/branch bifurcation.A white dot represents the initial guess for the solver, the red dots show the true solution, and the yellow dots show the ZGV turning point.The black dots represent previously computed solutions.

Figure 3 .
Figure3.A comparison between the dispersion phase velocity curves calculated by the present method (solid colored lines), with[23] (•), and by[3] (□).The ground profile comprises four layers of depth 4.0 km each situated over an elastic half-space, Table1, referenced as Profile C in Kumar and Naskar[23].

Figure 4 .
Figure 4. Dispersion curves of the present method for ground profile, Table2in[6], where → indicates a ZGV (zero group velocity) location and the circles indicate frequencies (28.5 Hz and 38.8 Hz) and wavenumbers where modal shapes are displayed for each (see Figure5).(a) Propagating wavenumber plotted against frequency for ground profile, Table2.Zero group velocity detected at 28.4 Hz and 38.7 Hz; (b) Phase velocity plotted against frequency for ground profile, Table2.Zero group velocity detected at 28.4 Hz and 38.7 Hz.

Figure 5 .
Modal shapes evaluated at two frequencies related to the Latur earthquake ground profile, Table

Figure 7 .
Figure 7.A Ricker wavelet time profile, which applies to a uniform pressure over a rigid disk in Example 2.

Figure 8 .
For a Ricker wavelet load over a rigid disk, responses at 5, 10, 20, 40, and 60 m from center of disk load.Upper (a) and lower (b) solutions compared for two ground profiles.(a) Ground profile represented in Table 2, without the top 5 m layer; (b) Ground profile represented in Table 2 .
Figure 1.A layered half-space model geometry.The parameters of the model are shear and pressure wave velocities C s , C p , respectively, in each layer, with mass density ρ i and layer thickness h i .