Love-Type Waves in Multilayered Elastic Media Containing Voids: Haskell Matrix Method

Propagation of Love-type waves is investigated in a model of multilayered elastic solid half-space containing voids. Using Haskell matrix method, the dispersion relations are derived for two situations of topmost boundary surface of the model: (a) when it is stress-free and (b) when it is rigid. For both the situations, there exist two fronts of Love-type waves that are dispersive in nature and propagating with distinct speeds. One of these wave fronts is analogous to that obtained for multilayered Cauchy elastic half-space, while the other front is new and appeared due to the presence of voids in the model. Conditions of propagation for both the fronts of Love-type waves are derived analytically for 2-layered model. Comparison of dispersion curves corresponding to both these fronts is also depicted graphically for 2- and 3-layered models with stress-free boundary surface. Effect of layer thickness and presence of void parameters is studied on the speeds of wave fronts and depicted graphically for a specific 2-layered model.

for various purposes like detection of viruses and bacteria, DNA characterization, etc. Thus, it can be said that a same physical phenomenon can sometimes be harmful and at other times it can be useful like Love waves. Love waves are SH -waves, which can propagate in a medium that contains a thin elastic layer overlying an elastic half-space (see Love [13]). Although, there is vast literature available on the propagation of Love waves, but still they are of utmost interest due to their important applications in diverse fields. As a matter of reference, some notable works pertaining to Love waves are contained in Dai and Kuang [3], Dey et al. [4], Kaur et al. [10], Liu et al. [12], Singh [21] and many more.
An elastic material containing small vacuous pores distributed evenly throughout is known as an elastic material with voids. The development of non-linear and linear theories of elastic material containing voids was made by Nunziato and Cowin [15] and Cowin and Nunziato [2] based on the idea of distributed bodies stemmed from Goodman and Cowin [7]. This theory rests on the concept of a material for which (i) the bulk density is written as product of the density field of the matrix material and the void volume fractional field, (ii) the reference configuration is assumed to be free from strains and stresses. Due to such representation of bulk density of the material, an additional degree of freedom called change in void volume fraction field is introduced in the theory. Hence, there are four independent kinematical variables including three of displacement field. The transmission of loads across a surface element is characterized not only by the force stress tensor, but also by the generalized force known as equilibrated stress. The theory of elastic material with voids is thus an extension of the classical theory of elasticity. The behavior of plane harmonic waves in a linear elastic solid containing voids was first investigated by Puri and Cowin [17], who showed that there exist two dilatational waves consisting of a predominantly dilatational wave of classical elasticity and other is predominantly a wave carrying the change in void volume fraction. Various problems of waves and vibrations in different types of elastic materials containing voids have been attempted by several researchers in the past. Some notable works among them are given by Chandrasekharaiah [1], Singh et al. [19], Singh and Tomar [20], Tomar [23], Tomar and Ogden [24], and Wright [25] including several others.
There are various interesting problems in the field of seismology or elsewhere that involve the propagation of disturbances in layered elastic half-space. For example, Earth's crust is made up of several elastic layers, each layer exhibiting different material properties and resting on the upper mantle. It is doubtless to say that the material inside Earth's crust may be porous and it may be best simulated with elastic material containing voids. This has motivated us to investigate the propagation of surface waves in multilayered elastic media containing voids. Recently, authors have studied the propagation of Rayleigh-like surface waves through the multilayered elastic media containing voids using Haskell matrix method. It has been shown that the presence of voids in the model affects the speed of Rayleigh-like wave significantly for small values of wavenumber (see Khurana et al. [11]). In the fundamental papers by Thomson [22] and Haskell [8], they have developed a matrix method in order to derive the dispersion equations for Rayleigh and Love waves traveling through multilayered elastic media. In the literature, this method is called as Thomson-Haskell matrix method, and it has been found that it is very useful in deriving the dispersion equations efficiently which was otherwise formidable for models consisting of even more than two layers. Some computational programs based on this matrix formulation have been reported in literature that conveniently plot the dispersion curves for the models consisting of considerable number of elastic layers (see Rosenbaum [18]). An alternate matrix formulation was introduced by Dunkin [5] to handle the numerical difficulties that arise while solving the dispersion equations for multilayered elastic models at high frequencies. Recently, using finite difference scheme, Panja and Mandal [16] studied the propagation of Love waves in multilayered viscoelastic orthotropic medium and investigated the effect of viscoelastic parameters and initial stress on the speed of Love waves.
In the present problem, we have discussed the dispersion of Love-type surface waves in a multilayered isotropic linear elastic solid containing uniformly distributed voids. These voids are vacuous pores and contain nothing of energetic or physical significance. The model comprises of n layers (where the last layer is assumed to be a half-space) of different thicknesses resting on each other in welded contact. We consider two different situations of the model depending on the nature of topmost boundary surface, namely, (i) mechanically stress-free and (ii) rigid. For a two-dimensional problem, a standard procedure to derive dispersion relation is as follows. First appropriate boundary conditions are set up at each interface as well as at the topmost boundary that are the continuity of displacement, change in void volume fraction, shear stress and equilibrated stress at each interface, while vanishing of both the stresses (for stress-free boundary surface) or both the kinematic variables (for rigid boundary) at the boundary surface. In each case, there are 4n − 2 boundary conditions to be satisfied, which yield a homogeneous system of 4n − 2 equations in same number of unknown constants, to be determined. The condition for the existence of non-trivial solution of this system is the vanishing of determinant of the corresponding coefficient matrix. The determinental equation so obtained is the dispersion relation, whose roots provide us the speed(s) of propagating waves and other characteristics inherited. But this procedure of finding roots is formidable even for more than two layers. To overcome this problem, Haskell [8] proposed a matrix method which is described as below. Firstly, a linear relation is developed between the requisite physical quantities at the top and bottom of each layer. By imposing boundary conditions, a relationship is then developed between the physical quantities at the topmost boundary surface with those at the last interface. Finally, the required dispersion relations for both the models are obtained by imposing the radiation condition for the propagation of surface waves. It has been concluded that two fronts of Love-type surface waves exist that are dispersive in nature and propagating with distinct speeds. One of these fronts of Love-type waves is independent of the presence of void parameters in the model and analogous to Love-type wave propagating through multilayered elastic half-space (particularly for stress-free boundary surface). However, the other front is new and appeared due to the presence of voids in the model. Further, these dispersion relations are then reduced for 2-layered and 3-layered models. A particular numerical data has been borrowed from the literature to perform computations when the boundary of topmost layer is taken to be stress-free. Dispersion curves of both the wave fronts are depicted and compared for 2-layered and 3-layered models. In particular, for n = 2, the effects of presence of voids and the thickness of layer on the corresponding wave speeds have been studied. The boundary conditions imposed in the 2-layered model have also been verified graphically. Further, the conditions of propagation for both the fronts of Love-type waves are derived analytically for the 2-layered model and the results have been verified numerically.

