1 Introduction

The theoretical analysis of seismic waves in multi-layered elastic media is of great practical importance due to its varied and wide applications in many engineering and applied sciences, namely, seismology, geophysics, rock mechanics, soil mechanics, petroleum engineering, structural engineering, and hydrogeology. Numerous works are available on the experiments and theories functioned to study the seismic wave regulation phenomena such as radiation, propagation, reflection, and refraction through deformable solids. The responses of seismic wave propagation in some complex media are very useful in the exploration of natural resources buried inside the Earth's surface, e.g., oils, gases, other useful hydrocarbons and minerals. The books of renowned authors such as Love ^[1], Ewing et al. ^[2], Biot ^[3], and Gubbins ^[4] contain relevant basic information of seismic waves.

In the earlier era, a number of authors investigated seismic waves in perfectly elastic media. However, the Earth cannot be dealt with exactly in this manner. From the second law of thermodynamics, the dissipation of energy can be found in all physical deformable solids in the context of internal friction mechanism. This loss of energy is responsible for the attenuation of seismic waves in imperfect elastic bodies ^[5]. The deficiency of perfect elasticity also decreases the wave amplitude because of anelastic attenuation. Imperfect elastic bodies can be considered to have properties intermediate between those of elastic and viscous bodies, and they are called viscoelastic bodies. Materials, such as coal tar, sediments, and salt, are examples of viscoelastic materials. In the asthenosphere, the materials are highly viscous, mechanically weak, and ductility deforming. The rocks in the lithosphere are materially cooler and behave like viscoelastic materials ^[6]. The most dynamical Earth processes occur in these zones and are accountable for the earthquake. Due to these facts, the problems of seismic waves and vibrations in viscoelastic media are receiving more attention from many researchers ^[7-17]. Červený ^[18] reviewed the inhomogeneous harmonic plane wave propagation through anisotropic viscoelastic media. Manolis and Shaw ^[19] derived the existence of harmonic waves in mild stochasticity exhibiting heterogeneous viscoelastic media. Recently, Singh et al. ^[20] modelled the propagation of Stoneley waves through the common interface of two different homogeneous viscoelastic semi-infinite media.

The real Earth materials such as soil or sand are found at each level and made of loosely connected mesoparticle grains or platelets. Weiskopf ^[21] presented a mathematical analysis for stresses in a semi-infinite substance, in which he introduced a constant for soils. He pointed out that due to mutual slipping of the granular particles, the shear modulus or rigidity in the sandy medium is much smaller in comparison with that in an ideal elastic solid material. Taking into consideration the importance of dry sandy materials in the composition of Earth, a number of researchers ^[22-25] studied the torsional wave propagation in media considering different aspects of the problem. Reflection and refraction of shear elastic waves at a plane interface between two sandy half-spaces in welded contact were investigated by Pal et al. ^[26]. Tomar and Kaur ^[27] developed a problem of reflection and transmission on the incident of SH-type plane wave at a cyclic corrugated interface of dry sandy and anisotropic elastic media. Some recent reviews on the propagation of torsional wave in the dry sandy medium as a Gibson half-space are presented by Shekhar and Parvez ^[28] and Vishwakarma et al. ^[29].

The wide variation and rapid changes in material properties of the Earth with the depth are typically due to the heterogeneity which affects the propagation of seismic surface waves significantly. Researchers and seismologists mostly favour the heterogeneous coupled field structure to analyze the underground response of seismic surface waves. Some exemplary works on heterogeneous media were acknowledged by several authors including Bullen ^[30] and Birch ^[31]. Dey et al. ^[32] accomplished the possibility of torsional wave propagation in linear, quadratic and exponential functioned heterogeneous elastic media. Gupta et al. ^[33] considered linear heterogeneity in a crustal layer for the study of torsional surface wave. Ke et al. ^[34-35] contemplated the Love wave propagation in linear and exponential heterogeneous fluid saturated porous media. Subsequently, several researchers ^[36-39] investigated seismic surface waves in various heterogeneous media.

