## 4 Solutions for the Equations  
### Solution for the Soil  
By introducing the concept of local coordinate system and the Laplace transform technique, Equation (1) can be solved to derive the following equation. The derivation procedure can be seen in Supplementary Materials.  
\[U_{ij}=\left[A_{ij}\mathrm{I}_{0}\big{(}q_{ij}r\big{)}+B_{ij}\mathrm{K}_{0} \big{(}q_{ij}r\big{)}\right]\cdot\left[C_{ij}\sin\big{(}h_{ij}z^{\prime}\big{)} +D_{ij}\cos\big{(}h_{ij}z^{\prime}\big{)}\right] \tag{3}\]  
where \(U_{ij}\) is the Laplace transform of \(u_{ij}\); \(\mathrm{I}_{0}\big{(}q_{ij}r\big{)}\) and \(\mathrm{K}_{0}\big{(}q_{ij}r\big{)}\) are modified Bessel functions; \(z^{\prime}=z-h_{i}\); \(A_{ij}\), \(B_{ij}\), \(C_{ij}\), and \(D_{ij}\) are undetermined constants; \(q_{ij}\) and \(h_{ij}\) satisfy the relation  
\[q_{ij}^{2}=\frac{\left\{\eta_{ij}^{2}+\mathrm{i}\Big{[}D_{\mathrm{vi}j}\Big{(} \eta_{ij}^{2}-2\Big{)}+2D_{\mathrm{s}ij}\Big{]}\Big{\}}\right\}h_{ij}^{2}-\left( \omega/v_{ij}\right)^{2}}{1+\mathrm{i}D_{\mathrm{s}ij}} \tag{4}\]  
where \(\eta_{ij}=c_{ij}/v_{ij}\); \(c_{ij}\) and \(v_{ij}\) are longitudinal wave velocity and shear wave velocity of the soil, respectively; \(D_{\mathrm{s}ij}=G^{\prime}_{ij}/G_{ij}\) and \(D_{\mathrm{vi}j}=\lambda^{\prime}_{ij}/\lambda_{ij}\) denote the hysteretic-type damping coefficients.  
The boundary conditions at the upper and lower surfaces of the soil layer can be respectively expressed as  
\[\left.\frac{\partial u_{ij}(r,z^{\prime},t)}{\partial z^{\prime}}\right|_{z^{ \prime}=0}=\left.\left[\frac{k_{(i+1)}u_{ij}(r,z^{\prime},t)}{E_{ij}}+\frac{ \delta_{(i+1)j}}{E_{ij}}\frac{\partial u_{ij}(r,z^{\prime},t)}{\partial t} \right]\right|_{z^{\prime}=0} \tag{5}\]\[\left.\frac{\partial u_{ij}(r,z^{\prime},t)}{\partial z}\right|_{z^{\prime}=l_{i}}= \left.-\left[\frac{k_{ij}u_{ij}(r,z^{\prime},t)}{E_{ij}}+\frac{\delta_{ij}}{E_{ ij}}\frac{\partial u_{ij}(r,z^{\prime},t)}{\partial t}\right]\right|_{z^{\prime}=l_{i}} \tag{6}\]  
Specifically, the stiffness and damping coefficient at the lower surface of the first soil layer tend to infinity according to Assumption (2), while those at the upper surface of the \(n\)th soil layer can be considered to be zero, as it is a free surface.  
Applying the Laplace transform technique to Equations (5) and (6) and together with Equation (3), it can be derived that  
\[\tan\bigl{(}h_{ij}l_{i}\bigr{)}=\frac{\left(\overline{K_{ij}}+\overline{K^{ \prime}_{ij}}\right)h_{ij}l_{i}}{\left(h_{ij}l_{i}\right)^{2}-\overline{K_{ij} }\overline{K^{\prime}_{ij}}} \tag{7}\]  
where \(\overline{K_{ij}}=\frac{k_{ij}+i\omega k_{ij}}{E_{ij}}l_{i}\) and \(\overline{K^{\prime}_{ij}}=\frac{k_{(i+1)j}+i\omega k_{(i+1)j}l_{i}}{E_{ij}} l_{i}\) denote the dimensionless complex stiffness.  
Equation (7) can be solved by means of the bisection method to derive a series of eigenvalues (denoted by \(h_{ijk}\), where \(k=1,\,2,\,\ldots\)). On this basis, \(q_{ijk}\) (\(k=1,\,2,\,\ldots\)) can be obtained according to Equation (4). Further, Equation (3) can be rewritten as  
\[U_{ij}=\sum_{k=1}^{\infty}\left[A_{ijk}\mathrm{I}_{0}\left(q_{ijk}r\right)+B_ {ijk}\mathrm{K}_{0}\left(q_{ijk}r\right)\right]\cos\left(h_{ijk}z^{\prime}- \varphi_{ijk}\right) \tag{8}\]  
where \(\varphi_{ijk}=\arctan\left[\overline{K^{\prime}_{ij}}\right/\left(h_{ijk}l_{i} \right)\right]\). Specifically, the displacement of the outermost zone (i.e., the \(m\)th zone) decays to zero at infinity. So, \(A_{ijk}=0\), when \(j=m\).  
The shear stress of the \(j\)th zone acting on the (\(j-1\))th zone can be expressed as  
\[\tau_{ij}=G_{ij}^{*}\sum_{k=1}^{\infty}\left[A_{ijk}q_{ijk}\mathrm{I}_{1} \left(q_{ijk}r_{ij}\right)-B_{ijk}q_{ijk}\mathrm{K}_{1}\left(q_{ijk}r_{ij} \right)\right]\cos\left(h_{ijk}z^{\prime}-\varphi_{ijk}\right) \tag{9}\]  
where \(G_{ij}^{*}=G_{ij}\bigl{(}1+\mathrm{i}D_{ijl}\bigr{)}\); \(\mathrm{I}_{1}\left(q_{ijk}r_{ij}\right)\) and \(\mathrm{K}_{1}\left(q_{ijk}r_{ij}\right)\) represent modified Bessel functions.