Basic Equations
Following Cowin and Nunziato [2], the constitutive relations and equations of motion without body force densities for a linear isotropic elastic material with voids are given by and where τ ij is the stress tensor, h i is the equilibrated stress vector, g is the intrinsic equilibrated body force, u (= u i ) denotes the displacement vector, φ is the change in void volume fraction, e ij is the strain tensor, δ ij denotes the Kronecker delta, and i, j ∈ {x, y, z}. The well known relation between strain tensor and displacement vector is given by The material parameters λ and μ are the well-known Lame's constants, whereas α, β, ξ , ω are the void parameters. A comma in the subscript denotes the partial spatial derivative, while the superposed dot denotes the partial temporal derivative. Here, x, y, z are the spatial variables, t is the time variable, ρ denotes the bulk density, and χ represents the equilibrated inertia. In addition to the above, these material moduli must satisfy the following inequalities for the strain energy density function to be positive definite.
It can be noticed that the term, namely, ωφ occurring in (1) and (2) is due to mechanical dissipation of the material. Material exhibiting such phenomenon is known as Kelvin-Voigt or simply Voigt material in the literature. Viscoelastic material or material possessing slow rate deformations lies under the category of Voigt material. In the present work, we shall consider the material possessing no mechanical dissipation, that is, non-Voigt material. Hence, this term will be ignored from equations (1) and (2) for further formulation (see also Iesan [9] and Singh and Tomar [20]).