We consider the torsional surface wave propagation in an initially stressed heterogeneous viscoelastic layer of finite thickness sandwiched between a layer of finite thickness and a half-space of heterogeneous dry sandy media. The heterogeneities are the functions of depth, and the initial stress arising in the viscoelastic layer is due to the horizontal compression and horizontal tension. Besides these, the bonded viscoelastic layer possesses both elastic and viscous (linear anelastic) properties which may be modelled as an infinite number of possible configurations of elastic springs and viscous dashpots (Kelvin-Voigt model). We obtain a complex frequency equation using the suitable boundary conditions for the present model. Both equations (real and imaginary parts of the complex frequency equation) consist of the wave number and the velocity of the wave along with other parameters. The real part of this frequency equation describes the dispersion phenomena, whereas the imaginary part describes damping phenomena of the wave. Hence, we may conclude that the velocity associated with the dispersion stands as the 'phase velocity' of wave, and the velocity associated with the damping may be termed as the 'damped velocity'. The affected behaviour of phase velocity (V_p) and damped velocity (V_d) with respect to the real wave number (k₁H₁) on the propagation of torsional wave is discussed in this study for various parameters, namely, the dissipation factor, the attenuation coefficient, the sandy parameters, the initial stress, the heterogeneity parameters, and the thickness ratio parameter. Mathematical equations reveal that the damping nature of the surface wave found in the present study is exclusively due to the viscoelastic coefficient (dissipation factor/viscosity) of the medium. It simply means that the damping of wave occurs, as long as the medium is viscoelastic. The damping phenomenon may be observed in the form of energy loss during viscous lubrication between moving particles of the viscoelastic medium.

2 Formulation and basic assumptions of problem

In the cylindrical coordinate system (r, θ, z), let us assume a geometry of composite structure in which a heterogeneous dry sandy half-space M₃ is overloaded by two layers, namely, an initially stressed viscoelastic layer M₂ and a heterogeneous dry sandy layer M₁, as shown in Fig. 1. The thickness of uppermost layer M₁ is H₂, and that of sandwiched layer M₂ is H₁. The origin of the cylindrical coordinate system (r, θ, z) is located at the interface separating the layer M₂ from the half-space M₃, and the z-axis is pointed positively downwards in the half-space. The direction of the torsional wave propagation is along the radial and causing displacements in the azimuthal direction only, hence the rate of change along the azimuthal direction vanishes, i.e., .

We consider transformations for the rigidity and density of the uppermost dry sandy layer,

(1)

for the sandwiched viscoelastic layer,

(2)

and for the dry sandy half-space,

(3)

where η₁>1 and η₂>1 are sandy parameters of uppermost layer and lower half-space, respectively, and υ₀ is the viscosity of sandwiched viscoelastic layer. The heterogeneities α >0, β >0, and γ >0 of the uppermost layer, sandwiched layer, and lower half-space, respectively, have dimensions equal to the length. If η₁=1 and η₂=1, then the uppermost layer and lower half-space correspond to elastic solid media. If υ₀=0, then the sandwiched viscoelastic layer also corresponds to elastic solid media.

Let (u_r, u_θ, u_z), (v_r, v_θ, v_z), and (w_r, w_θ, w_z) be vectors of displacement components of the uppermost layer, sandwiched layer, and lower half-space, respectively. Then, by the characteristic of torsional waves,

(4)

Torsional waves are horizontally polarized surface waves, giving the clockwise and anticlockwise twist to the medium about the direction of motion of wave. Therefore, to find the solution to displacement components for wave changing harmonically in the time along the radial direction, we assume solutions like the following equation in all the three media:

(5)

where k is the wave number, ω is the angular frequency, and J₁(kr) is the first order Bessel's function of the first kind.

2.1 Dynamics of uppermost layer M₁

In the absence of any body force, the only non-vanishing equation of motion for propagation of torsional wave in the dry sandy layer is given by ^[3]

(6)

