Head Waves in Modiﬁed Weiskopf Sandy Medium

: In this study, a modiﬁed anisotropic elastic Weiskopf model for a sandy medium was introduced, satisfying the highest rotational symmetry compatible with a Weiskopf sandiness parameter greater than the unity. The developed approach was applied for computing and comparing head SP waves propagating along the free surface of the halfspace. The comparison revealed a substantial discrepancy in the SP wave pulses detected at the points of observation. The developed models for the generation and detection of SP waves are based on the FEA (ﬁnite element analysis) of the inner Lamb problem.


Weiskopf Sandy Model
The constitutive equation for modeling sandy materials was proposed by Weiskopf (1945) [1]. The model reflects an experimental observation [2][3][4][5][6][7] that sandy materials, or more generally, granular materials, exhibit a ratio between Young's and shear moduli larger than that predicted by Hook's law: where E and µ are Young's and shear moduli, respectively, ξ is a parameter, known as the "soil constant" [8], and ν is Poisson's ratio. Another useful representation for Equation (1), found by Weiskopf, is as follows: where χ ≥ 1 is a parameter known as "sandiness"; at χ = 1, the sandy medium becomes elastic and isotropic. Equation (1) and its natural generalization is discussed in more detail in the next section. The concept of constitutive Equation (1) proposed in [1] was applied to analyze (i) the soil pressure on retaining walls [5,6,9]; (ii) the propagation of the horizontally polarized SH and Love waves in a medium containing a sandy layer [10][11][12][13][14]; (iii) the propagation of the polarized, in the vertical (sagittal) plane, Rayleigh and Lamb waves [15][16][17]; and (iv) the propagation of the edge waves [18].