The Model and Problem
We desire a model of n − 1 layers having finite thicknesses consisting of elastic medium with voids lying over an elastic solid half-space containing voids. Here, the model will be referred to as n-layered elastic solid half-space containing voids. With reference to a rectangular Cartesian co-ordinate system, we consider that each layer consisting of uniform elastic solid having uniform thickness with evenly distributed small voids is parallel to the x − y plane. The interfaces are numbered from (0) th to (n − 1) th , viz, the (m) th interface represents the interface between m th and (m + 1) th layers. We assume that the layers are in welded contact at each interface. Also, the z-axis points vertically downward into the medium with the plane z = 0 coincident with the (m − 1) th interface as shown in Fig. 1. Throughout the formulation, the subscript 'm' corresponds to the physical quantity in the m th layer (1 ≤ m ≤ n). We consider a two-dimensional problem in x − z plane for the propagation of Love-type waves, therefore, we shall take where (u m , v m , w m ) are the components of displacement vector u m and φ m denotes the change in the void volume fraction in the m th layer. With these considerations, the requisite Here, (Y z ) m denotes the transverse component of shear stress, that is, the y-component of the stress vector acting across the plane to which z-axis is normal, (h z ) m denotes the normal component of equilibrated stress vector, and (g) m is the intrinsic equilibrated body force. Using (3) into the field equation (2) 1 , we obtain the following scalar equations From (3), we note that φ m is considered to be the function of spatial co-ordinates x and z, therefore, first and last equations in (5) imply that the coefficient β m vanishes. It can be noted from equation (2), that there is no coupling between the displacement and the change in void volume fraction field, hence the corresponding effects will be independent. Thus, the field equations given in (2) reduce to the following uncoupled equations For the propagation of Love-type surface waves in the positive direction of x-axis with speed c, the time-harmonic plane wave solutions of above equations are taken in the form given as where f m and g m are functions of spatial co-ordinate z; k is the wavenumber connected to the angular frequency p through the speed c by the relation p = kc. And the angular frequency p is related to the linear frequency f through the relation p = 2πf .
Inserting the above expressions of v m and φ m into equations given in (6), we obtain where and prime ( ) denotes the derivative with respect to z. Here, b m corresponds to the speed of classical distortional wave and γ m denotes the speed of wave associated with change in void volume fraction. The general solutions of equations in (6) are given by where A m , B m , C m , and D m are arbitrary constants. Using (10) and (11) into the expressions of requisite stress components (Y z ) m and (h z ) m given in (4), we obtain The quantities given in (10)- (13) can be written as follows

Boundary Conditions
We have assumed that any two adjacent layers of the model are in welded contact with each other, therefore, the boundary conditions at each interface will be the continuity of following physical quantities: (i) the transverse component of displacement, (ii) the change in void volume fraction, (iii) the transverse component of shear stress, and (iv) the normal component of equilibrated stress. We shall derive a relationship between these physical quantities at (m − 1) th interface with those at the (m) th interface, that is, the top and the bottom of m th layer in the model. For this, we substitute z = 0 and z = d m in the equations (14)- (17) to obtain the following matrix equations Here, , with X m being non-singular matrix. The superscript 't ' denotes the transpose of that matrix. For the sake of convenience, the subscript 'm' has been dropped from the quantities on the left side of above equations. The superscripts 'm − 1' and 'm', respectively refer to the physical quantity at (m − 1) th and (m) th interfaces. The non-zero entries of matrices X m and Y m are given by where P m = kd m r am , Q m = kd m r bm , and d m denotes the thickness of m th layer. Eliminating the column matrix of arbitrary constants, that is, M m from equations (18) and (19), we obtain where is fourth-order square matrix. According to Meehan [14], the matrix T m is referred as Transfer matrix because it relates and transfers the characteristics between the top and the bottom of the m th layer. The non-zero entries of matrix T m are given as Writing the equation (20) for (m − 1) th layer, we obtain Owing to equations (20), (22) and imposing the boundary conditions of physical quantities at (m − 1) th interface, we obtain the following linear relationship between the quantities at Using the formulation of transfer matrix repeatedly for the remaining layers, the physical quantities at the last interface, that is, (n − 1) th interface are connected to those at the top interface, that is, (0) th interface through the following matrix equation where is fourth-order square matrix. It is also referred as Propagator matrix in the literature (see [14]). It can be seen that the matrix P propagates the information between the topmost layer and the last layer in the considered model. Owing to equation (18) for m = n and the equation (24), the column matrix of arbitrary constants is given by