where τ_rθ and τ_θz are stress components,

(7)

Using Eq. (7), the equation of motion (6) can be written as

(8)

Taking as the solution to Eq. (8) and substituting it into Eq. (8) with Eq. (1), we have

(9)

Taking substitution into Eq. (9), we obtain

(10)

where .

Therefore,

(11)

The displacement component u_θ of uppermost layer can be written as

(12)

2.2 Dynamics of sandwiched layer M₂

The only non-vanishing equation of motion for torsional wave propagation in the viscoelastic layer under the initial stress is given by ^[3]

(13)

where the stress components S_rθ and S_θz are

(14)

From Eqs. (13) and (14), we obtain

(15)

Substituting and Eq. (2) into Eq. (15), we get

(16)

where μ₂=μ₀₂+ iωυ₀.

By using substitution , Eq. (16) takes the form of

(17)

Here, .

Therefore,

(18)

Hence, the displacement component v_θ of viscoelastic layer is

(19)

2.3 Dynamics of lower half-space M₃

The only non-vanishing dynamical equation for torsional wave propagation in the dry sandy half-space without body force is given by ^[3]

(20)

where σ_rθ and σ_θz are stress components,

(21)

Substituting Eq. (21) into Eq. (20) yields

(22)

Let us assume as a solution to the above equation. Therefore, Eq. (22) with Eq. (3) provides

(23)

Substitute into Eq. (23), and f₃(z) satisfies the following differential equation:

(24)

where .

Therefore,

(25)

The appropriate solution for the lower half-space is

(26)

3 Boundary conditions

The following boundary conditions for the proposed model must be satisfied:

(ⅰ) At z=0, the displacement and stress components of the sandwiched layer and half-space are continuous, respectively,

(27)

(28)

(ⅱ) At z=-H₁, the displacement and stress components of the uppermost and sandwiched layers are continuous, respectively,

(29)

(30)

(ⅲ) At z=-(H₁+H₂), the uppermost layer is stress free,

(31)

4 Dispersion and damping equations

Substituting the solutions given in Eqs. (12), (19), and (26) into the boundary conditions, we get the following homogeneous algebraic system of equations for A₁, A₂, B₁, B₂, and D₁:

(32)

(33)

(34)

(35)

(36)

Eliminating constants A₁, A₂, B₁, B₂, and D₁ from the above equations, we obtain a closed form expression for a torsional surface wave,

(37)

where , and .

The above expression of gives the wave velocity profile of torsional waves in a heterogeneous viscoelastic layer bonded between a layer and a half-space of heterogeneous dry sandy media. Equation (37) is complex because of the dissipation of viscoelasticity. In fact, the wave number k is complex and thus may be written as

(38)

where k₁ and k₂ are real, and is the attenuation coefficient which is dimensionless. Thus, the phase velocity c of torsional wave can be evaluated by the relationship,

(39)

The dimensionless function is called the quality factor, and the inverse function of the quality factor is

(40)

which is known as the dissipation factor. The dissipation factor Q₀^-1 also represents the tangent of the angle by which the strain lags the stress in the viscoelastic layer.

To separate the real and imaginary parts of complex frequency equation (37), we assume

(41)

where r₁, r₂, r₃, x₁, x₂, x₃, y₁, y₂, y₃, θ₁, θ₂, and θ₃ are defined in Appendix A. Therefore, the separation of velocity equation (37) into real and imaginary parts yields the two following equations,

(42)

(43)

where E₁=x₂(T₁-Q₀^-1T₂) -y₂(Q₀^-1T₁+T₂), E₂=x₂(Q₀^-1T₁+T₂) +y₂(T₁-Q₀^-1T₂), and the values of T₁, T₂, T₃, and T₄ are defined in Appendix B.

The real part of , i.e., Eq. (42) deals with the dispersion phenomenon of torsional wave propagation and is known as the dispersion equation which originates the dispersion curves, which is the phase velocity against the real wave number, whereas the imaginary part of , i.e., Eq. (43) provides the damping characteristics of torsional wave propagation and is known as the damping equation which originates the damping curves, which is the damped velocity against the real wave number.