Head SP Waves
Head SP waves, discovered by Jeffreys (1926), play an important role in transmitting wave energy in the close vicinity of the epicentral zone of an earthquake or undersurface blast loading [8,[19][20][21]. SP waves are also the subject of extensive studies on the acoustics of layered composites [22,23]. In this respect, also see applications of energy dissipation for SP head waves along with Rayleigh, Rayleigh-Lamb, and other types of SAWs [24][25][26][27][28].
The considered head SP waves (bulk S wave being converted into P wave at the free surface) propagating along a free surface of a substrate arise at the reflection of bulk S waves from the free surface of a substrate at the critical angle. Figure 1 displays the application of Snell's law for the incident and reflected waves.
where the upper indices i and r refer to the incident and reflected waves, correspondingly, α i S = α r S ; α r P is the angle of the reflected P wave; and c S and c P are the corresponding velocities, c S = √ µ/ρ and c P = (λ + 2µ)/ρ; herein, λ and µ are Lame constants, and ρ is the material density. It is assumed that µ > 0 and µ > 3λ/2, which ensures positive definiteness of the elasticity tensor [19,24]. s 2023, 12, x FOR PEER REVIEW

Head SP Waves
Head SP waves, discovered by Jeffreys (1926), play an import wave energy in the close vicinity of the epicentral zone of an earth blast loading [8,[19][20][21]. SP waves are also the subject of extensive of layered composites [22,23]. In this respect, also see application for SP head waves along with Rayleigh, Rayleigh-Lamb, and other The considered head SP waves (bulk S wave being converted surface) propagating along a free surface of a substrate arise at waves from the free surface of a substrate at the critical angle. Fi plication of Snell's law for the incident and reflected waves. A scheme for the head SP wave origination: S1 wave propagat giving rise to the reflected SS1 and SP1 waves; S2 wave propagates at the c the reflected SS2 wave and SP2 head wave. Assuming α r P = 90 o yields the critical angle of incidence for the S wave: At the angle α i S,critical , the incident S wave gives rise to the P wave that travels along the free surface; such a wave is known as the SP head wave [19,20].
By applying ray methods [29,30], the critical distance from the epicenter, where the SP head wave starts to propagate, was found: where h is the depth of the buried pulse. The principal scheme, in a case accounting for the Earth's curvature, requires a modification; see Figure 2.
Thus, in order to account for the Earth's curvature, we require angles of the S wave incident that are slightly smaller than the critical angle defined by Equation (4). Some other problems related to finding the dynamic displacement field from the buried pulse loading are analyzed in [31][32][33][34]. While the curvature effects are important in geophysical large-scale models, the current research is confined by a simpler model with a plain boundary.  propagates at the subcritical angle, giving rise to the reflected SS1 wave and SP1 quasi head S2 wave propagates at the critical angle, giving rise to the reflected SS2 wave and SP2 tru wave.
Thus, in order to account for the Earth's curvature, we require angles of the incident that are slightly smaller than the critical angle defined by Equation (4) other problems related to finding the dynamic displacement field from the buried loading are analyzed in [31][32][33][34]. While the curvature effects are important in geop large-scale models, the current research is confined by a simpler model with boundary.
It should be noted that head waves have been intensively studied in quite number of works [19,20] related to seismic wave analyses of nearby earthquake physical prospecting, and explosion seismology over the past 100 years or so. In exploration, it is customary to call the SP waves the refracted or refractive waves, a method based on them is called the refractive method [20]. In this respect, Moho (1910) appears to have been the first to use head waves in earthquake records to in layering of the continental Earth's crust. While his analysis was based on geom optic considerations, thus limiting to high frequency waves, it was later revealed that his analysis remains valid in waves of an arbitrary frequency.
And the final remark concerns the validity of the ray methods, which are used in the discussed works [19,20,22]. It is assumed that (i) the velocity gradie much smaller than the circular frequency, which is not the case in most real situatio the radii of the curvature of the interfaces are significantly greater than the wave which is applicable to the nearby regions to the epicenter; and (iii) there is no pos to derive an exact formula relating the frequency and degree of accuracy of the res a given number of members of the ray series [19,20,22]. Actually, the accuracy formulas depends not only on the parameters of the media but also on the direction ray and other conditions. These shortcomings eminent to the ray method will b come in the next sections, where the FE method is applied to constructing the so for the analysis of SP wave arrivals at the points of observation located on a free su

Current Research
Herein, presumably for the first time, a comparative analysis of propagating head waves in isotropic and sandy media is performed with FE (finite element) mo coupled with constructing analytical expressions for the (apparent) velocities o waves. The applied FE models utilize solutions to the inner Lamb problem of the pulse, allowing for the generation of desired SP head waves along with P and S bu Rayleigh waves. In this respect, please see [35][36][37][38]. It should be noted that head waves have been intensively studied in quite a large number of works [19,20] related to seismic wave analyses of nearby earthquakes, geophysical prospecting, and explosion seismology over the past 100 years or so. In seismic exploration, it is customary to call the SP waves the refracted or refractive waves, and the method based on them is called the refractive method [20]. In this respect, Mohorovicic (1910) appears to have been the first to use head waves in earthquake records to infer the layering of the continental Earth's crust. While his analysis was based on geometrical optic considerations, thus limiting to high frequency waves, it was later revealed [19,20] that his analysis remains valid in waves of an arbitrary frequency.
And the final remark concerns the validity of the ray methods, which are mainly used in the discussed works [19,20,22]. It is assumed that (i) the velocity gradients are much smaller than the circular frequency, which is not the case in most real situations; (ii) the radii of the curvature of the interfaces are significantly greater than the wavelength, which is applicable to the nearby regions to the epicenter; and (iii) there is no possibility to derive an exact formula relating the frequency and degree of accuracy of the results for a given number of members of the ray series [19,20,22]. Actually, the accuracy of ray formulas depends not only on the parameters of the media but also on the direction of the ray and other conditions. These shortcomings eminent to the ray method will be overcome in the next sections, where the FE method is applied to constructing the solutions for the analysis of SP wave arrivals at the points of observation located on a free surface.

Current Research
Herein, presumably for the first time, a comparative analysis of propagating true SP head waves in isotropic and sandy media is performed with FE (finite element) modeling coupled with constructing analytical expressions for the (apparent) velocities of head waves. The applied FE models utilize solutions to the inner Lamb problem of the buried pulse, allowing for the generation of desired SP head waves along with P and S bulk and Rayleigh waves. In this respect, please see [35][36][37][38].
The analysis reveals a substantial discrepancy in both velocities of the propagating bulk and head waves along with magnitudes of the SP head waves, indicating (i) a slower apparent SP head wave velocity in sandy media than in the equivalent isotropic one; (ii) decreasing velocity values with an increasing sandiness parameter value; (iii) an increasing critical distance d SP with the increase in the sandiness parameter value; and (iv) an increase in the SP head wave magnitudes compared with the magnitudes of the same waves in the equivalent isotropic medium. The latter fact is of special importance in estimating SP wave magnitudes in the vicinity of the epicenters of shallow earthquakes and undersurface blast loadings.
Thus, the current research targeting the analysis of bulk and SP wave arrivals at a free surface of a Weiskopf halfspace under the action of an impulse delta-like force load, which is applied in the interior of the halfspace, becomes, apparently, the first study of such a problem. It should be emphasized that this research differs from a number of previous studies related to the analysis of the so-called inner Lamb problem of a point force load inside a homogeneous isotropic halfspace [36], and studies of surface acoustic wave propagation in layered media, including ones containing Weiskopf layers [39]. Another remark concerns the performed analysis of the displacement fields in two principally different Weiskopf media, differing by the Weiskopf parameter, and the latter specifies the so-called "sandiness" [1] and is of special importance in geophysics [5]. Again, as the literature review shows, the current research is presumably the first study related to the analysis of head waves appearing in Weiskopf media.

The Modified Weiskopf Model
The original Weiskopf theory [1] was applied to solve static problems for a sandy halfspace, satisfying the equations of state Equation (1). Considering these equations in the framework of anisotropic elasticity reveals that the Weiskopf medium is eventually an anisotropic one characterized by three independent elastic constants: E, ν, µ. Such a three-constant medium possesses cubic symmetry with the elasticity tensor C in the principal crystallographic axes, written as [40] Following [40], only the upper triangular part is shown in the expressions of Equation (6). Introducing three other material parameters, yields η > −2λ, λ, µ > 0 where relations in Equation (8) follow from the positive definite condition imposed on elasticity tensor Equation (6). Indeed, the positive definite condition for the elasticity tensor corresponding to an arbitrary elastic anisotropy reads as where A is an arbitrary non-vanishing second-order tensor, and C stands for the fourthorder elasticity tensor, which is symmetric with respect to its outer pairs of indices. Note also that the positive definite condition implies the strong ellipticity condition of the elasticity tensor; however, the reverse implication is not true. The latter expressions are convenient in comparing the components of the elasticity tensor for the considered cubic medium with ones for an isotropic medium, for which an additional relation holds: Now, the Weiskopf soil and sandy constants [1] in Equation (1), written in terms of parameters η, λ, and µ, become note that for an isotropic medium at η = λ + 2µ, Equation (11) yields ξ = 2(1 + ν) and χ = 1. Applying Equation (11), the elastic modulus η can be expressed in terms of other elastic moduli and the corresponding soil and sandy parameters: Equation (12) can be used for computing the elastic modulus η for the variation in the parameter χ.

Inner Plain Lamb Problem
To analyze head SP waves, the plane (strain) inner Lamb problem with a buried vertical pulse load was considered (Figure 3).
Applying Equation (11), the elastic modulus  can be expressed in terms of other elastic moduli and the corresponding soil and sandy parameters: Equation (12) can be used for computing the elastic modulus  for the variation in the parameter  .

Inner Plain Lamb Problem
To analyze head SP waves, the plane (strain) inner Lamb problem with a buried vertical pulse load was considered (Figure 3). According to the theory of head waves, in the first part, the wave traveling past the SP wave travels at the speed of S c ; in the second part, the wave travels at the speed of P c . Taking into account the expression in Equation (4)   According to the theory of head waves, in the first part, the wave traveling past the SP wave travels at the speed of c S ; in the second part, the wave travels at the speed of c P . Taking into account the expression in Equation (4) for the critical angle and an easy verified relation, tan(arcsinα) = α √ 1 − α 2 (13) the arrival time of the SP wave to a point of observation on the free surface becomes Now, Equation (14) allows us to define the apparent velocity of the SP head wave: Equation (15) shows that the apparent velocity of the SP head wave depends upon the relative depth of the origin and distance from the epicenter.

Boundary Conditions
For the FEA (finite element analysis), the following boundary conditions were imposed. Regarding edge A (Figure 3), the traction free boundary conditions are as follows: where t ν is the surface traction field and ν is the normal unit outward to the upper boundary ∂Ω A . Regarding edge B (Figure 3), symmetry boundary conditions are as follows: where a comma denotes the corresponding partial derivative. Along edges C and D (Figure 3), the non-reflecting (absorbing) boundary conditions based on the perfectly matched layers (PML) were specified; see [41,42]. The PML ensures efficient absorbing even at oblique angles of incidence [43,44].

Initial Conditions 22
The initial conditions at t = 0 were assumed as corresponding to the state of rest:

Force Loading
The applied force representing the buried delta pulse can be represented as a multiplication of the two delta functions: where p 0 is the force magnitude, n is the unit vector directed along the x 2 axis (Figure 3), and x 0 is the point on edge B where the force was applied.

Numerical Implementation
As was pointed out earlier, the governing equation of motion is solved by the FE method for spatial discretization and the FD (finite difference) numerical scheme for the time-domain discretization. These were implemented with the explicit Lax-Wendroff energy-preserving numerical scheme [45,46]. To achieve a conditionally stable numerical algorithm, the Courant-Friedrichs-Lewi condition was imposed on the time increment [46].
where ∆t CFL stands for the so-called Courant-Friedrichs-Lewi time increment, ensuring the numerical stability condition for the initially conditionally stable numerical scheme; ∆x min is the minimal linear size of the finite element used; and c P max is the maximal P wave velocity, if multiple materials are considered, or if the considered single material has varied physical properties. Thus, all the numerical computations were (mechanical) energy-preserving and numerically stable, which is extremely important for modeling and analyzing shock wave propagation [47,48]. It should also be noted that the hourglass control was used for increasing element stiffness and avoiding the appearance of nonphysical self-balanced stress fields, associated with the hourglass type deformation of the adjacent finite elements [45]. Another remark concerns the use of four-node linear quadrilateral elements with reduced integration, which are reasonably accurate in a FEA related to solving plane problems of non-stationary dynamics [49].
To minimize the non-physical oscillations behind the wave fronts, the time parameter T was chosen to satisfy the following empirical condition [36]: where ∆t Courant is the Courant number [50]. The magnitudes of the displacement fields, which arrived at the free surface, are shown in Figure 4 for a point located from the epicenter at 6h (h is the depth of the buried pulse). The plots characterizing the arrivals of the SP waves at other points of observation are analogous.
quadrilateral elements with reduced integration, which are reasonably accurate in a FEA related to solving plane problems of non-stationary dynamics [49].
To minimize the non-physical oscillations behind the wave fronts, the time parameter T was chosen to satisfy the following empirical condition [36]: where Courant t  is the Courant number [50]. The magnitudes of the displacement fields, which arrived at the free surface, are shown in Figure 4 for a point located from the epicenter at 6h (h is the depth of the buried pulse). The plots characterizing the arrivals of the SP waves at other points of observation are analogous.  Analyzing the peaks related to the arrival of P and SP waves in both sandy media shows that the peaks arrive at the same time and have almost identical magnitudes. Actually, the velocities of P waves traveling along the horizontal pass are exactly identical, while the velocities of S waves traveling along the inclined pass are substantially different due to a difference in shear moduli; this gives an explanation for the notable difference in the arrival times of the S wave signals. Another remark concerns the apparent lack of dispersion of SP waves in the considered media; see Figure 4. This phenomenon is also discussed in [8,[19][20][21].

Concluding Remarks
As the presented literature review shows, the problem of bulk wave propagation in a laminated anisotropic halfspace is analyzed apparently for the first time; previous studies are mainly related to the bulk wave propagation in either isotropic halfspaces, e.g., [36], or to the dispersive SAW (surface acoustic waves) in anisotropic layered media [39,[51][52][53][54]. Analyzing the peaks related to the arrival of P and SP waves in both sandy media shows that the peaks arrive at the same time and have almost identical magnitudes. Actually, the velocities of P waves traveling along the horizontal pass are exactly identical, while the velocities of S waves traveling along the inclined pass are substantially different due to a difference in shear moduli; this gives an explanation for the notable difference in the arrival times of the S wave signals. Another remark concerns the apparent lack of dispersion of SP waves in the considered media; see Figure 4. This phenomenon is also discussed in [8,[19][20][21].

Concluding Remarks
As the presented literature review shows, the problem of bulk wave propagation in a laminated anisotropic halfspace is analyzed apparently for the first time; previous studies are mainly related to the bulk wave propagation in either isotropic halfspaces, e.g., [36], or to the dispersive SAW (surface acoustic waves) in anisotropic layered media [39,[51][52][53][54].
The performed analysis reveals the timings and amplitudes of both bulk and SP head waves, which either arrive at the free surface (bulk P and S waves) or travel along the free surface (SP head waves). One of the most important findings is associated with a noticeable discrepancy between the arrival times of the considered waves propagating in the Weiskopf media with different sandiness parameters. Moreover, it is found that the wave amplitudes are also prone to the sandiness parameter value; it is actually observed that in a vicinity of the SP wave arrivals, the wave amplitudes differ enough to distinguish one type of the Weiskopf material from another.
Another remark concerns the adopted FE modeling techniques, which are based on the explicit Lax-Wendroff energy-preserving numerical scheme [45,46]. To achieve a conditionally stable numerical algorithm, the Courant-Friedrichs-Lewi condition was imposed on the time increments. Thus, all the numerical computations were (mechanical) energy-preserving and numerically stable, which is extremely important for modeling and analyzing shock wave propagation [47,48].
In works on the analysis of head waves [55][56][57], Snell's law is quite often used to obtain the so-called critical angle, which is necessary for the appearance of SP head waves. In this respect, it was found by both theoretical and experimental studies that Snell's law is valid for not only geometrical optic approximations, when the elastic wavelength is (infinitely) small compared with the distances from the interface [58], but this law remains applicable in cases of low frequencies, when the wavelength becomes sufficiently large [19].