Dispersion Relation
Equation (25) is applicable to body waves propagating through a multilayered half-space consisting of evenly distributed voids. To study the propagation of Love-type surface wave in multilayered elastic half-space containing voids, the radiation condition must be satisfied, that is, the disturbance must go on diminishing with increasing depth co-ordinate. Mathematically, the displacement component (v n ) and the change in void volume fraction (φ n ) in the half-space must go on vanishing as z → ∞. Owing to equations (10) and (11) for n th layer, we must have B n = D n = 0, provided the expressions r an and r bn are purely imaginary. That is, these quantities must be of the form r an = −ιR a and r bn = −ιR b , where are real and positive, for which we must have c < {a n , b n }. Under these assumptions, the matrix equation (25) is expressed as where N 2l = −ι(kμ n r bn ) −1 P 3l , N 4l = −ι(kα n r an ) −1 P 4l , l = 1, 2, 3, 4.
In the following, we shall derive the dispersion equations for two different types of topmost boundary surface in the model.

Stress-Free Boundary Surface
Here, the topmost boundary surface, that is, (0) th interface is considered to be free from mechanical stresses. Therefore, the appropriate boundary conditions shall be the vanishing of force stress tensor and equilibrated stress tensor. Mathematically, these conditions at the (0) th interface are written as With these considerations, equation (26) reduces to the linear homogeneous system of equations in v 0 and φ 0 , given by Owing to the expressions of elements of matrix P , it can be noticed that the elements P 12 , P 21 , N 22 , and N 41 vanish. Thus, the condition for non-trivial solution of the above homogeneous system of equations yields This is the dispersion relation for Love-type surface wave propagating in the n-layered model with stress-free boundary. Using (27) into this relation, we obtain the following two dispersion equations kμ n r bn P 11 − ιP 31 = 0, kα n r an P 22 − ιP 42 = 0.
Each of these dispersion equations shall provide a wave front of propagating Love-type wave in the model. Thus, we conclude that there exist two fronts for Love-type wave progressing through the considered model. It can also be noticed that (i) both the wave fronts are dispersive in nature, (ii) the wave front given by equation (28) 1 is independent of the void parameters in the model and completely aligns with the dispersion relation for Love-type wave propagating through the multilayered elastic half-space (see Haskell [8]), and (iii) the second front given by equation (28) 2 is new and appeared due to the presence of voids in the model. It is new in the sense that this wave front is obtained in addition to the classical wave front of Love wave propagation. The additional wave front in the material models has also been found earlier in the studies by Dey et al. [4] and Kaur et al. [10].

Rigid Boundary Surface
When the topmost boundary surface is assumed to be rigid, then the boundary conditions would be the vanishing of transverse displacement and the change in void volume fraction, that is, at the (0) th interface, we have Proceeding on similar way as in case of stress-free boundary, we shall obtain the following two dispersion equations kμ n r bn P 13 − ιP 33 = 0, kα n r an P 24 − ιP 44 = 0.
Here, again the first front of Love-type waves represented by equation (29) 1 is found to be independent of void parameters, while the second front given by equation (29)

Particular Cases
In this section, dispersion equations for Love-type waves have been deduced from the present formulation under special cases. These special models are investigated in detail as follows.

2-Layered Model
Here, we desire to have a two-layered model consisting of a layer of uniform thickness d 1 resting over a half-space such that the topmost boundary surface coincides with the plane z = −d 1 . We shall study the propagation of Love-type waves for two cases, namely (i) stress-free boundary and (ii) rigid boundary in the 2-layered model. These cases can be easily reduced from the general formulation developed under Sects. 5.1 and 5.2. For this purpose, we shall put n = 2, and the propagator matrix P reduces to P = T 1 . The entries of transfer matrix T 1 can be obtained from (21).

Stress-Free Boundary
Substituting n = 2 into equations (28), we obtain following two dispersion equations tan Q 1 = ι μ 2 r b2 μ 1 r b1 , The dispersion equation (30) 1 represents first front of Love-type waves, which is independent of void parameters and exactly matches with the well known dispersion relation of classical Love wave (see Ewing et al. [6]). The dispersion equation (30) 2 corresponds to the second front of Love-type waves, which is new and has appeared due to the presence of voids in the model. This equation is in complete agreement with the equation (41) in the corresponding model obtained by Dey et al. [4].