5 Particular case

If the uppermost layer is omitted, and both the intermediate layer and half-space are considered as stress free and homogeneous elastic, then the dispersion equation (42) reduces to the classical equation of Love wave ^[1],

(44)

and the damping equation (43) vanishes identically.

6 Numerical computation and discussion

To appraise the influenced characteristics of torsional wave propagation due to different parameters like heterogeneity parameters 1/(αk₁), 1/(βk₁), and 1/(γk₁), the dissipation factor Q₀^-1, the attenuation coefficient Ω, sandy parameters η₁ and η₂, the initial stress parameter ω =P₀/(2μ₀₂), and the thickness ratio of layers H =H₂/H₁ involved in the expression , we need numerical examples which best explain the theoretical study graphically. In this order, we take numerical data μ₀₁=3.23× 10¹⁰ N/m² and ρ₀₁=2 802 kg/m³ for the uppermost dry sandy layer, μ₀₂=7.1× 10¹⁰ N/m² and ρ₀₂=3 321 kg/m³ for the sandwiched viscoelastic layer, and μ₀₃=29.17× 10¹⁰ N/m² and ρ₀₃=5 563 kg/m³ for the dry sandy lower half-space from Gubbins ^[4]. All the numerical computations in this section are done by the dimensionless process. The fixed values of dimensionless parameters are provided in Table 1. The effects of all affecting dimensionless parameters on the dimensionless phase velocity c/c₂=V_p and the damped velocity c/c₂=V_d with respect to the dimensionless real wave number k₁H₁ of a torsional wave are shown in Figs. 2-10. An overview of all the figures establishes that the phase velocity of torsional wave always decreases with the increase in the real wave number, while the damped velocity increases for the very short frequency range of the real wave number and then falls down.

The variational effects of heterogeneity associated with the uppermost layer on the phase velocity and damped velocity of torsional wave are manifested in Figs. 2(a) and 2(b), respectively. Curve 1 of Fig. 2 shows the absence of heterogeneity in the uppermost layer (1/(αk₁) =0). It is clear from Fig. 2 that, as the magnitude of heterogeneity associated with the uppermost layer increases, the phase velocity also increases, while the damped velocity decreases in the frequency region 0.1 < k₁H₁ < 0.7 and increases after a point, nearly at k₁H₁=0.7. The impact of heterogeneity associated with the uppermost layer on the phase velocity is slightly significant at the higher frequency region compared with the lower frequency region. It has almost negligible effects on the damped velocity at the two point, nearly k₁H₁=0.1 and k₁H₁=0.7. The figure also suggests that the phase velocity of torsional wave is larger in the heterogeneous uppermost layer than the homogeneous uppermost layer, and the damped velocity is larger in the homogeneous uppermost layer than the heterogeneous uppermost layer for the frequency region k₁H₁ < 0.7. However, for the frequency region k₁H₁>0.7, it is reversible.

In Figs. 3(a) and 3 (b), the effects of heterogeneity associated with the bonded layer on the phase velocity and damped velocity are reflected, respectively. Curve 1 of Fig. 3 displays the absence of heterogeneity in the bonded layer (1/(βk₁) =0). The meticulous observation of figures concludes that, as the magnitude of heterogeneity associated with the bonded layer increases, the phase velocity decreases, whereas it increases the damped velocity. It can also be seen that the effect of heterogeneity associated with the bonded layer on the phase velocity is almost negligible at the higher frequency region compared with the lower frequency region, whereas it has a very significant effect on the damped velocity at the lower frequency region compared with the higher frequency region. It is very clear from the figure that the phase velocity of torsional wave is larger in the bonded homogeneous layer than the bonded heterogeneous layer. Moreover, the damped velocity is larger in the bonded heterogeneous layer than the bonded homogeneous layer.

