## 2 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. [PERSON]. [39] first introduced the PMM to study the propagation of the stress wave through unfilled joints satisfying purely elastic deformational behavior. [PERSON]. [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. 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_{\rm 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. \(\theta\) denotes the incident and the reflection angle of the incident P wave and \(\varphi\) 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 \(\rho_{\rm r}\), P- and S-wave velocity of \(V_{\rm p}\)\({}^{\rm r}\) and \(V_{\rm s}\)\({}^{\rm r}\). Deformation of each filled medium is characterized by the Kelvin viscoelastic model, which has the density of \(\rho_{\rm f}\), P- and S-wave velocity of \(V_{\rm p}\)\({}^{\rm r}\) and \(V_{\rm s}\)\({}^{\rm r}\), and P- and S-wave quality factor of \(Q_{\rm p}\) and \(Q_{\rm s}\), where \(V_{\rm p}\)\({}^{\rm r}\) and \(V_{\rm s}\)\({}^{\rm r}\) denote the purely elastic wave velocity when the viscosity is not considered. The interface between the rock and the filled  
Figure 1: (**a**) A P wave impinges on a rock mass containing joints filled with a viscoelastic medium, (**b**) the schematic view of the sandwiched element. medium is assumed to be welded, in which both stresses and displacements are continuous. Two local coordinates \(x_{\rm{r}}\)-\(o_{\rm{r}}\)-\(z_{\rm{f}}\) and \(x_{\rm{i+1}}\)-\(o_{\rm{i+1}}\)-\(z_{\rm{i+1}}\) are applied as shown in Figure 1b.