(A) Conditions of propagation:
The conditions of existence for both the fronts of Love-type waves propagating in the considered model are investigated using the expressions of displacement component and the change in void volume fraction as given in (10) and (11). As discussed under Sect. 5 that the Love-type surface waves diminish with increasing depth in the half-space, therefore, the expressions r a2 and r b2 must be purely imaginary. Using equation (9), we can rewrite these expressions as Here, c I and c I I refer to the speeds of first and the second front respectively of Love-type waves progressing in the model. For the solutions of v 1 and φ 1 to be bounded in the layer, the quantities r b1 and r a1 must be real. This yields the conditions of existence for propagation of Love-type waves in the model as , as has been defined earlier in (9).
We also notice that the second front exists provided We may conclude that (i) the first wave front propagates with speed c I , only if c I is greater than the speed of classical shear wave in the layer and less than the speed of classical shear wave in the half-space, (ii) the second front propagates with speed c I I , only if c I I remains confined between a 1 and a 2 . These conditions of existence for Love-type waves have been verified numerically under Sect. 7.

(B) Surface response:
We shall obtain the expressions of field variables and stresses in 2-layered model having stress-free boundary surface. Owing to (10)- (13), the expressions of transverse displacement, change in void volume fraction, transverse shear stress, and the normal component of equilibrated stress, respectively, in the layer and in the half-space are given as follows In the layer: In the half-space: Now, we impose the boundary conditions at the topmost surface and at the interface between the layer and the half-space, given by The relation amongst the arbitrary constants A 1 , B 1 , A 2 , C 1 , D 1 , and C 2 can be derived by inserting the expressions given in (33) and (34) into these boundary conditions. We note that the arbitrary constants B 1 and A 2 can be written in the terms of constant A 1 , whereas the arbitrary constants D 1 and C 2 can be expressed in terms of the constant C 1 as follows Using (35) into (33) and (34), the plots for the variation of transverse displacement (v), the change in void volume fraction (φ), the shear stress component (Y z ) and the equilibrated stress (h z ) have been depicted against the depth co-ordinate (z) under Sect. 7. The continuity of these components at the interface between the layer and the half-space have been verified graphically by specifying the values of constants A 1 and C 1 . Also, the vanishing of stresses at the topmost boundary surface of the model has been observed through the graphs.

(C) Limiting Cases:
We shall study the behavior of Love-type surface waves propagating in the 2-layered model for the limiting cases of angular frequency (p) and thickness of the layer (d 1 ).
(i) Limiting low frequency: In this case, we assume that p → 0 and then investigate the behavior of both the wave fronts whose dispersion equations are given by (30). We observe that as p → 0, the speed of first wave front tends to the speed of shear wave in the half-space.
To investigate the speed of second wave front under the limit p → 0, we see that in view of the restriction provided in (32), this limiting case cannot be studied. This is because the second wave front exists only after a certain non-zero value of frequency p. We conclude that (ii) Limiting high frequency: Here, we assume that the angular frequency p → ∞. And the dispersion equations (30) 1 and (30) 2 of the first and second wave fronts, respectively, yield Hence, we conclude that under limiting high frequency, the speed of first wave front tends to the speed of shear wave in the layer, whereas the speed of second wave front approaches to the speed of wave associated with the void volume fraction in the layer.
It can be noted from equation (9) that a 2 is a function of angular frequency (p). This indicates that for the limiting case of thickness, the speed of first wave front tends to the speed of shear wave in the half-space, whereas the speed of second wave front approaches to a particular value of a 2 for a fixed value of frequency (p). Thus, for different values of frequency, the limiting value of the speed of second wave front will be different. These three limiting cases have also been verified graphically in the subsequent Sect. 7.

(D) Reduced Cases:
We have already developed the dispersion equations of Love-type surface waves propagating in the 2-layered model with stress-free boundary surface. Also, some characteristics of the two wave fronts thus obtained have been investigated. From this general formulation, we shall obtain a few reduced models by either considering or neglecting the presence of voids in the model.