Figures 4(a) and 4(b) portray the effects of heterogeneity associated with the lower half-space on the phase velocity and the damped velocity, respectively. Curve 1 of Fig. 4 illuminates the absence of heterogeneity in the half-space (1/(γk₁) =0). Minute observation of the figure depicts that the phase velocity and damped velocity of torsional wave increase as the magnitude of heterogeneity associated with the lower half-space increases. It is also established through the figure that the phase and damped velocities of torsional wave are larger in the heterogeneous half-space than the homogeneous half-space. Furthermore, the phase velocity is affected more significantly by heterogeneity associated with the uppermost layer compared with other heterogeneities of the medium, and the damped velocity is affected more significantly by heterogeneity associated with the bonded layer compared with other heterogeneities of the medium.

The effects of dissipation factor associated with the bonded layer on the phase velocity and damped velocity of torsional wave are shown in Figs. 5(a) and 5(b), respectively. Curve 1 in Fig. 5 indicates the case when the bonded layer is viscous free (Q₀^-1=0). It can be observed from Fig. 5(a) that, the phase velocity also increases with the increase in the magnitude of dissipation factor, and the increasing effect becomes more significant for the larger magnitude of dissipation factor as well as for the higher frequency region compared with the lower frequency region. Figure 5(b) reflects that, as the magnitude of dissipation factor increases, the damped velocity increases for the frequency region k₁H₁ < 0.6 and decreases for the frequency region k₁H₁>0.9. As the magnitude of dissipation factor decreases, the damping nature of damping curves also decreases. In the presence of viscosity in the viscoelastic layer, the damping nature of damping curves is more than that of their absence. Hence, it can be concluded that the viscosity of viscoelastic layer is responsible for the damping nature of damped velocity curves.

Figures 6(a) and 6(b) describe the effects of attenuation coefficient arising due to the complex wave number on the phase velocity and damped velocity of torsional wave, respectively. It is clear from Fig. 6 that the phase velocity increases and the damped velocity decreases with an increase in the magnitude of attenuation coefficient. The effect of attenuation coefficient on the phase velocity is found negligible for the lower frequency region compared with the higher frequency region. It is also noticed from the figure that, in the case of complex wave number (δ≠0), the phase velocity is larger than the case of real wave number (δ = 0), and the damped nature in the damped velocity curve is more for the real wave number than the complex wave number.

Figures 7(a) and 7(b) reveal the effects of sandiness associated with the uppermost layer on the phase velocity and damped velocity of torsional wave, respectively. Figures 8(a) and 8(b) show the effects of sandiness associated with the half-space on the phase velocity and damped velocity of torsional wave, respectively. Curves 1 of Figs. 7(a) and 7(b) represent the elastic uppermost layer (η₁=1), and Curves 1 of Figs. 8(a) and 8(b) represent the elastic lower half-space (η₂=1). The meticulous inspection of Figs. 7(a) and 8(a) delineates that the phase velocity of torsional wave increases with the increase in the sandiness associated with the uppermost layer or the half-space. However, Figs. 7(b) and 8(b) suggest that, the damped velocity of torsional wave increases with the increase in the sandiness associated with the half-space, and it decreases with the increase in the sandiness associated with the uppermost layer for the frequency region k₁H₁ < 1.5. It is also very clear from Figs. 7 and 8 that the increasing effect of sandiness in the half-space is more prominent than the sandiness in the uppermost layer on the phase velocity as well as the damped velocity of the wave. However, for a particular real wave number, the phase and damped velocities are found larger for a particular magnitude of sandiness associated with the uppermost layer than the sandiness associated with the half-space. Moreover, the phase velocity of the wave is found larger in the sandy layer and sandy half-space than the elastic layer and elastic half-space, respectively. However, the damped velocity is found larger only in the elastic upper layer than the sandy upper layer, while for the half-space, it is reversible.

