Stress Wave Propagation through Rock Joints Filled with Viscoelastic Medium Considering Di ﬀ erent Water Contents

: A rock mass often contains joints ﬁlled with a viscoelastic medium of which seismic response is signiﬁcant to geophysical exploration and seismic engineering design. Using the propagator matrix method, an analytical model was established to characterize the seismic response of viscoelastic ﬁlled joints. Stress wave propagation through a single joint highly depended on the water content and thickness of the ﬁlling as well as the frequency and incident angle of the incident wave. The increase in the water content enhanced the viscosity (depicted by quality factor) of the ﬁlled joint, which could promote equivalent joint sti ﬀ ness and energy dissipation with double e ﬀ ects on stress wave propagation. There existed multiple reﬂections when the stress wave propagated through a set of ﬁlled joints. The dimensionless joint spacing was the main controlling factor in the seismic response of the multiple ﬁlled joints. As it increased, the transmission coe ﬃ cient ﬁrst increased, then it decreased instead, and at last it basically kept invariant. The e ﬀ ect of multiple reﬂections was weakened by increasing the water content, which further inﬂuenced the variation of the transmission coe ﬃ cient. The water content of the joint ﬁlling should be paid more attention in practical applications.


Introduction
The study on the seismic response of the rock mass is of great significance to geophysical exploration and seismic engineering design. Usually, a rock mass contains joints that are discontinuous interfaces and that significantly influence mechanical behaviors of the rock mass [1][2][3][4][5][6][7]. Due to the weathering effect, the natural joint was often filled with loose and soft materials such as sand with pore space containing normal air with high moisture content which significantly controlled the joint mechanical behaviors [1,2,[8][9][10]. It was found that the thickness and water content of the filling had great influences on the seismic response of the filled joint [11][12][13]. It was found that saturated soils behaved viscoelastic in response to the dynamic loads, which can be characterized by the Kelvin viscoelastic model [14,15]. Previous studies also showed that the viscosity had an important effect on the dynamic response of the joint and could also dissipate wave energy [16][17][18][19][20][21][22][23]. Meanwhile, when the density of the filled material could not be neglected compared with the rock density, the initial mass of the filled joint could affect wave propagation [20]. Therefore, it is crucial to investigate the seismic response of the joints filled with a viscoelastic medium.
Many analytical models were proposed to characterize the seismic response of the filled joints. Ref. [24] applied the displacement discontinuity model (DDM) to study wave propagation through a filled rock joint. The model assumed that the stress was continuous, but the displacement was not when the stress wave impinged on the filled joint which was controlled by the fracture compliance [25]. Based on the DDM, [26] proposed a three-phase medium model for the filled rock joint considering the interaction of the three-phase system consisting of gas, liquid and granular solid. Ref. [20] proposed a displacement and stress discontinuity model (DSDM) to study wave propagation through filled joints considering the initial mass effect of the filled material. This model assumes that not only the displacement but also the stress is discontinuous in the front and rear interface of the filled joint. However, the above studies mainly have two limitations: (1) Wave reflection in the filled layer was neglected when it needs to be considered as the ratio of filled thickness to the wavelength is relatively large; (2) There are many parameters in these analytical models which are difficult to determine by laboratory tests, such as the viscosity and stiffness of the filled joint. To overcome these problems, a simple and feasible way is to model the filled joint as a layer of independent continuum solid medium. Properties of the filled joint are controlled by those of the filled medium such as the density, filled thickness (aperture), wave velocity and quality factor (Q) that can be readily determined. Furthermore, multiple reflections in the filled layer can be taken into account. It was found that a joint or an interface between two solids can be represented as a thin-layer interface [27]. Ref. [28] pointed out that a thin-layer interface or a zero-thickness interface for a joint should be essentially the same from the physical point of view. The thin-layer medium model (TLMM) was also proved to be reasonable and feasible to characterize stress wave propagation through filled joints in previous studies [29]. They found that predicted results of stress wave propagation through joints by the purely elastic TLMM are nearly the same as the zero-thickness interface model if the ratio of the wavelength to the filled thickness is large enough. Zhu et al. in 2012 [30] investigated the propagation of the stress waves normally impinging on a single rock joint filled with a viscoelastic medium. By modeling the interface between two solids as a thin viscoelastic layer with stiffness and inertia term, wave propagation through a single filled joint was addressed by [31]. However, thus far, analytical models based on the TLMM were rarely reported to investigate seismic response of viscoelastic filled joints as the water content of the filling varied.
The purpose of this paper is to establish an analytical model on the seismic response of the rock joints filled with a viscoelastic medium, considering the water content, the filled thickness (aperture), viscosity of the filled material (quality factor), the joint spacing (related to joint density) as well as the incident angle and frequency of the incident wave.

Method
When the stress wave impinges on a viscoelastic filled joint, only one portion of the stress wave can pass through the filled joint while the other portion will be reflected. Meanwhile, wave conversion happens for the case of the oblique incidence at the joint [18], as shown in Figure 1. There are multiple reflections at joint interfaces [32,33]. The overall transmitted wave, after propagating through a set of filled joints, is the superposed wave consisting of the first-arriving transmitted wave and the later-arriving transmitted waves induced by multiple reflections at each joint. In previous studies, three types of analytical methods were mainly proposed to solve the multiple-reflection problem, i.e., the method of characteristics (MC) [33,34]; the virtual wave source method (VWSM) [35,36]; the recursive method (RM) [21][22][23]30]. MC was mainly applied to solve the 1-D problem of stress wave propagation through joints. For the VWSM, it is necessary to clarify the analytical solution of the wave propagation through a single joint. The propagator matrix method (PMM) is one of the RMs, that was widely applied to solve wave propagation problems of stratified medium [37,38]. The algorithm needs not to identify the analytical solution of wave propagation through a single interface or joint. Ref. [39] first introduced the PMM to study the propagation of the stress wave through unfilled joints satisfying purely elastic deformational behavior. Ref. [21] extended the PMM from the purely elastic case to the viscoelastic case and studied the seismic response of the rock mass containing viscoelastic unfilled joints. Combining with the TLMM, the PMM is also used to derive formulas of wave propagation through rock mass containing joints filled with viscoelastic medium in this paper.
formulas of wave propagation through rock mass containing joints filled with viscoelastic medium in this paper.
As shown in Figure 1, a 2-D plane is considered and both the local coordinates (xi, zi) and the global coordinates (X, Z) are applied. The rock mass is divided by the multiple parallel filled joints with identical joint spacing hr as shown in Figure 1a. The filled joints are numbered from 1 to n and the rock layers from 1 to n + 1, respectively, with Z increasing towards the lower layers. When a P wave impinges on the upper boundary, wave conversions, reflection and transmission occur. θ denotes the incident and the reflection angle of the incident P wave and φ represents the reflection angle of the converted SV wave. An element surrounded by the dash frame (see Figure 1a) is chosen to derive the stress and displacement relations between two adjacent rock layers. The sandwiched element consists of two layers rocks and one filled layer with a thickness of hf, which are numbered as i, f and i + 1 from the upper to the low layer, respectively. Each rock layer is assumed to be a linearly elastic, homogeneous and isotropic medium with the density of ρr, P-and S-wave velocity of Vp r and Vs r . Deformation of each filled medium is characterized by the Kelvin viscoelastic model, which has the density of ρf, P-and S-wave velocity of Vp p and Vs p , and P-and S-wave quality factor of Qp and Qs, where Vp p and Vs p denote the purely elastic wave velocity when the viscosity is not considered. The interface between the rock and the filled h f + - As shown in Figure 1, a 2-D plane is considered and both the local coordinates (x i , z i ) and the global coordinates (X, Z) are applied. The rock mass is divided by the multiple parallel filled joints with identical joint spacing h r as shown in Figure 1a. The filled joints are numbered from 1 to n and the rock layers from 1 to n + 1, respectively, with Z increasing towards the lower layers. When a P wave impinges on the upper boundary, wave conversions, reflection and transmission occur. θ denotes the incident and the reflection angle of the incident P wave and ϕ represents the reflection angle of the converted SV wave. An element surrounded by the dash frame (see Figure 1a) is chosen to derive the stress and displacement relations between two adjacent rock layers. The sandwiched element consists of two layers rocks and one filled layer with a thickness of h f , which are numbered as i, f and i + 1 from the upper to the low layer, respectively.
Each rock layer is assumed to be a linearly elastic, homogeneous and isotropic medium with the density of ρ r , P-and S-wave velocity of V p r and V s r . Deformation of each filled medium is characterized by the Kelvin viscoelastic model, which has the density of ρ f , P-and S-wave velocity of V p p and V s p , and P-and S-wave quality factor of Q p and Q s , where V p p and V s p denote the purely elastic wave velocity when the viscosity is not considered. The interface between the rock and the filled medium is assumed to be welded, in which both stresses and displacements are continuous. Two local coordinates x f -o f -z f and x i+1 -o i+1 -z i+1 are applied as shown in Figure 1b. The symbols "+" and "−" represent the upper boundary and the lower boundary of a layer (e.g., in the (i + 1) th layer, "+" represents the upper boundary (z i+1 = 0), whereas "−" represents the low boundary (z i+1 = h r ). For the linearly elastic rock layer, the 2-D wave equation can be written as: where u i denotes the particle displacement and σ ij denotes the stress, λ r and µ r are Lamé constants for rock and λ r = ρ r (V p r ) 2 -2ρ r (V s r ) 2 , µ r = ρ r (V s r ) 2 .
For the Kelvin viscoelastic filled layer, the 2-D wave equation can be written as: where λ p and µ p are the purely elastic Lamé constants of the filled layer when the viscosity is not considered, λ p = ρ f (V p p ) 2 -2ρ f (V s p ) 2 , µ p = ρ f (V s p ) 2 and λ f and µ f denote the viscoelastic modulus for the filled medium. Equations (1)-(8) are transformed into frequency domains by Fourier transform. For the rock layer, where, variables with the symbol "~" over them denotes the Fourier transform. For the filled layer, Appl. Sci. 2020, 10, 4797

of 19
where ω is angular frequency and i = √ −1. The quality factors related to energy attenuation can be defined as [40]: When Q p →∞ and Q s →∞, Equations (13)-(16) represent the pure elastic case without viscosity energy loss. When Q p →0 and Q s →0, the filled medium has a very large viscosity. Based on Equations (15)-(17), we can obtain corresponding complex equivalent Lamé constants and an equivalent wave velocity of the filled medium as follows: Substituting the Equation (18) into Equations (15) and (16), we can obtain the following equations: It can be seen that Equations (20) and (21) have the same form as those of purely elastic wave equations. Based on the Helmholtz decomposition theory [41], the 2-D dynamic equation can be expressed as two displacement potential functions. Assume that the ϕ i (the displacement potential of the P wave) and ψ 2i (the second component of the displacement potential of the S wave) are the potential functions of the ith layer in the frequency domain which can be expressed by following two Equations [21,39].
ϕ i (x, z, ω) = U p i e i(−k xi x+k zpi z) + D pi e i(−k xi x−k zpi z) (22) ψ 2i (x, z, ω) = U si e i(−k xi x+k zsi z) + D si e i(−k xi x−k zsi z) where k xi , k zpi are the P-wave numbers of the ith layer in the x and z directions, and k zsi is the S-wave number of the ith layer in the z direction. As shown in Figure 1b, the sandwiched element consists of two types of materials in turn. When the wave propagates in the rock layer, the wave numbers k xi , k zpi and k zsi can be expressed as k x r = k p r sin θ, k zp r = k p r cos θ, k r zs = (ω/V r s ) 2 − k r 2 x , where k p r = ω/V p r is the wave number of the P wave. Similarly, when the wave propagates in the filled layer, the wave numbers k xi , k zpi and k zsi can be expressed as k is the wave number of the P-wave propagation through the filled medium.
The particle displacement can be written as: Appl. Sci. 2020, 10, 4797 6 of 19 Substituting Equations (22) and (23) into Equations (24) and (25), we can form: Based on the constitutive Equations (20) and (21) and Equations (26) and (27), we can also form: Equations (26)- (29) can be written as the matrix form: For the rock layer, For the filled layer, and Equation (30) establishes the relation between stress and displacement and the potential amplitudes in each layer. Note that the matrix Q i is related to z coordinate. As shown in Figure 1b, two local coordinate systems (x f -o f -z f and x i + 1 -o i+1 -z i+1 ) are adopted, respectively. For the upper boundary of a layer, z = 0 and Q i can be expressed as Q i + . For the low boundary of a layer, z = h i and Q i can be expressed as Q i − . Therefore, we can form: The relation of the potential amplitudes between two adjacent rock layers can be derived from the sandwiched element as shown in Figure 1b. As mentioned before, the stress and the displacement at the interface between the rock layer and the filled layer are continuous. In addition, based on Equation (30), we can easily form: Appl. Sci. 2020, 10, 4797 8 of 19 From (31)-(33), we can form: Note that R i can be regarded as the stiffness matrix of the filled joint which is related to the thickness, density, P-and S-wave velocity, and viscosity (quality factor) of the filled medium. It has the similar form compared with that of the unfilled joints [21,39], which additionally considers the effect of the wave impedance of the filling.
The relation of Equation (34) can be extended from the first layer to the nth layer by chain rule: When the wave impinges on the lower boundary of the first layer, the propagation distance is 0, therefore: Meanwhile, after the wave propagating through n filled joints, the wave propagates into the (n + 1)th layers in fact. The displacement and stress vector at the upper boundary of the (n + 1) rock layer should be expressed.
From Equations (37) and (38), we can form: Substituting Equations (36) and (39) into (35), we can form: Appl. Sci. 2020, 10, 4797 9 of 19 For the (n + 1)th layer, there only exists down-going P-and S-waves, therefore U p = 0 and U s = 0. Because only the P-wave impinges on the first layer, D s is equal to 0. Based on Equation (40), the reflected and transmitted waves can be calculated. We can obtain reflection-transmission coefficients that are defined as ratios of reflected and transmitted wave amplitudes to the incident wave amplitude. Specifically, the reflection-transmission coefficients can be written as R p→p = max(abs(U p 1 )), R p→s = max(abs(U s 1 )), T p→p = max(abs(D p n+1 )) and T p→s = max(abs(D s n+1 )). For these symbols, the first subscript letter represents the incident wave and the second subscript letter stands for the transmitted or reflected wave. Details can be found in [21]. During calculation, the wave has to be first transferred into frequency domain signals by the fast Fourier transform (FFT) to calculate the reflected and transmitted waves. Finally, these waves in the frequency domain are transferred into signals in the time domain by the inverse fast Fourier transform (IFFT). In the following section, the experimental incident wave is first used as the input incident wave for model validation and then one cycle of sinusoidal wave is used as the input incident wave for parametric studies. According to [21], the start time to oscillate needed to have enough delay rather than from 0 to eliminate the aliasing effect as much as possible in the processing of FFT. Meanwhile, some amplitude compensations are also needed.

Experimental Calibration
In this section, the analytical model is calibrated according to the experimental data from [12]. In their study, a layer of quartz sand with a fixed thickness of 2 mm but varied water content of 0%, 2%, 5%, 10%, and 15% (saturated) was sandwiched between two norite bars to simulate filled joints with different water contents. The stress wave was induced by a split Hopkinson rock bar (SHRB) system. The quartz sand had a porosity of around 0.4 and material density of 2620 kg·m −3 (solid density). The density of fillings with different water contents can be calculated and listed in Table 1. The norite had a P-wave velocity of 6000 m·s −1 and density of 2900 kg·m −3 [12,42]. The S-wave velocity of the norite was taken as 2/3 of its P-wave velocity (4000 m·s −1 ). The experimental incident wave was adopted as the incident wave to predict transmitted waves after propagating through filled joints with different water contents. The start time of the wave signal to oscillate has a delay of around 0.0012 s compared with that in [12]. Table 1 shows the calibrated parameters in cases of different water contents. Except the density of the filled medium, other parameters all decrease with increasing water content. According to Section 2, the decrease in the quality factor indicates the increase in the viscosity of the filled joint.  Figure 2 shows the experimental incident wave, transmitted waves through a single filled joint with different water contents and predicted transmitted waves. The amplitude of the experimental transmitted wave obviously decreases with increasing water content. Predicted transmitted waves are basically consistent with experimental results except that their durations are a bit larger than those in the experiment. Figure 3 shows a comparison of the transmission coefficient between the experimental and predicted results. We can see that the experimental transmission coefficient decreases with increasing water content, which agrees well with the experimental result. From

Parametric Studies
In this section, reflection-transmission coefficients when a P wave obliquely impinges on a single filled joint and set parallel filled joints with identical spacing were investigated, respectively, considering the incident angle θ and the frequency f of the incident wave, as well as the filled thickness hf and the quality factor of the filled material Qp and Qs. One of these parameters was studied while other parameters were fixed. Parameters of the rock material are taken as those of the norite introduced in Section 3.1. Parameters of the filled medium for different water contents are taken according to Table 1.

Parametric Studies
In this section, reflection-transmission coefficients when a P wave obliquely impinges on a single filled joint and set parallel filled joints with identical spacing were investigated, respectively, considering the incident angle θ and the frequency f of the incident wave, as well as the filled thickness hf and the quality factor of the filled material Qp and Qs. One of these parameters was studied while other parameters were fixed. Parameters of the rock material are taken as those of the norite introduced in Section 3.1. Parameters of the filled medium for different water contents are taken according to Table 1.

Parametric Studies
In this section, reflection-transmission coefficients when a P wave obliquely impinges on a single filled joint and set parallel filled joints with identical spacing were investigated, respectively, considering the incident angle θ and the frequency f of the incident wave, as well as the filled thickness h f and the quality factor of the filled material Q p and Q s . One of these parameters was studied while other parameters were fixed. Parameters of the rock material are taken as those of the norite introduced in Section 3.1. Parameters of the filled medium for different water contents are taken according to Table 1.

Case of One Single Filled Joint
When a P wave obliquely impinges on a single filled joint, there are reflected P and SV waves as well as transmitted P and SV waves as shown in Figure 4.

Case of One Single Filled Joint
When a P wave obliquely impinges on a single filled joint, there are reflected P and SV waves as well as transmitted P and SV waves as shown in Figure 4.  Figure 5 shows reflection-transmission coefficients versus the filled thickness hf considering different water contents of the filled medium. It displays that Rp→p first increases abruptly with the increase in hf, then its increase tends to be gentle while Tp→p decreases basically in a negative exponential form. For the conversion waves, Rp→sv first increases sharply with hf increasing then it basically increases in a linear manner while Tp→sv basically keeps invariant. In general, Rp→sv varies more violently with the increase in hf compared with Tp→sv. As the water content increases, Tp→p decreases evenly, and corresponding curves become divergent; Rp→p and Rp→sv have an increasing trend on the whole; Tp→sv curves for different water contents basically overlap with each other.    Figure 5 shows reflection-transmission coefficients versus the filled thickness h f considering different water contents of the filled medium. It displays that R p→p first increases abruptly with the increase in h f , then its increase tends to be gentle while T p→p decreases basically in a negative exponential form. For the conversion waves, R p→sv first increases sharply with h f increasing then it basically increases in a linear manner while T p→sv basically keeps invariant. In general, R p→sv varies more violently with the increase in h f compared with T p→sv . As the water content increases, T p→p decreases evenly, and corresponding curves become divergent; R p→p and R p→sv have an increasing trend on the whole; T p→sv curves for different water contents basically overlap with each other.

Case of One Single Filled Joint
When a P wave obliquely impinges on a single filled joint, there are reflected P and SV waves as well as transmitted P and SV waves as shown in Figure 4.  Figure 5 shows reflection-transmission coefficients versus the filled thickness hf considering different water contents of the filled medium. It displays that Rp→p first increases abruptly with the increase in hf, then its increase tends to be gentle while Tp→p decreases basically in a negative exponential form. For the conversion waves, Rp→sv first increases sharply with hf increasing then it basically increases in a linear manner while Tp→sv basically keeps invariant. In general, Rp→sv varies more violently with the increase in hf compared with Tp→sv. As the water content increases, Tp→p decreases evenly, and corresponding curves become divergent; Rp→p and Rp→sv have an increasing trend on the whole; Tp→sv curves for different water contents basically overlap with each other.   Figure 6 shows the reflection -transmission coefficients varying with wave frequency f. We can see that T p→p decreases basically in a negative exponential form with increasing f. R p→p and R p→sv first increase abruptly with increasing f, then they only have slight changes. Compared with the other three types of coefficients, T p→sv is relatively insensitive to f. Cases of different water contents present similar variation of reflection-transmission coefficients versus f except for the magnitude difference. As the water content increases, T p→p decreases evenly and corresponding curves become divergent; Both R p→p and R p→sv generally have an increasing trend of which the variation is smaller than that of T p→p ; T p→sv curves for different water contents basically overlap with each other which indicates that it is independent of the water content of the filled medium.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 21 three types of coefficients, Tp→sv is relatively insensitive to f. Cases of different water contents present similar variation of reflection-transmission coefficients versus f except for the magnitude difference.
As the water content increases, Tp→p decreases evenly and corresponding curves become divergent; Both Rp→p and Rp→sv generally have an increasing trend of which the variation is smaller than that of Tp→p; Tp→sv curves for different water contents basically overlap with each other which indicates that it is independent of the water content of the filled medium.  Figure 7 shows the reflection-transmission coefficients versus the incident angle θ when a P wave obliquely impinges on a single filled joint considering the different water contents. It shows that Tp→p does not change much with increasing θ until θ approaches the critical angle 90°, and it suddenly decreases to 0. When θ is smaller than around 40°, Rp→p first decreases with increasing θ. As θ is larger than around 40°, it gradually increases to a peak and then it deceases instead. When θ is larger than around 83°, Rp→p increases rapidly. Tp→sv changes gently with the increase in θ. Rp→sv first increases to the maximum value at θ of approximately 40° and then it decreases with increasing θ instead. Generally, Tp→sv is much smaller than Tp→p. As the water content increases, Tp→p decreases evenly and related curves become divergent when θ is smaller than 83°; When the water content increases from 0% to 5%, both Rp→p and Rp→sv have an obvious increase while they have very limited changes as the water content increases from 5% to 10%; Tp→sv is basically independent on the water content of the filled medium.  Figure 7 shows the reflection-transmission coefficients versus the incident angle θ when a P wave obliquely impinges on a single filled joint considering the different water contents. It shows that T p→p does not change much with increasing θ until θ approaches the critical angle 90 • , and it suddenly decreases to 0. When θ is smaller than around 40 • , R p→p first decreases with increasing θ. As θ is larger than around 40 • , it gradually increases to a peak and then it deceases instead. When θ is larger than around 83 • , R p→p increases rapidly. T p→sv changes gently with the increase in θ. R p→sv first increases to the maximum value at θ of approximately 40 • and then it decreases with increasing θ instead. Generally, T p→sv is much smaller than T p→p . As the water content increases, T p→p decreases evenly and related curves become divergent when θ is smaller than 83 • ; When the water content increases from 0% to 5%, both R p→p and R p→sv have an obvious increase while they have very limited changes as the water content increases from 5% to 10%; T p→sv is basically independent on the water content of the filled medium. Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 21 Figure 7. Reflection-transmission coefficients versus the incident angle θ when the P wave impinges on a single filled joint with a thickness of 10 mm. The incident wave has the frequency of 50 Hz. Figure 8 shows the reflection-transmission coefficients versus the quality factor when a P wave impinges on a single filled joint with a thickness of 40 mm. For convenience, a variable Q is introduced to express Qp and Qs as Qp = Q and Qs = 2/3Q, respectively. To only investigate the influence of the quality factor of the filled medium, the density, thickness and wave velocity of the filled medium were fixed as constant values (Vp p = 185 m•s -1 , Vs p = 123 m•s -1 and ρf = 1807 kg•m -3 ). It shows that Tp→p decreases abruptly with Q increasing, after achieving a minimum value, it increases with an increment of Q instead and eventually it basically keeps invariant. Tp→sv does not change much with increasing Q. Both Rp→p and Rp→sv increase rapidly with Q increasing and then they increase gently, and they basically keep invariant as Q is larger than 20. Generally, Q has a more significant effect on Tp→p, Rp→p and Rp→sv than Tp→sv.  Figure 8 shows the reflection -transmission coefficients versus the quality factor when a P wave impinges on a single filled joint with a thickness of 40 mm. For convenience, a variable Q is introduced to express Q p and Q s as Q p = Q and Q s = 2/3Q, respectively. To only investigate the influence of the quality factor of the filled medium, the density, thickness and wave velocity of the filled medium were fixed as constant values (V p p = 185 m·s −1 , V s p = 123 m·s −1 and ρ f = 1807 kg·m −3 ). It shows that T p→p decreases abruptly with Q increasing, after achieving a minimum value, it increases with an increment of Q instead and eventually it basically keeps invariant. T p→sv does not change much with increasing Q. Both R p→p and R p→sv increase rapidly with Q increasing and then they increase gently, and they basically keep invariant as Q is larger than 20. Generally, Q has a more significant effect on T p→p , R p→p and R p→sv than T p→sv .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 21 Figure 8. Reflection-transmission coefficients of a P wave impinging on a single filled joint versus the quality factor Q. The incident P wave has the frequency of 50 Hz and an incident angle of 20°.
Due to the viscosity of the filled medium, the wave energy dissipated after propagating through the filled joint, which can be defined as follows [21,40]: Due to the viscosity of the filled medium, the wave energy dissipated after propagating through the filled joint, which can be defined as follows [21,40]: Figure 9 shows the energy loss versus the quality factor Q corresponding to Figure 8. It displays that the energy loss E loss first increases abruptly with the increase in Q. After achieving a peak, E loss decreases rapidly with an increment of Q. When Q is larger than 25, E loss changes gently with Q increasing. Figure 8. Reflection-transmission coefficients of a P wave impinging on a single filled joint versus the quality factor Q. The incident P wave has the frequency of 50 Hz and an incident angle of 20°.
Due to the viscosity of the filled medium, the wave energy dissipated after propagating through the filled joint, which can be defined as follows [21,40]: Figure 9 shows the energy loss versus the quality factor Q corresponding to Figure 8. It displays that the energy loss Eloss first increases abruptly with the increase in Q. After achieving a peak, Eloss decreases rapidly with an increment of Q. When Q is larger than 25, Eloss changes gently with Q increasing.

Case of Multiple Parallel Filled Joints
When the stress wave propagates though a rock mass containing a set of joints, multiple reflections occur at each joint [33,34,43]. In this section, involved filled joints are assumed to have the same properties such as filled thickness, density, wave velocity, quality factor of the filled medium, etc. The rock mass is assumed to contain a set of joints with the identical joint spacing. A cycle of sinusoidal P wave with unit amplitude and frequency of 50 Hz was used to be the incident wave. The start time of the wave signal to oscillate has a delay of around 0.1 s. According to [21], the actual central frequency of the incident wave is 42.5 Hz for one cycle of sinusoidal wave with an input frequency of 50 Hz. It is reasonable to use the central frequency of 42.5 Hz to represent the frequency of the incident wave. The ratio between the joint spacing and wavelength have an important effect on the transmission coefficient because of the multiple reflections at filled joints [21,[33][34][35]. Therefore, it is convenient to study the dimensionless joint spacing which is defined as ξ = h r /χ = (h r f )/V p r , where χ is the P-wave wavelength. Figure 10 shows transmitted waveforms (P-wave component) after a P wave obliquely propagates through four parallel filled joints with an identical filled thickness h f of 40 mm and joint spacing h r of 50 m, considering different water contents. It displays that each transmitted wave is made up of several cycles of waveforms arriving at different times induced by multiple reflections at joints. We can clearly see the wave superposition presenting multiple peaks and valleys in turn. As the water content increases, the transmitted wave decreases in the amplitude and the component of the later-arriving wave gradually fades away. When the water content is 15%, it is difficult to observe the later-arriving transmitted wave (after 0.18 s).
spacing hr of 50 m, considering different water contents. It displays that each transmitted wave is made up of several cycles of waveforms arriving at different times induced by multiple reflections at joints. We can clearly see the wave superposition presenting multiple peaks and valleys in turn. As the water content increases, the transmitted wave decreases in the amplitude and the component of the later-arriving wave gradually fades away. When the water content is 15%, it is difficult to observe the later-arriving transmitted wave (after 0.18 s).  Figure 11 shows that variation of transmission coefficients Tp→p and Tp→sv versus ξ after the P wave obliquely propagates through four parallel joints considering different water contents. It shows that there exists two threshold lines ξ1 and ξ2 dividing each curve into three parts, which is similar to results in published papers [21,[33][34][35]44]. Tp→p first decreases with increasing ξ (ξ < ξ1) and then it decreases instead (ξ1 < ξ < ξ2) and at last it basically keeps invariant (ξ > ξ2). With the water content increasing, Tp→p decreases and the position for the first threshold point moves to the left while the second one moves to the right. The increasing part of the curve (ξ < ξ1) shortens while the decreasing part (ξ1 < ξ< ξ2) widens. Tp→sv first increases with increasing ξ, then it decreases after achieving a peak in general. When the water content is low, there exist fluctuations on the curve. As the water content becomes higher, these fluctuations gradually fade away and the curve becomes more and more flat.  Figure 11 shows that variation of transmission coefficients T p→p and T p→sv versus ξ after the P wave obliquely propagates through four parallel joints considering different water contents. It shows that there exists two threshold lines ξ 1 and ξ 2 dividing each curve into three parts, which is similar to results in published papers [21,[33][34][35]44]. T p→p first decreases with increasing ξ (ξ < ξ 1 ) and then it decreases instead (ξ 1 < ξ < ξ 2 ) and at last it basically keeps invariant (ξ > ξ 2 ). With the water content increasing, T p→p decreases and the position for the first threshold point moves to the left while the second one moves to the right. The increasing part of the curve (ξ < ξ 1 ) shortens while the decreasing part (ξ 1 < ξ< ξ 2 ) widens. T p→sv first increases with increasing ξ, then it decreases after achieving a peak in general. When the water content is low, there exist fluctuations on the curve. As the water content becomes higher, these fluctuations gradually fade away and the curve becomes more and more flat. From Figures 10 and 11, it indicates that the effects of multiple reflections at joints on the transmitted wave and transmission coefficient were significantly reduced when the water content increased.
From Figures 10 and 11, it indicates that the effects of multiple reflections at joints on the transmitted wave and transmission coefficient were significantly reduced when the water content increased.

Discussion
It was found that not only the displacement but also the stress at the front and rear interfaces of the filled joint were discontinuous when loaded by the stress wave [20,45]. This attributes to the inertial effect of the mass of the filled medium. Both the stress and the displacement at the interface between the rock and filled material are assumed to be continuous in the current study. It can be shown that the sandwiched element surrounded by two interfaces as shown in Figure 1b can characterize the stress and displacement discontinuities at the filled joint according to Equation (34

Discussion
It was found that not only the displacement but also the stress at the front and rear interfaces of the filled joint were discontinuous when loaded by the stress wave [20,45]. This attributes to the inertial effect of the mass of the filled medium. Both the stress and the displacement at the interface between the rock and filled material are assumed to be continuous in the current study. It can be shown that the sandwiched element surrounded by two interfaces as shown in Figure 1b can characterize the stress and displacement discontinuities at the filled joint according to Equation (34). Therefore, our model can simulate the discontinuous response of stress and the displacement at the filled joint, which has similar physics to DSDM by [20]. Additionally, our model can automatically consider the wave reflection inside the filled layer while DSDM loses the capability. However, the influence of the wave reflection is prominent and cannot be neglected when the ratio of the filled thickness to the wavelength is large enough. As described in Section 3.1, our analytical model was suitable to characterize the seismic response of the filled joint considering different water contents. Properties of the filled joints can be simply described by several common parameters such as the filled thickness, the density, the wave velocity, and the quality factors of the filled medium. These parameters have definite physical meanings and can be obtained relatively more readily compared with those of DSDM [20], i.e., the viscous stiffness.
The water content has a significant effect on the dynamic mechanical characteristics of the filled joint. As shown in Table 1, the increase in the water content can increase the density but reduce the wave velocity and quality factor of the filled medium. As a result, the transmission of the stress wave was significantly influenced. Especially, the variation of the quality factor stands for the variation of the viscosity. As derived in Equations (7) and (8), the viscoelastic wave equations have complex forms composed of a real part and imaginary part. The increase in the imaginary part has double effects, not only can it increase the equivalent dynamic modulus to let more wave energy pass through the filled joint (positive effect), but it can also cause energy loss (negative effect) [21]. For T p→p , the negative effect first plays the domain role and the positive effect plays the domain role in turn when the quality factor increases. For R p→p , R p→sv and T p→sv , the positive effect always plays the domain role. When the quality factor is large, the viscous effect of the filled medium is not prominent. In this case, the reflection-transmission coefficients are independent of the quality factor. We can now understand the variation of the reflection-transmission coefficients versus the quality factor as shown in Figure 8 as well as the energy loss with quality factor in Figure 9. The transmitted waves after propagating through multiple parallel joints are the superposed waves consisting of the first-arriving wave and the later-arriving waves caused by multiple reflections. When the water content increases, the wave reflection increases while the wave transmission decreases (see Figures 5-7). According to this, the multiple reflections at the joint should increase with increasing water content. However, the multiple reflection decreases with increasing water content instead (see Figures 10 and 11). It can be inferred that this phenomenon may result from the increase in the viscosity of the filled medium with the water content increasing, which significantly dissipated the energy of waves by multiple reflections at the filled joints. Thus, the change of the quality factor (viscosity) of the filled medium with increasing the water content has a significant influence on the seismic response of the filled joints which should be paid more attention in practical applications.

Conclusions
We established an analytical model to characterize stress wave propagation through rock joints filled with a viscoelastic medium considering different water contents. Some important conclusions were reached. The seismic response of the single filled joint highly depended on the water content, thickness of the filling as well as the frequency and incident angle of the incident wave, etc. The increase in the water content enhanced the viscosity (depicted by quality factor) of the filled joint, which could promote equivalent joint stiffness and energy dissipation with double effects on stress wave propagation. When the stress wave propagated through a set of filled joints, there existed multiple reflections at the joints. The dimensionless joint spacing was the main controlling factor in the seismic response of the multiple parallel filled joints, which exhibited three-stage variation. The transmission coefficient first increased with the dimensionless joint spacing, then it decreased instead, and at last it basically kept invariant. The increase in the water content could weaken the effect of multiple reflections. It mainly resulted from the increase in the viscosity of the filled joint which significantly dissipated the energy of waves by multiple reflections at the filled joints. The water content of the joint filling should be paid more attention in practical applications.