Elastic layer with voids resting over a classical elastic half-space:
Here, we shall assume a model consisting of an elastic layer with voids overlying a classical elastic half-space. This model can be achieved by ignoring the presence of void parameters, namely, α 2 and ξ 2 from the half-space. In this case, the dispersion equation of first wave front given by (30) 1 shall remain unchanged as it is independent of void parameters, however, the dispersion equation (30) 2 corresponding to the second wave front reduces to Thus, it can be concluded that in the absence of void parameters from the half-space, there still exist two wave fronts of Love-type surface waves.

Classical elastic layered half-space:
Neglecting the presence of voids from both the layer and the half-space of the considered model, one shall be left with the model consisting of an elastic layer resting over an elastic half-space. This model can be achieved from the above considered case by ignoring the void parameters, namely, α 1 and ξ 1 from the layer too. This yields in the vanishing of second wave front and hence, the dispersion equations given through (30) lead to the existence of only first wave front of Love-type wave in the model. The dispersion equation of this only wave front, that is, equation (30) 1 exactly matches with the dispersion equation of classical Love wave given in Ewing et al. [6].
The dispersion equation of this model can also be achieved directly from (24). For this purpose, we substitute n = 2 into the matrix equation (24) to obtain Upon neglecting the presence of voids completely, we note that the second and fourth equations of this system become redundant. And the entries of so reduced matrix T 1 are found to be in complete agreement with those of equation (9.5) relevant to the corresponding problem of Haskell [8].

Rigid Boundary
To obtain the dispersion equations in 2-layered model with rigid boundary surface, we substitute n = 2 into equations (29). The dispersion equations for this model reduce to tan Q 1 = ι μ 1 r b1 μ 2 r b2 , The first front of Love-type waves progressing in the 2-layered model with rigid boundary is given by equation (40) 1 , which is independent of the presence of void parameters. Whereas, the second wave front is given by equation (40) 2 and it depends solely upon the void parameters in the model.

3-Layered Model
We shall study the problem in a model consisting of an elastic half-space with voids overlaid by two uniform elastic layers containing voids that are in welded contact with each other and of finite thicknesses. For this case, we substitute n = 3 and the propagator matrix can be written as P = T 2 T 1 . The entries of transfer matrices T 1 and T 2 can be obtained from equation (21). Further, depending upon the nature of topmost boundary surface in the model, we shall discuss the following two cases.
It has been found that there exist two fronts for Love-type waves propagating in the 3layered model with stress-free boundary surface. It is noted that the dispersion equation (41) 1 corresponding to the first wave front is independent of void parameters present in the model. However, equation (41) 2 depicts that the second wave front is influenced by the presence of voids in the model. Also, it is noted that both the wave fronts are dispersive in nature.
It is clear that both the wave fronts are dispersive in nature. It can also be seen that the first wave front is not influenced by the voids present in the model, while the second wave front is affected by the presence of void parameters.