The curves plotted in Figs. 9(a) and 9(b) elucidate the effects of initial stress (acting on the bonded layer) on the phase velocity and damped velocity of torsional wave, respectively. In particular, Curves 1 and 2 represent that the bonded layer is under the tensile initial stress (negative magnitude of ω), Curve 3 represents a stress free medium (ω =0), and Curves 4 and 5 show that the medium is under the compressive initial stress (positive magnitude of ω). It is clear from Fig. 9 that the phase velocity and damped velocity decrease with an increase in the magnitude of initial stress. The effect of initial stress is significant on the phase velocity as well as the damped velocity for the lower frequency region compared with the higher frequency region. More expressly, the phase velocity and damped velocity increase with an increase in the magnitude of tensile initial stress, whereas the phase velocity and damped velocity decrease with an increase in the compressive initial stress acting on the bonded viscoelastic layer.

In Figs. 10(a) and 10(b), the impacts of thickness ratio of uppermost layer and bonded layer are elucidated on the phase velocity and damped velocity of torsional wave, respectively. It has been found from Fig. 10 that, the phase velocity decreases patently as the value of thickness ratio increases, however, the damped velocity increases for the lower frequency range of k₁H₁ < 0.7 and is found to be decreased for the higher frequency region of k₁H₁>0.7. At a point, nearly k₁H₁=0.7, the damped velocity is almost negligibly affected by the magnitude of thickness ratio.

7 Conclusions

Within the framework of a heterogeneous layered Earth model of viscoelastic and sandy materials, an analytical study has been carried out for the torsional wave propagation. A closed form velocity profile equation is obtained for this model and separated into two velocity equations, namely, the dispersion and damped equations. Both the velocity equations are in well agreement with the classical Love wave condition. The conspicuous discussion on the affected nature of phase and damped velocities due to the heterogeneities, viscosity, sandiness, external stress, layers thickness and complex wave number is collected in the numerical and discussion section. Some salient consequences of the present investigation are given below.

(ⅰ) The heterogeneities in the said model have considerable effects on both phase and damped velocities. The presence of heterogeneity in the uppermost layer or in the half-space enhances the phase velocity of the wave, whereas in the bonded layer, it diminishes the phase velocity of the wave. On the other hand, the damped velocity of torsional wave increases in the presence of heterogeneity in the bonded layer or in the half-space; whereas its presence in the uppermost layer enhances the damped velocity for the higher frequency region and diminishes the damped velocity for the lower frequency region.

(ⅱ) The dissipation factor plays a vital role in the propagation of torsional wave and is largely responsible for the damping nature of damped velocity of the wave. It enhances both the phase velocity and damping velocity of the torsional wave.

(ⅲ) The complexity of the wave number (i.e., the attenuation factor arises due to the complex wave number) enhances the phase velocity and diminishes the damped velocity.

(ⅳ) The presence of sandiness in the uppermost layer enhances the phase velocity as well as the damped velocity of wave. However, in the half-space, it enhances the phase velocity and diminishes the damped velocity of torsional wave.

(ⅴ) The tensile initial stress presence in the bonded layer enhances both phase and damped velocities, while the presence of compressive initial stress diminishes both phase and damped velocities.

(ⅵ) As we increase the thickness of the uppermost layer or decrease the thickness of the bonded layer, the phase velocity diminishes, while the damped velocity enhances for the lower frequency region and diminishes for the higher frequency region.

It is very clear from the study that the presences of heterogeneities, viscosity, sandiness, external stress, layers thickness, and complex wave number affect the propagation of torsional wave significantly. It is expected that the results of the present study may be useful in the further theoretical or observational studies of torsional or twisting-type surface waves propagation in the more realistic models of the viscoelastic and sandy materials in the Earth. Furthermore, due to possible applications in geophysics, the obtained results are very useful in many fields of applied sciences, for example, in the analysis of seismic data, soil mechanics, composite materials, civil engineering, and earthquake science.

Appendix A Appendix B Acknowledgements The authors would like to convey their deep sense of indebtedness to Indian Institute of Technology (Indian School of Mines) Dhanbad for providing all necessary resources for this work.