Numerical Computations
The various characteristics of both the fronts of Love-type surface waves are discussed numerically in this Section. The behavior of both the wave fronts under various circumstances for the 2-and 3-layered models with stress-free boundary surface are considered. Dispersion equations corresponding to these respective models have been derived in Sects. 6.1.1 and 6.2.1 through equations (30) and (41). For graphical illustrations of the theoretical results obtained in the preceding sections, MatLab software has been used. A particular set of values has been assigned to various material parameters in the model that are given in Table 1. Values of the material parameters for topmost layer and for the half-space have been borrowed from Puri and Cowin [17] and Kaur et al. [10]. The values of material parameters   Fig. 2 that the speed of first wave front lies within the range given by 2160.25 m/sec < c I < 2500.02 m/sec for the considered frequency range. Hence, the condition of existence for the first wave front in 2-layered model given by (31) 1 has been verified numerically. It is further noted that for limiting low frequency (f → 0), the value of c I approaches b 2 which is the speed of shear wave in the half-space. For the limiting high frequency case (f → ∞), the value of c I tends to b 1 , that is, the speed of shear wave in the layer. These limiting cases are in complete agreement with the analytical results obtained in (36) and (37). Also, for the 3layered model, the analogous condition of existence has been noticed through the figure. The value of wave speed c I lies in the range given by 2160.25 m/sec = b 1 < c I < b 3 = 2500.02 m/sec for the considered frequency and is independent of the shear wave speed (b 2 ) in the intermediate layer. It can be analyzed that the maximum dispersion occurs in the range 10 2 Hz < f < 10 4 Hz. Outside this range of linear frequency, the value of wave speed does not show significant variation for both the models. However, the numerical value of wave speed for the first front in 3-layered model remains less than that in the 2-layered model throughout the considered frequency range, but worth notable in the range 50 Hz < f < 3000 Hz only. Figure 3 shows the comparison of dispersion curves for the second front of Love-type surface waves in the considered 2-layered and 3-layered models. The values of the speed of second wave front c I I for the given range of frequency (1258.23 Hz < f < 10 5 Hz) have been calculated from the dispersion equations (30) 2 and (41) 2 . The second wave front ceases to propagate at a frequency less than 1258.23 Hz as discussed analytically in Sect. 6.1.1 and also depicted graphically in the subsequent figure. At the limiting high frequency (f → ∞), the value of c I I approaches to 5577.73 m/sec, that is, the speed of wave associated with void volume fraction (γ 1 ) in the layer. Hence, the limiting case discussed through equation (37) has been verified numerically. Further, it is concluded that the wave speed of second front in the model with n = 2 is greater than that in the model with n = 3 throughout the considered range of frequency, but it is numerically notable in the range 1258.23 Hz < f < 6500 Hz. For both the models, similar pattern of dispersion curves is observed. As the value of frequency increases, the value of speed of second wave front decreases monotonically to attain a certain minimum value towards the end of the considered frequency range. However, as compared to the speed of first front of Love-type surface waves, the speed of second wave front for both the models is enormously high. This indicates that the second front of Lovetype surface waves, which has appeared due to the presence of voids in the model, travels at a faster speed. This wave front reaches the observing station prior to the first wave front and hence may act as an alarming signal for a possible earthquake, if the model pertains to seismic waves propagating through the layered model of Earth. Figure 4 illustrates the condition of existence for the second front of Love-type surface waves progressing in the 2-layered model for the frequency range, 1000 Hz < f < 12000 Hz. This condition of propagation has already been derived analytically in Sect. 6.1.1. It is evident from the figure that this wave front ceases to propagate for the value of frequency less than 1258.23 Hz. This validates the condition of existence for the second wave front obtained in (31) 2 , as it is computed that max 1 2π Further, it is observed that the speed of second wave front (c I I ), which is depicted in blue dashed curve in the figure remains confined between the red solid curve (a 1 ) and the pink solid curve (a 2 ). Next, we shall study the effect of thickness of layer on the speed of both fronts of Lovetype surface waves in 2-layered model at various different frequencies. Figure 5 depicts the behavior of wave speed of the first front c I in the given range of thickness of the layer (0 m < d 1 < 10 m) at three specific frequencies, that is f = 0.1, 1 and 5 kHz. For each considered frequency, it has been concluded that with increase in the thickness of layer, the value of wave speed c I gradually decreases to a certain value of d 1 (= d * ) and thereafter, the speed remains almost constant. However, the value of d * enhances as the value of frequency decreases. Also, it is evident graphically that as d 1 → 0, the value of speed of first wave front c I → b 2 (= 2500.02 m/sec), which is the speed of shear wave in the half-space. This observation is also derived analytically through equation (38) 1 . For a fixed value of d 1 , it has been found that the value of speed c I decreases with increase in the frequency. Figure 6 shows the variation of wave speed of second front c I I with the thickness of the layer d 1 at three different values of frequency, f = 2, 5 and 10 kHz. For each of these fixed values of f , the value of c I I decreases with increase in the value of d 1 for given range  Further, we have investigated the effect of void parameters, namely, α i and ξ i (i = 1, 2) on the speed of second wave front for the considered 2-layered model through Figs. 7 and 8. The variation of the wave speed c I I with respect to void parameters of the layer and of the half-space has been shown graphically. In Fig. 7, c I I has been plotted against the void parameter α in the given range, 5.5 GPam 2 < α < 8.5 GPam 2 at a particular value of frequency, that is, f = 2 kHz. While plotting against the void parameter of the layer, namely α 1 , the value of void parameter of the half-space is kept constant, that is, α 2 = 8 GPam 2 . It can be seen that the value of c I I increases monotonically with increase in the value of α 1 in the given range. Similarly, the plot of c I I is depicted with respect to α 2 for the fixed value of α 1 (= 7 GPam 2 ) and same pattern is observed. Figure 8 depicts the variation of the speed of second wave front c I I (m/sec) with the void parameter ξ of the layer and of the half-space in the given range, 9 GPa < ξ < 12 GPa. The value of linear frequency f is taken as 2 kHz. While plotting the curve corresponding to void parameter of the layer, that is, ξ 1 , the value of ξ 2 is kept constant (= 12 GPa). It has been found that the value of wave speed increases monotonically as the value of ξ 1 increases. Similar pattern is observed for the curve corresponding to ξ 2 , fixing the value of ξ 1 = 11 GPa. From Figs. 7 and 8, it can be clearly seen that the effect of voids parameters, namely, α i and ξ i on the wave speed of second front is appreciable. Thus, we believe that second front of Love-type surface waves may be of great importance to the seismologists and other researchers in the relevant area of research.    Figure 10 depicts the variation of the change in void volume fraction (|φ|) against the depth parameter (z) at different values of linear frequency, namely, f = 1.5, 2, 5 and 10 kHz. The pattern of all the curves is almost same as in the case of displacement component. Also, the continuity of |φ| at the interface of the layer and the half-space, that is, z = 0 has been verified numerically. Thus, it can be concluded from Figs. 9 and 10 that both the kinematical variables present in this theory vary similarly with respect to (i) depth co-ordinate (z) and (ii) linear frequency f in the considered 2-layered model. The absolute values of v and φ decrease continuously with increase in the depth parameter. This has also been accounted for the theoretical development in order to satisfy the radiation condition, that is, the amplitude of surface wave diminishes with depth parameter.
In Fig. 11, the variation of transverse shear stress component (|Y z |) is depicted with depth co-ordinate (z) for the same set of values of frequency that have been considered in the case of displacement component. While, Fig. 12 illustrates the variation of z-component of equilibrated stress (|h z |) with depth co-ordinate (z) for the same set of frequency taken in the case of |φ|. The considered range of depth for both the plots is −1 m < z < 2 m. From these figures, we conclude that (i) Both the force stress and equilibrated stress vanish at the boundary surface.
(ii) Both of these stresses are continuous at the interface between the layer and the halfspace.
(iii) At a specific value of frequency, all the stress curves follow similar pattern, that is, starting from zero value at the boundary surface (z = −1) they attain a maximum value while reaching the interface (z = 0). However, after attaining peak at the interface, the value of the stress decreases monotonically with increase in the depth parameter in the half-space.
(iv) The stress curves are steeper for higher values of frequency in the layer, while the curves tend to become flat for smaller values of frequency.  (v) In the half-space, a sharp decay is observed across the interface for higher values of frequency. However, for smaller values of f this decay is comparatively slow.
Conclusions (i) and (ii) have been accounted for the corresponding boundary conditions at the topmost surface, that is, free from mechanical stresses, and at the interface between the layer and the half-space.

Conclusions
Haskell matrix method has been employed to study the propagation of Love-type surface waves in multilayered elastic solid half-space containing voids. Two different models are considered depending upon the nature of topmost boundary, namely, stress-free and rigid boundary surface. By imposing appropriate boundary conditions, the dispersion equations have been derived for these models. Numerical computations have been performed for 2layered and 3-layered models with stress-free boundary surface. We may conclude that 1. There exist two fronts of Love-type surface waves propagating with distinct speeds that are dispersive in nature. 2. One of the wave fronts is independent of void parameters, while the other wave front is new and appeared due to the presence of voids in the model. 3. The first wave front in the stress-free boundary surface model is analogous to that obtained earlier by Haskell [8]. 4. It has been noticed that the speed of second wave front is significantly higher than that of the first wave front. Due to the faster speed of second front, it is expected to arrive much earlier than the first wave front at the observing station. The second wave front, if captured, can be considered as an alarming signal for the occurrence of a possible earthquake. Since the second front has appeared purely due to the presence of voids, hence the study of surface waves in material with voids can play an important role in earthquake prediction. 5. For the 2-layered model, the second wave front ceases to propagate below a certain value (non-zero) of frequency. 6. The presence of intermediate layer in the 3-layered model lowers down the speed of both the wave fronts, however the effect is significant on the second front of Love -type wave. 7. With the increase in the thickness of layer resting over the half-space, the speeds of both the wave fronts decrease. 8. The presence of void parameters, namely, α i and ξ i (i = 1, 2) in the model considerably affects the speed of second front of Love-type surface waves.