Shear Wave Splitting and Polarization in Anisotropic Fluid-Infiltrating Porous Media: A Numerical Study

Abstract

The triggering and spreading of volumetric waves in soils, namely pressure (P) and shear (S) waves, developing from a point source of a dynamic load, are analyzed. Wave polarization and shear wave splitting are innovatively reproduced via a three-dimensional Finite Element research code upgraded to account for fast dynamic regimes in fully saturated porous media. The mathematical–numerical model adopts a u-v-p formulation enhanced by introducing Taylor–Hood mixed finite elements and the stability features of the solution are considered by analyzing different implemented time integration strategies. Particularly, the phenomena have been studied and reconstructed by numerically generating different types of medium anisotropy accounting for (i) an anisotropic solid skeleton, (ii) an anisotropic permeability tensor, and (iii) a Biot’s effective stress coefficient tensor. Additionally, deviatoric-volumetric coupling effects have been emphasized by specifically modifying the structural anisotropy. A series of analyses are conducted to validate the model and prove the effectiveness of the results, from the directionality of polarized vibrations, the anisotropy-induced splitting, up to the spreading of surface waves.

1. Introduction

The propagation of shear (S) and compressive/pressure (P) waves in poroelastic materials are phenomena relevant for seismology [1,2,3,4,5,6], geotechnical earthquake engineering [7,8,9,10,11], reservoir management [12], and biomechanics of bones and tissues [13]. The shear wave is a polarized transversal wave that propagates within the solid skeleton perpendicularly to the direction of the wave propagation and it often propagates slower than the pressure wave. When a shear wave bumps into an anisotropic medium, it may split into two or more waves with different speeds and orientations, a phenomenon known as shear wave splitting [14]. In particular, if we consider a transversely isotropic medium hit by a shear wave, the wave splits into two orthogonal polarized shear waves propagating at different velocities and orientations according to the material symmetry axis which may not coincide with the initial propagation direction.

The polarization of three-component shear wavetrains carry unique information about the internal structure of the rock: specifically, commonly observed shear-wave splitting may contain information about the orientation of crack distributions.

This information cannot usually be recovered from shear waves recorded at the free surface (i.e., Rayleigh ones) because of interference with the interaction of the shear wave with the surface, even for nearly vertical incidence. However, shear-wave splitting in synthetic three-component vertical seismic profiles in some cases may be interpreted directly in terms of the direction of the strike of vertical cracks [12]. The evolution of such fluid-saturated micro-cracks under changing conditions can be modeled by anisotropic poro-elasticity (APE). Numerical modeling with APE approximately matches a huge range of phenomena, including the evolution of shear-wave splitting during earthquake propagation, and enhanced oil recovery operations [15,16]. APE assumes, and recent observations of shear-wave splitting confirm, that the fluid-saturated cracks in the crust and (probably) upper mantle are so closely spaced that the cracked rocks are highly compliant critical systems with self-organized criticality [17].

Estimating the orientations of cracks, and hence of preferred directions of flow, by seismic investigations could be of crucial importance to production and reservoir engineering [12], as well as evaluating the influence of cracks on the effective elasticity of fractured rocks [18,19]. Additionally, polarization anomalies in seismic shear wavetrains, diagnostic of propagation through anisotropic media, have been observed in dilatancy zones in seismic regions. Stress-induced dilatancy will open cracks with preferred orientations, which will be effectively anisotropic to short-period seismic waves. The polarization anomalies are due to the shear waves splitting, in propagation through anisotropic media, into components with different polarizations and different velocities. This polarization writes characteristic signatures into the shear wavetrains [20].

Different numerical techniques can be used for analyzing wave propagation and shear wave splitting within a porous medium, such as finite difference/finite volume method [21], pseudospectral method [22], and finite element [23]. From the theoretical point of view, pioneer works by Biot are to be recalled [1], whereas in recent years the approaches by [5,10,24,25] can be of reference when analyzing localization and softening effects in wave propagation, or [26] for the contribution of anisotropy.

Here, the triggering and spreading of compressive waves in soils from a point source of a dynamic load are analyzed. Our focus is on the on polarization and shear wave splitting due to the anisotropy of the permeability tensor, the anisotropy of the solid skeleton, as well as to the novel Biot’s tensor that leads to an aniostropic effect on the hydro-mechanical coupling within the effective stress principle. The adopted multi-field dynamics poromechanics [27] Finite Element code is an upgraded dynamic version of a previous static one [28] able to describe the coupled hydro-mechanical behaviour of geomaterials. Particularly, two different solvers have been implemented to perform the time-space integration: the first considers Taylor–Hood elements together with an implicit Euler scheme, the second adopts equal-order elements and a semi-explicit-implicit scheme. A stabilization procedure based on the pressure projection method is used to circumvent the lack of inf-sup condition and resolve the sharp pore pressure gradient [29,30,31,32,33,34]. Comparisons are performed by taking into account literature examples and efficiency features are analyzed. The novel characteristics of the model lie in the synthesis of the compressibility of the fluid phase and the anisotropic permeability, as well as in the extension of the effective stress principle that predicts the anisotropy coupling of the partial stress of solid skeleton and that of pore fluid. One- and two-dimensional analyses have been conducted to validate the model against literature results, whereas a series of three-dimensional simulations are used to show the wave propagation of porous media that exhibits different types and combination of anisotropy.

As for notations and symbols, bold-faced letters denote tensors (including vectors which are rank-one tensors); the symbol ’$·$’ denotes a single contraction of adjacent indices of two tensors (e.g., $\mathit{a}·\mathit{b}=a_i{b}_i$ or $\mathit{c}·\mathit{d}=c_{ij}{d}_{jk}$ ); the symbol ‘:’ denotes a double contraction of adjacent indices of tensor of rank two or higher ( e.g., $\mathit{c}:{ϵ}^{\mathbf{e}}$ = ${\mathit{C}}_{ijkl}{ϵ}_{kl}^e$ ); the symbol ‘⊗’ denotes a juxtaposition of two vectors (e.g., $\mathit{a}\otimes \mathit{b}=a_i{b}_j$) or two symmetric second order tensors (e.g., ${(\mathit{\alpha }\otimes \mathit{\beta })}_{ijkl}={\alpha }_{ij}{\beta }_{kl}$). Moreover, ${(\mathit{\alpha }\oplus \mathit{\beta })}_{ijkl}={\alpha }_{jl}{\beta }_{ik}$ and ${(\mathit{\alpha }\ominus \mathit{\beta })}_{ijkl}={\alpha }_{il}{\beta }_{jk}$. We also define identity tensors ${\left(\mathit{I}\right)}_{ij}={\delta }_{ij}$, ${\left({\mathit{I}}^4\right)}_{ijkl}={\delta }_{ik}{\delta }_{jl}$, and ${\left({\mathit{I}}_{\mathrm{sym}}^4\right)}_{ijkl}=\frac{1}{2}({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{kj})$, where ${\delta }_{ij}$ is the Kronecker delta.

2. Mathematical Model

This section provides a brief review on the field equations for dynamic poromechanics problems, the finite element discretization, and the monolithic and semi-implicit scheme used to simulate the wave propagation in porous media (Section 2.1). This review is followed by the constitutive laws that replicate the anisotropy of the solid skeleton elasticity, and those of the hydraulic responses (permeability), of hydromechanical coupling mechanisms (tensorial Biot’s coefficient).

2.1. Governing Equations

The porous media are here modeled following the mixture theory where the existence of an effective medium of a size suitable for the Representative Elementary Volume (REV) is assumed. A fully saturated porous material is then represented by a corresponding effective medium in which the specific distributions of the fluid and solid constituents inside the REV are homogenized such that the volume of each material point is occupied by a fraction of solid and fluid constituents (Figure 1). The balance of mass and linear momentum equations in the dynamic regime read

$$\begin{array}{c}\hfill {\dot{\left({\rho }^{\alpha }\right)}}_{\alpha }+{\rho }^{\alpha }{\nabla }^x·v_{\alpha }=0;\end{array}$$

(1)

$$\begin{array}{c}\hfill {\rho }^{\alpha }{\dot{\left(v_{\alpha }\right)}}_{\alpha }={\nabla }^x·{\sigma }^{\alpha }+{\rho }^{\alpha }g+h^{\alpha },\end{array}$$

(2)

where: ${\rho }^{\alpha }=n^{\alpha }{\rho }_{\alpha }$ is the partial mass density of $\alpha$-phase (i.e., solid S or fluid F) computed as a product of the volume fraction $n_{\alpha }$ and the intrinsic mass density of $\alpha$ constituent (The Stanford notation is adopted here to indicate the partial quantities (superscript index) and intrinsic quantities (subscript index))—with: ${\sum }_{\alpha }{\rho }^{\alpha }=\rho$. $v_{\alpha }$ being the velocity vector, ${\sigma }^{\alpha }$ the partial Cauchy stress tensor of $\alpha$ phase, $g$ is the gravitational acceleration, and $h^{\alpha }$ is the volume specific local interaction force between the phases, so that: ${\sum }_{\alpha }{h}^{\alpha }=0$. Furthermore, the spatial gradient operator error and the material time derivative ${\dot{(*)}}_{\alpha }$ are referred to as:

$$\begin{array}{c}\hfill {\nabla }^x(*)=\frac{\partial (*)}{\partial x};\end{array}$$

(3)

$$\begin{array}{c}\hfill {\dot{(*)}}_{\alpha }=\frac{d_{\alpha }(*)}{dt}=\frac{\partial (*)}{\partial t}+{\nabla }^x(*)·v_{\alpha },\end{array}$$

(4)

No thermal effects nor mass exchanges between phases are taken into account; the bulk modulus of the $\alpha$-phase can be defined as:

$$K_{\alpha }={\rho }_{\alpha }\frac{d{p}_{\alpha }}{d{\rho }_{\alpha }}$$

(5)

where $p_{\alpha }$ is the intrinsic Cauchy pressure of the $\alpha$-phase.

2.2. Constitutive Laws

2.2.1. Effective Stress Law for Anisotropic Elasticity

It is assumed to consider a homogeneous, linear elastic, anisotropic porous medium so that the stress–strain law can be written as:

$$\sigma =C:ϵ$$

(6)

where $\mathit{c}$ is the fourth-order elastic constitutive tensor of the mixture and its components depend on the material symmetries; when adopting Voigt notation:

$$\overline{\sigma }={\left[\begin{array}{cccccc}{\sigma }_{xx}& {\sigma }_{yy}& {\sigma }_{zz}& {\sigma }_{xy}& {\sigma }_{yz}& {\sigma }_{zx}\end{array}\right]}^T;\overline{ϵ}={\left[\begin{array}{cccccc}{ϵ}_{xx}& {ϵ}_{yy}& {ϵ}_{zz}& {ϵ}_{xy}& {ϵ}_{yz}& {ϵ}_{zx}\end{array}\right]}^T;$$

(7)

a transversely isotropic material is characterized by the 6 × 6 compliance matrix:

$$\overline{\mathit{S}}=\left[\begin{array}{cccccc}1/E_h& -{\nu }_{hh}/E_h& -{\nu }_{vh}/E_v& 0& 0& 0\\ -{\nu }_{hh}/E_h& 1/E_h& -{\nu }_{vh}/E_v& 0& 0& 0\\ -{\nu }_{hv}/E_h& -{\nu }_{hv}/E_h& 1/E_v& 0& 0& 0\\ 0& 0& 0& 1/G_{hh}& 0& 0\\ 0& 0& 0& 0& 1/G_{hv}& 0\\ 0& 0& 0& 0& 0& 1/G_{vh}\end{array}\right]$$

(8)

The model exhibits symmetry along any direction h belonging to the isotropy plane and different material properties along the normal direction v. The independent constants are: $E_h,E_v,G_{hv},{\nu }_{hh},{\nu }_{hv}$, being the other terms given by:

$${\nu }_{hv}/E_h={\nu }_{vh}/E_v;G_{hh}=E_{hh}/\left(2(1+{\nu }_{hh})\right).$$

(9)

Generally, the symmetry axes for an anisotropic material with respect to the reference system are not coaxial with the global axes so that the elastic matrix has all non-zero elements; hence, to express stress and strain measures into the reference system, the $6\times 6$ transformation matrix is adopted:

$$\overline{\mathit{M}}=\left[\begin{array}{cccccc}{r}_{11}^2& r_{21}^2& r_{31}^2& 2r_{11}{r}_{21}& 2r_{21}{r}_{31}& 2r_{11}{r}_{31}\\ r_{12}^2& r_{22}^2& r_{32}^2& 2r_{12}{r}_{22}& 2r_{22}{r}_{32}& 2r_{12}{r}_{32}\\ r_{13}^2& r_{23}^2& r_{33}^2& 2r_{13}{r}_{23}& 2r_{23}{r}_{33}& 2r_{13}{r}_{33}\\ r_{11}{r}_{12}& r_{21}{r}_{22}& r_{31}{r}_{32}& r_{11}{r}_{22}+r_{21}{r}_{12}& r_{21}{r}_{32}+r_{31}{r}_{22}& r_{11}{r}_{32}+r_{31}{r}_{12}\\ r_{12}{r}_{13}& r_{22}{r}_{23}& r_{32}{r}_{33}& r_{12}{r}_{23}+r_{22}{r}_{13}& r_{22}{r}_{33}+r_{32}{r}_{23}& r_{12}{r}_{33}+r_{32}{r}_{13}\\ r_{11}{r}_{13}& r_{21}{r}_{23}& r_{31}{r}_{33}& r_{11}{r}_{23}+r_{21}{r}_{13}& r_{21}{r}_{33}+r_{31}{r}_{23}& r_{11}{r}_{33}+r_{31}{r}_{13}\end{array}\right];\mathit{R}=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]$$

(10)

so that:

$${\overline{\mathit{C}}}_{glob}=\overline{\mathit{M}}{\overline{\mathit{C}}}_{loc}{\overline{\mathit{M}}}^T,or{\overline{\mathit{S}}}_{glob}={\overline{\mathit{M}}}^{-T}{\overline{\mathit{S}}}_{loc}{\overline{\mathit{M}}}^{-1};$$

(11)

with the bar above the symbols indicating the matrix version of the tensor quantities.

2.2.2. Tensorial Nature of Biot’s Coefficient

Here, the effective stress principle is adopted, where the Biot’s effective stress coefficient is not a scalar but a second-order tensor with different eigenvalues (see Cowin and Doty [26], Carroll [35]):

$${\sigma }_E^S=\sigma +\mathit{A}{p}_F;\mathit{A}=\left(\mathit{I}-\mathit{C}{\mathit{S}}_S\right):\mathbf{1},$$

(12)

where ${\sigma }_E^S$ is the effective stress tensor, $\sigma ={\sigma }^S+{\sigma }^F$ is the total stress tensor sum of the two partial stresses, $p_F$ is the Cauchy pore pressure, $\mathbf{1}$ and $\mathit{I}$ are the unit and the identity tensors, respectively, whereas ${\mathit{S}}_s$ refers to the solid skeleton. Equation (12) represents an extension of the classical effective stress concept, based on the assumption of isotropy for both porous material and solid phase. The expression comes from analyzing the constitutive stress–strain behaviour of the REV by assuming structural (anisotropic porous geometry), intrinsic (anisotropic solid material) or both anisotropies. As described in Cowin and Doty [26], the Biot’s effective stress tensor could have six nonzero components if the material symmetry of compliance tensor ${\mathit{S}}_S$ is less than transversely isotropic and/or its axis of symmetry are not coaxial with axis of the transversely isotropic model of the porous material $\mathit{C}$. If both porous material and its solid skeleton are isotropic, we obtain that:

$$\mathit{A}=B\mathbf{1};B=\left[1-(K/K_S)\right]$$

(13)

where B is the classic Biot’s effective stress coefficient [36], equal to one if incompressibility of the solid skeleton is assumed.

2.2.3. Darcy’s Law

By assuming laminar flow, the generalized form of Darcy’s law is considered, relating the fluid mass flow rate to the pore pressure gradient:

$$n^F{w}_F=-{\mathit{K}}^F({\nabla }^x{p}_F-{\rho }_{F}g)$$

(14)

where $w_F=v_F-w_S$ is the relative velocity and $w^F=n^F{w}_F$ the Darcy velocity. The quantity ${\mathit{K}}^F$ is the second-order specific permeability tensor:

$${\mathit{K}}^F={\mathit{k}}^F/{\gamma }_F={\widehat{\mathit{K}}}^F/{\mu }_F$$

(15)

where ${\gamma }_F={\rho }_{F}g$ is the fluid intrinsic unit weight, ${\mu }_F$ is the fluid viscosity, ${\mathit{k}}^F$ is the hydraulic conductivity tensor, and ${\widehat{\mathit{K}}}^F$ is the intrinsic permeability. The latter is a second-order tensor of which the eigenvalues and spectral directions can be determined from experimental filtration tests or through 3D tomography images and geometrical analyses as used in Sun et al. [37].

3. Numerical Implementation

3.1. Galerkin Form

The subsequent set of coupled partial differential equations (PDEs) is hence obtained:

$$\begin{array}{c}\hfill {\dot{\left({\mathit{u}}_S\right)}}_S=v_S;\end{array}$$

(16)

$$\begin{array}{c}\hfill {\rho }^S{\dot{\left(v_S\right)}}_S={\nabla }^x\left[{\sigma }_E^S-\left(\mathit{A}-n^F\mathbf{1}\right)p_F\right]+{\rho }^{S}g+\frac{{\left(n^F\right)}^2{\gamma }_F}{{\mathit{K}}^F}{w}_F-p_F{\nabla }^x\left(n^F\right);\end{array}$$

(17)

$$\begin{array}{c}\hfill {\rho }^F{\dot{\left(v_F\right)}}_S+{\rho }_F{\nabla }^x\left(v_F\right)w_F=-n^F{\nabla }^x\left(p_F\right)+{\rho }^{F}g-\frac{{\left(n^F\right)}^2{\gamma }_F}{{\mathit{K}}^F}{w}_F;\end{array}$$

(18)

$$\begin{array}{c}\hfill \frac{n^S}{K_S}{\dot{\left(p_S\right)}}_S+{\nabla }^x{v}_S+\frac{n^F}{K_F}{\dot{\left(p_F\right)}}_S+\frac{1}{{\rho }_F}{\nabla }^x\left({\rho }^F{w}_F\right)=0.\end{array}$$

(19)

where the unknowns of the coupled problem are $u_S$, $v_S$, $v_F$, and $p_F$. This set of equations represents the three-field formulation (it is recalled that $v_S$ are secondary unknowns and the first equation is written to compute them). Furthermore, if we consider that the porous material is subjected to low and relatively small frequencies ($\le$30 Hz in the geomechanical problems), we can assume that ${\dot{\left(u_S\right)}}_S\approx {\dot{\left(u_F\right)}}_S$, and ${\dot{\left(w_F\right)}}_S\approx 0$ so we can rearrange the equations and obtain the classical two-field formulation for consolidation problems with $u_S$ and $p_F$. The governing equations become:

$$\begin{array}{c}\hfill {\dot{\left(u_S\right)}}_S=v_S;\end{array}$$

(20)

$$\begin{array}{c}\hfill \rho {\dot{\left(v_S\right)}}_S={\nabla }^x\left({\sigma }_E^S-\mathit{A}{p}_F\right)+\rho g;\end{array}$$

(21)

$$\begin{array}{c}\hfill \frac{n^S}{K_S}{\dot{\left(p_S\right)}}_S+{\nabla }^x{v}_S+\frac{n^F}{K_F}{\dot{\left(p_F\right)}}_S+\frac{1}{{\rho }_F}{\nabla }^x\left({\rho }^F{w}_F\right)=0.\end{array}$$

(22)

By considering a finite domain $\mathsf{\Omega }$ of the mixture with its boundary $\mathsf{\Gamma }$ and a set of independent test functions, the weak formulation for the previous PDEs can be written as:

$$\begin{array}{cc}\hfill {\mathcal{G}}_{v_S}\left({\mathit{u}}_S,v_S,\delta {\mathit{u}}_S\right)=& {\int }_{\mathsf{\Omega }}\delta {\mathit{u}}_S·\left[{\dot{\left({\mathit{u}}_S\right)}}_S-v_S\right]dv=0;\hfill \\ \hfill {\mathcal{G}}_{{\mathit{u}}_S}\left({\mathit{u}}_S,v_S,v_F,p_F,\delta {\mathit{u}}_S\right)=& {\int }_{\mathsf{\Omega }}{\nabla }^x\left(\delta {\mathit{u}}_S\right)·\left[{\sigma }_E^S-\left(\mathit{A}-n^F\mathit{I}\right)p_F\right]dv+{\int }_{\mathsf{\Omega }}\delta {\mathit{u}}_S·\left[-\frac{{\left(n^F\right)}^2{\gamma }_F}{{\mathit{K}}^F}{w}_F+p_F{\nabla }^x\left(n^F\right)\right]dv\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & +{\int }_{\mathsf{\Omega }}\delta {\mathit{u}}_S·\left\{{\rho }^S\left[{\dot{\left(v_S\right)}}_S-g\right]\right\}dv-{\int }_{\partial {\mathsf{\Omega }}_S}\delta {\mathit{u}}_S·{\mathit{t}}^{S}da=0;\hfill \\ \hfill {\mathcal{G}}_{v_F}\left({\mathit{u}}_S,v_S,v_F,p_F,\delta v_F\right)=& -{\int }_{\mathsf{\Omega }}{\nabla }^x\left(\delta v_F\right)n^F{p}_{F}dv-{\int }_{\partial {\mathsf{\Omega }}_F}\delta v_F·{\mathit{t}}^{F}da\hfill \end{array}$$

(24)

$$\begin{array}{c}\hfill +{\int }_{\mathsf{\Omega }}\delta v_F·\left\{{\rho }_F\left[{\dot{\left(v_F\right)}}_S-g+\left({\nabla }^x{v}_F+\frac{n^F{\gamma }_F}{{\mathit{K}}^F}\right)w_F\right]\right\}dv=0;\end{array}$$

(25)

$$\begin{array}{cc}\hfill {\mathcal{H}}_{p_F}\left(v_S,v_F,p_F,\delta p\right)=& {\int }_{\mathsf{\Omega }}\delta p\mathit{A}:{\nabla }^x{v}_{S}dv+{\int }_{\mathsf{\Omega }}\delta p\mathsf{\Lambda }{\dot{\left(p_F\right)}}_{S}dv-{\int }_{\mathsf{\Omega }}{\nabla }^x\left(\delta p\right)·\left(n^F{w}_F\right)dv+{\int }_{\partial {\mathsf{\Omega }}_p}\delta p·\overline{v}da=0.\hfill \end{array}$$

(26)

where $\delta {\mathit{u}}_S$, $\delta {\mathit{u}}_F$, $\delta v_F$, and $\delta p$ are the test functions used as weight. The quantities ${\mathit{t}}^S$ and ${\mathit{t}}^F$ are the surface tractions of the single phases acting on the Neumann boundary and $\overline{v}$ is the volume flux of the fluid going through the boundary of the mixture. The compressibility modulus $\mathsf{\Lambda }$ is computed as:

$$\mathsf{\Lambda }=\frac{n^S}{K_S}+\frac{n^F}{K_F}-\mathit{A}:\mathit{C}:\mathit{A}=\frac{n^S}{K_S}+\frac{n^F}{K_F}-\mathbf{1}:{\mathit{C}}_S\mathit{S}{\mathit{c}}_S:\mathbf{1};$$

(27)

this expression is the same as Equation (8.5) in Cowin [38].

In order to write the Galerkin formulation in a compact way, Equations (23)–(26) and all the field variables are collected within vectors:

$${\mathcal{G}}_{\mathit{u}}=\left[\begin{array}{c}{\mathcal{G}}_{v_S}\\ {\mathcal{G}}_{{\mathit{u}}_S}\\ {\mathcal{G}}_{v_F}\\ {\mathcal{G}}_{{\mathit{p}}_F}\end{array}\right],\mathit{u}=\left[\begin{array}{c}{\mathit{u}}_S\\ v_S\\ v_F\\ {\mathit{p}}_F\end{array}\right],\delta \mathit{u}=\left[\begin{array}{c}\delta {\mathit{u}}_S\\ \delta {\mathit{u}}_S\\ \delta v_F\\ \delta {\mathit{p}}_F\end{array}\right],{\dot{\left(\mathit{u}\right)}}_S=\left[\begin{array}{c}{\dot{\left({\mathit{u}}_S\right)}}_S\\ {\dot{\left(v_S\right)}}_S\\ {\dot{\left(v_F\right)}}_S\\ {\dot{\left({\mathit{p}}_F\right)}}_S\end{array}\right],{\mathit{u}}_0=\left[\begin{array}{c}{\mathit{u}}_{S0}\\ v_{S0}\\ v_{F0}\\ {\mathit{p}}_{F0}\end{array}\right].$$

(28)

The spatial discretization is carried out via the Finite Element Method, referring to the following discrete unknown variables and test functions, respectively:

$$u^h(x,t)={\overline{u}}^h+\sum _{i=1}^{N_u}{N}_{ui}\left(x\right)u_i\left(t\right),\delta u^h\left(x\right)=\sum _{i=1}^{M_u}{M}_{ui}\left(x\right)\delta u_i\left(t\right);$$

(29)

where ${\overline{u}}^h$ represents the solution on the Dirichlet boundary, $N_u$ and $M_u$ are the numbers of nodes, and $N_{ui}$ and $M_{ui}$ indicate the shape functions at node i depending only on the position $x$. The final three-field variational problem can be re-written via the Petrov–Galerkin formulation ${\mathcal{G}}_u(\delta u^h,u^h)=0$ using the previous test and trial functions.

3.2. Matrix Form and Time Discretization

By omitting the superscript h for vectors $\mathit{u}$ and $\delta \mathit{u}$ and indicating with ${\mathit{f}}_S$, ${\mathit{f}}_F$ and ${\mathit{f}}_p$ the space-discrete Neumann boundary terms, the subsequent matrix formulation can be obtained:

$${\mathcal{G}}_u^h=\left[\begin{array}{cccc}\mathit{I}& 0& 0& 0\\ 0& {\mathit{M}}_{22}& 0& 0\\ 0& 0& {\mathit{M}}_{33}& 0\\ 0& 0& 0& {\mathit{M}}_{44}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{u}}}_S\\ {\dot{v}}_S\\ {\dot{v}}_F\\ {\dot{\mathit{p}}}_F\end{array}\right]+\left[\begin{array}{cccc}0& -\mathit{I}& 0& 0\\ {\mathit{K}}_{21}& {\mathit{K}}_{22}& {\mathit{K}}_{23}& {\mathit{K}}_{24}\\ 0& {\mathit{K}}_{32}& {\mathit{K}}_{33}& {\mathit{K}}_{34}\\ 0& {\mathit{K}}_{42}& {\mathit{K}}_{43}& 0\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_S\\ v_S\\ v_F\\ {\mathit{p}}_F\end{array}\right]-\left[\begin{array}{c}0\\ {\mathit{a}}_S+{\mathit{f}}_S\\ {\mathit{a}}_F+{\mathit{f}}_F\\ {\mathit{f}}_p\end{array}\right]=0$$

(30)

The system (30) is composed by a set of pure differential equations where the time evolution for all the primary variables can be written as a first-order ordinary differential equation:

$${\mathcal{G}}_{\mathit{u}}^h=\mathit{M}\dot{\mathit{u}}+\mathit{K}\mathit{y}-\mathit{f}=0;\dot{\mathit{u}}=g(t,\mathit{u}),and\mathit{y}\left(t_0\right)={\mathit{y}}_0$$

(31)

This formulation can be solved via implicit time integration schemes or explicit ones. In the latter case, the stability constraint must be taken into consideration when choosing the time step size. In the case where both solid and fluid constituents are incompressible, the system (30) becomes a differential algebraic one, where the fourth relation turns into a volume balance equation for the mixture acting as an algebraic incompressibility constraint. In order to solve a saddle point problem, the time integration schemes must be upgraded by special numerical techniques to maintain the stability of the numerical solution. Following Markert et al. [39], two numerical time integration strategies have been adopted: the first one is an implicit monolithic scheme together with mixed order interpolation for the primary variables; the second one is a semi-explicit/implicit splitting scheme together with an equal order interpolation for the unknowns.

3.2.1. Implicit Monolithic Schemes

Different types of implicit schemes exist and it is possible to subdivide them mainly into two classes: one step and multi-step methods. Thanks to the lower computational cost and memory storage of the variables, only the first class is considered here. In Algorithm 1, the pseudo code of the generalized trapezoidal scheme inside a Newton–Raphson procedure is shown. By considering the integration variable $\theta =1$, the Backward Euler scheme is obtained; this method is first order accurate but small time steps are mandatory in order to reduce the artificial numerical damping leading to wrong solution (see Jansen et al. [40]). For $\theta =1/2$, the second-order accurate Crank–Nicholson scheme comes out, and with $\theta =0$ the Explicit Forward Euler scheme appears.

Note that the last equation of the coupled system (30) is modified by inserting the expression of Darcy’s velocity derived from the fluid momentum balance in order to improve the numerical stability of the scheme.

Algorithm 1: Newton–Raphson Algorithm.

Initialization: ${\mathit{u}}_n={\left[{\mathit{u}}_{Sn}{v}_{Sn}{v}_{Fn}{\mathit{p}}_n\right]}^T=0,{\mathit{F}}_n={\left[0{\mathit{f}}_{Sn}{\mathit{f}}_{Fn}{\mathit{q}}_n\right]}^T=0$

for $n=1:n_{end}$

${\mathit{F}}_{n+1}={\mathit{F}}_n+\theta \Delta \mathit{F},{\mathit{u}}_{n+1}={\mathit{u}}_n,$

for $i=1:i_{max}$

$\Delta {\mathit{u}}_{n+1}^i={\mathit{u}}_{n+1}^i-{\mathit{u}}_n$

Compute matrix:${\mathit{M}}_{ii},{\mathit{K}}_{ii}$

Compute residual:

${\mathit{R}}^{i+1}=\left[\begin{array}{c}0\\ {\mathit{r}}_{v_S}\\ {\mathit{r}}_{v_F}\\ {\mathit{r}}_{p_F}\end{array}\right]={\mathit{F}}_{n+1}+\left[\begin{array}{c}0\\ {\mathit{a}}_2\\ {\mathit{a}}_3\\ {\mathit{a}}_4\end{array}\right]-\left[\begin{array}{cccc}\mathit{I}& 0& 0& 0\\ 0& {\mathit{M}}_{22}& 0& 0\\ 0& 0& {\mathit{M}}_{33}& 0\\ 0& 0& {\mathit{M}}_{43}& {\mathit{M}}_{44}\end{array}\right]\frac{\Delta {\mathit{u}}_{n+1}^i}{\Delta t}-\left[\begin{array}{cccc}0& -\mathit{I}& 0& 0\\ {\mathit{K}}_{21}& {\mathit{K}}_{22}& {\mathit{K}}_{23}& {\mathit{K}}_{24}\\ 0& {\mathit{K}}_{32}& {\mathit{K}}_{33}& {\mathit{K}}_{34}\\ 0& {\mathit{K}}_{42}& {\mathit{K}}_{43}& {\mathit{K}}_{44}\end{array}\right]\left[\begin{array}{c}\theta {\mathit{u}}_{n+1}^i+(1-\theta ){\mathit{u}}_n\end{array}\right]$

Check residual:

if${\mathit{R}}^{i+1}<toll$

break

end

Compute Jacobi matrix:

$\mathit{J}=\mathit{M}/\Delta t+\theta \mathit{K}$

Compute Y-increment:

$d\mathit{u}={\mathit{J}}^{-1}{\mathit{R}}^{i+1}$

Update solution:

${\mathit{u}}_{n+1}^{i+1}={\mathit{u}}_{n+1}^i+d\mathit{u}$

end

end

3.2.2. Semi-Explicit/Implicit Splitting Scheme

This scheme solves the coupled problem using a splitting procedure applied in the field of porous mechanics by Huang et al. [41,42]. Through this scheme, the system of equations is split into an implicit and a subsequent explicit step. This method is restricted by a critical time-step to guarantee stability and accuracy of the solution. The splitting separates the linear momentum balances from the mass balance of the mixture and decouples the displacement and velocity fields from the pore-fluid pressure field. To this aim, the time discretization of the governing equations is to be developed before the spatial one together with the definition of an intermediate velocity giving an approximation of the velocities of the phases in the next time step. By following the same procedure of Markert et al. [39], the time discretization and splitting procedure starts with the implicit discretization of the solid velocity. Equation (16) becomes through the trapezoidal rule:

$$\frac{(u_{Sn+1}-u_{Sn})}{\Delta t_n}=\frac{1}{2}(v_{Sn+1}+v_{Sn});$$

(32)

assuming that $n^{\alpha }\approx n_0^{\alpha }$, ${\nabla }^x{n}^{\alpha }\approx 0$ (small strain assumption); by explicitly considering the solid extra stress tensor ${\sigma }_{En}^S={\sigma }_E^S\left(u_{Sn}\right)$ and implicitly the pore-fluid pressure term $p_F$ together with the relative velocities through the intermediate velocities $v_S^{*}$ and $v_F^{*}$, the momentum balance equation for solid phase is rearranged and split as follows:

$$\begin{array}{c}\hfill {\rho }^S\frac{(v_S^{*}-v_{Sn})}{\Delta t_n}={\nabla }^x{\sigma }_{En}^S-\left(\mathit{A}-n^F\mathbf{1}\right){\nabla }^x\left(p_{Fn}\right)+{\rho }^{S}g+\frac{{\left(n^F\right)}^2{\gamma }_F}{{\mathit{K}}^F}{w}_F^{*};\end{array}$$

(33)

$$\begin{array}{c}\hfill {\rho }^S\frac{(v_{Sn+1}-v_S^{*})}{\Delta t_n}=-\left(\mathit{A}-n^F\mathbf{1}\right){\nabla }^x(p_{Fn+1}-p_{Fn}).\end{array}$$

(34)

The momentum balance for the fluid becomes:

$$\begin{array}{c}\hfill {\rho }^F\frac{(v_F^{*}-v_{Fn})}{\Delta t_n}+{\rho }^F\left({\nabla }^x{v}_F^{*}\right)w_F^{*}=-n^F{\nabla }^x\left(p_{Fn}\right)+{\rho }^{F}g-\frac{{\left(n^F\right)}^2{\gamma }_F}{{\mathit{K}}^F}{w}_F^{*};\end{array}$$

(35)

$$\begin{array}{c}\hfill {\rho }^F\frac{(v_{Fn+1}-v_F^{*})}{\Delta t_n}=-n^F{\nabla }^x(p_{Fn+1}-p_{Fn});\end{array}$$

(36)

where ${\rho }^F\left({\nabla }^x{v}_F^{*}\right)w_F^{*}$ is a convective term. The mass balance is expressed as:

$$\mathsf{\Lambda }\frac{p_{Fn+1}-p_{Fn}}{\Delta t}+\mathit{A}:{\nabla }^x\left(v_{Sn+1}\right)+\frac{1}{{\rho }_F}{\nabla }^x\left({\rho }^F{w}_{Fn+1}\right)=0.$$

(37)

Rebuilding the weak formulation, for the solid phase:

$$\begin{array}{cc}\hfill {\int }_{\mathsf{\Omega }}\delta {\mathit{u}}_S·\left[\frac{(v_S^{*}-v_{Sn})}{\Delta t_n}-g\right]n^S{\rho }_{S}dv-{\int }_{\mathsf{\Omega }}\delta {\mathit{u}}_S·\frac{{\left(n^F\right)}^2{\gamma }_F}{{\mathit{K}}^F}{w}_F^{*}dv& \hfill \end{array}$$

$$\begin{array}{cc}\hfill +{\int }_{\mathsf{\Omega }}{\nabla }^x\left(\delta u_S\right)·\left[{\sigma }_{En}^S-\left(\mathit{A}-n^F\mathbf{1}\right)p_{Fn}\right]dv-{\int }_{\partial {\mathsf{\Omega }}_S}\delta u_S·{\mathit{t}}_n^{S}da=0& ;\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {\int }_{\mathsf{\Omega }}\delta {\mathit{u}}_S·\left[{\rho }^S\frac{(v_{Sn+1}-v_S^{*})}{\Delta t_n}+\left(\mathit{A}-n^F\mathbf{1}\right){\nabla }^x(p_{Fn+1}-p_{Fn})\right]dv=0& ;\hfill \end{array}$$

(39)

for the fluid phase:

$$\begin{array}{cc}\hfill {\int }_{\mathsf{\Omega }}\delta v_F·\left[\frac{(v_F^{*}-v_{Fn})}{\Delta t_n}-g\right]n^F{\rho }_{F}dv+{\int }_{\mathsf{\Omega }}\delta v_F·n^F{\nabla }^x\left(p_{Fn}\right)dv+{\int }_{\mathsf{\Omega }}\delta v_F·\frac{n^F{\gamma }_F}{{\mathit{K}}^F}{w}_F^{*}dv=0& ;\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {\int }_{\mathsf{\Omega }}\delta v_F·\frac{(v_{Fn+1}-v_F^{*})}{\Delta t_n}{n}^F{\rho }_{F}dv-{\int }_{\mathsf{\Omega }}{\nabla }^x\left(\delta v_F\right)n^F(p_{Fn+1}-p_{Fn})dv-{\int }_{\partial {\mathsf{\Omega }}_F}\delta v_F·{\mathit{t}}_{n+1}^{F}da=0& ;\hfill \end{array}$$

(41)

with the mass balance:

$${\int }_{\mathsf{\Omega }}\delta p\mathit{A}:{\nabla }^x{v}_S^{*}dv+{\int }_{\mathsf{\Omega }}\delta p\mathsf{\Lambda }\frac{(p_{Fn+1}-p_{Fn})}{\Delta t_n}dv-{\int }_{\mathsf{\Omega }}{\nabla }^x\left(\delta p\right)·\left(n^F{w}_F\right)dv+{\int }_{\partial {\mathsf{\Omega }}_p}\delta p·{\overline{v}}_{n+1}da=0.$$

(42)

From the weak formulation the corresponding matrix equations are built following the scheme of Algorithm 2 and the expressions of all the matrices are plotted in Appendix A. After initialization, the procedure starts by explicitly computing the intermediate velocities $v^{*}$; then, the pore-fluid pressure is calculated and subsequently the velocity corrections and finally the solid displacements $u_{Sn+1}$. Clearly, each time-step requires an iteration check to be satisfied: $∥u_{n+1}^{i+1}-u_{n+1}^i∥<toll$, with the critical time-step given by (for a 3D linear FE):

$$\Delta t_{cr}=\frac{\Delta h_x\Delta h_y\Delta h_z}{c_p^x\Delta h_x+c_p^y\Delta h_y+c_p^z\Delta h_z}{c}_p^r=\sqrt{\frac{C_{rr}}{{\rho }^S}},r=x,y,z,$$

(43)

with $C_{rr}$ elastic matrix component of the porous material.

Algorithm 2: Prediction/Correction Algorithm.

Initialization: ${\mathit{u}}_{Sn}=0,v_n=0,{\mathit{p}}_n=0$

for $n=1:n_{end}$

${\mathit{f}}_n={\mathit{f}}_n+\Delta \mathit{f},v^P=v_n,{\mathit{p}}^P={\mathit{p}}_n,$

$\Delta {\mathit{u}}_S=\Delta t{v}_S^P,{\mathit{u}}_{Sn}+\Delta {\mathit{u}}_S$

for $i=1:i_{max}$

Compute matrix:$\mathit{M},{\mathit{K}}_{ii},{\overline{\mathit{K}}}_{ii}$

Compute prediction velocities:

$\left[\begin{array}{c}{v}_S^{*}\\ v_F^{*}\end{array}\right]={\left[\begin{array}{cc}\frac{{\mathit{M}}_{22}}{\Delta t}+{\mathit{K}}_{22}& -{\mathit{K}}_{23}\\ -{\mathit{K}}_{32}& \frac{{\mathit{M}}_{33}}{\Delta t}+{\mathit{K}}_{33}\end{array}\right]}^{-1}\left\{\begin{array}{c}\left[\begin{array}{cc}\frac{{\mathit{M}}_{22}}{\Delta t}& 0\\ 0& \frac{{\mathit{M}}_{33}}{\Delta t}\end{array}\right]\left[\begin{array}{c}{v}_{Sn}^P\\ v_{Fn}^P\end{array}\right]+\left[\begin{array}{cc}-{\mathit{K}}_{21}& {\mathit{K}}_{24}\\ 0& -{\overline{\mathit{K}}}_{34}\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_{Sn}^P\\ {\mathit{p}}_n^P\end{array}\right]+\left[\begin{array}{c}{\mathit{a}}_2\\ {\mathit{a}}_3\end{array}\right]+\left[\begin{array}{c}{\mathit{f}}_{Sn}\\ 0\end{array}\right]\end{array}\right\}$

Compute pore fluid pressure:

$\left[\begin{array}{c}0\\ {\mathit{p}}_{n+1}^{*}\end{array}\right]=\left[\begin{array}{c}0\\ {\mathit{p}}_n^P\end{array}\right]+{\left[\begin{array}{cc}\mathit{I}& 0\\ 0& {\overline{\mathit{K}}}_{44}\end{array}\right]}^{-1}\left\{\begin{array}{c}\left[\begin{array}{cc}0& 0\\ -{\mathit{K}}_{42}& {\mathit{K}}_{43}\end{array}\right]\left[\begin{array}{c}{v}_S^{*}\\ v_F^{*}\end{array}\right]-\left[\begin{array}{c}0\\ {\overline{\mathit{f}}}_{Pn+1}\end{array}\right]\end{array}\right\}$

Compute velocities correction:

$\left[\begin{array}{c}{v}_{Sn+1}\\ v_{Fn+1}\end{array}\right]=\left[\begin{array}{c}{v}_S^{*}\\ v_F^{*}\end{array}\right]+{\left[\begin{array}{cc}\frac{{\mathit{M}}_{22}}{\Delta t}& 0\\ 0& \frac{{\mathit{M}}_{33}}{\Delta t}\end{array}\right]}^{-1}\left\{\begin{array}{c}\left[\begin{array}{cc}0& {\overline{\mathit{K}}}_{24}\\ 0& -{\mathit{K}}_{34}\end{array}\right]\left[\begin{array}{c}0\\ {\mathit{p}}_{n+1}-{\mathit{p}}_n^P\end{array}\right]+\left[\begin{array}{c}0\\ {\mathit{f}}_{Pn+1}\end{array}\right]\end{array}\right\}$

Compute solid displacements:

${\mathit{u}}_{Sn+1}={\mathit{u}}_{Sn}+\frac{1}{2}\Delta t\left(v_{Sn+1}+v_{Sn}\right)$

${\mathit{r}}_i={\mathit{u}}_{Sn+1}-{\mathit{u}}^P$

if${\mathit{r}}_i<toll$

break

end

Update prediction variables:

$v^P=v_{n+1},{\mathit{u}}_S^P={\mathit{u}}_{Sn+1},\Delta {\mathit{u}}_S={\mathit{u}}_{Sn+1}-{\mathit{u}}_{Sn},$

end

end

4. Numerical Examples

In this section, some numerical analyses are described, accounting for different types of anisotropy affecting P- and S-wave motion within porous media.

Particularly, waves propagation and interaction are studied, influenced by (1) the induced anisotropy related to volumetric-deviatoric coupling effects and (2) the inherent anisotropy of transversely isotropic elasticity for the solid skeleton, the anisotropy of the permeability tensor, and the anisotropic hydro-mechanical coupling effect captured by the tensorial Biot’s approach. The numerical analyses have been developed via GeoMatFEM [28], a MATLAB research software suitable for coupled geo-mechanical simulations and now upgraded to a dynamic version. Two benchmarks are included to validate the implementation.

4.1. Benchmark Cases with Isotropic Elastic Materials

The code has been validated against two numerical examples:

• fully saturated soil column under harmonic load (cf. de Boer et al. [43]);
• wave propagation within a two-dimensional soil domain (cf. Markert et al. [39]).

The soil column model shown in Figure 2a is considered, with a vertical discretization of 10 Finite Elements/meter subjected to a vertical harmonic load Figure 2b. The material parameters are shown in

Table 1. Material parameters for benchmarks (1) and (2).

Parameter	Values	S.I. unit
E	14.52 × 10${}^6$	Pa
$\nu$	0.30
$n_0^F$	0.33
$k^F$	10${}^{-2}$	m/s
${\rho }_S$	2000	kg/m${}^3$
${\rho }_F$	1000	kg/m${}^3$

.

The results obtained by adopting different solvers together with the analytical solution are depicted in Figure 3. The Implicit Backward Euler scheme B.E. (coupled with mixed elements), the semi-explicit/implicit scheme S.E. (with linear equal order elements), both considering the uvp–formulation and the classical up–formulation), have been taken into account. For the two implicit schemes, a time step $dt=0.5\times {10}^{-3}\mathrm{s}$ has been adopted, while, for the semi-explicit/implicit scheme, a $dt=2.5\times {10}^{-4}\mathrm{s}$ has been assumed. All of the numerical solutions give equal results in terms of compaction (see Figure 3a) and pore pressure (see Figure 3c); slight differences are observable in the peak values of the fluid velocity on surface (see Figure 3b), but this is due to the choice of the exit error tolerance of the schemes, while, for the effective stress (see Figure 3d), the differences come from the fact that they are computed at Gauss points $0.0225$ m away from the reference surface.

As regards Benchmark (2), the F.E. domain together with the boundary conditions are shown in Figure 4a. The material parameters are the same as in the previous example (see Table 1) and the soil is subjected to an impulsive load Figure 4b, with $H(t-\tau )$ Heaviside function and $\tau =0.04s$ the duration of the impulse.

Figure 5 reports the main results by varying the numerical solver, i.e., the implicit Backward Euler scheme with mixed elements, considering both uv and uvp–formulation, the Semi-explicit/Implicit scheme with linear (L) and quadratic (Q) equal order elements. A composite time integration scheme between the trapezoidal rule and 2nd–order backward difference scheme (TR-BDF2) has been additionally adopted: this scheme is inserted in the Runge–Kutta method (s-stage DIRK) together with mixed elements. In Markert et al. [39], this scheme was applied for the first time in the field of multiphase porous mechanics, the propriety of such a method is that it satisfies all the stability requirements and it is second order accurate: for these reasons, its results have been taken as a reference solution. The same mesh discretization composed by 4 hexahedral elements per square meter is assumed for all the models; furthermore, for the three implicit schemes and for the semi-explicit one with linear elements, the time step $dt$ = 10⁻³ s has been adopted, while for the semi-explicit scheme with quadratic elements, $dt$ = $2.5$ 10⁻⁴ s has been assumed in order to satisfy the CFL stability condition.

Figure 5a shows the vertical displacements of node C; all the numerical solutions coincide, even considering pore pressure at node B (see Figure 5d) except for the lowest peak value. Figure 5c describes the quasi-elliptical (or “eight type”) motion of node A due to the Rayleigh wave generated by the impulsive load. As confirmed by Markert et al. [39], the uvp–formulation with a Backward Euler scheme provides stiff results due to the fact that this type of scheme possesses a strong artificial damping (see Jansen et al. [40]), while the up–formulation overestimates the displacements field due to the assumption of zero relative acceleration between the two phases. The results of the semi-explicit scheme are closer to the solution of TR-BDF2 scheme and the one with quadratic interpolation appears to be the best.

Figure 6 depicts the time sequence of displacements contour and deformed mesh for the semi-explicit scheme with quadratic finite elements. The slow pressure wave (P, generating radial compression) and the shear wave (S, with shear type deformation) propagating in the soil domain are accompanied by the Rayleigh wave (R–wave) moving at the surface of the medium.

4.2. Dynamic Poroelastic Responses of Isotropic Porous Media

3D analyses are needed to evaluate the seismic waves in the soil motion and thus to evaluate the predictive capabilities of the F.E. code. For comparison, we first simulate the wave propagation in a fully saturated, homogeneous, and isotropic poroelastic material and analyze the results. Only the semi-explicit-implicit procedure is used due to the robustness of this scheme to solve dynamic problems with a suitable time step size. Furthermore, the splitting procedure provides for a faster solver and requires less computational effort to solve the system of equations.

A solid square prism of soil (dotted lines of Figure 7) subjected to an impulsive pressure load applied on an area of 1 m2 on the top surface has been considered. Taking advantage of the symmetry of the problem, one quarter of the prism has been modeled only; bottom and lateral surfaces are assumed impermeable, frictionless and restrained along the normal direction, whereas the top surface permeable. The FE model is composed by 8820, $21\times 21\times 20$ 3D, linear and equal order elements, with material parameters listed in

Table 2. Material parameters for the 3D model.

Parameter	Values	S.I. Unit
E	12.0 × 10${}^6$	Pa
$\nu$	0.25
$n_0^F$	0.33
$k^F$	10.0${}^{-2}$	m/s
${\rho }_S$	2000.0	kg/m${}^3$
${\rho }_F$	1000.0	kg/m${}^3$
$K_F$	5.2 × 10${}^9$	Pa

. The external vertical impulsive load is the same as in the 2D benchmark case (Figure 4b), with equal duration; a time step size $\Delta t={10}^{-3}s$ is used in order to respect the CFL condition.

Three different soil compressibilities: $K_s=\infty$, $K_s=5.2\times {10}^9\mathrm{Pa}$ and $K_s=5.2\times {10}^7\mathrm{Pa}$, corresponding to three different Biot’s coefficients (Equation (13)): $B=1.0$, $B=0.998$ and $B=0.846$, have been taken into account. No appreciable variations in terms of soil displacement and velocities at node A (Figure 8a,b) are visible, with slight differences in terms of fluid velocity (Figure 8c). Far from the impulsive load, the smallest solid compressibility causes wider movements especially along the propagation direction of the impulsive load. When considering an isotropic material, the Rayleigh waves along X and Y presents the same shape and magnitude (see Figure 8d). In Figure 8e, the pore-fluid pressure evolution at nodes E and G is plotted: in case of incompressibility, the peak pressure occurs simultaneously with the external load and then dissipates during the analysis. In the case of larger compressibility, a delay in the peak appears.

By considering the highest compressibility ($K_s=5.2\times {10}^7\mathrm{Pa}$), the deformation states of Figure 9 show once again the triggering and evolution of shear and Rayleigh waves, already visible through 2D analysis.

4.3. Dynamic Poroelastic Responses with Transversely Isotropic Porous Media

The behavior of an anisotropic and fully saturated porous material replicated by an elastic transversely isotropic constitutive model for both solid material (intrinsic) and porous matrix (structural) is here considered (see Figure 10a), with material parameters shown in

Table 3. Transversely isotropic soil: material parameters for the 3D model.

Parameter	Values	S.I. Unit
$E_x,E_y$	9 × 10${}^6$	Pa
$E_z$	15 × 10${}^6$	Pa
${\nu }_{xy},{\nu }_{yx}$	0.25
${\nu }_{yz},{\nu }_{xz}$	0.21
${\nu }_{zx},{\nu }_{zy}$	0.35
$G_{xy}=\frac{E_x}{2(1+{\nu }_{xy})}$	3.6 × 10${}^6$	Pa
$G_{23},G_{31}$	6.0 × 10${}^6$	Pa
$n_0^F$	0.33
$k_x^F,k_y^F$	10${}^{-2}$	m/s
$k_z^F$	10${}^{-4}$	m/s
${\rho }_S$	2000	kg/m${}^3$
${\rho }_F$	1000	kg/m${}^3$
$K_m$	7.14 × 10${}^6$	Pa
$K_S$	3.57 × 10${}^9$	Pa
$K_F$	2.2 × 10${}^9$	Pa

.

4.3.1. Effect of Different Rotation in a Transversely Isotropic Symmetry Axis of Soil Material

Due to structural anisotropy, half of the full 3D domain (dashed blue line) has been taken into account. We assume to rotate by $\alpha$ the y-axis of the isotropy plane (solid phase and porous material) with respect to the y-axis of the global reference system (see Figure 10). A series of analyses have been performed changing from $\alpha =0{}^{\circ }$ to 90${}^{\circ }$, as shown in Figure 10b. From Equation (11), the rotation of the isotropic plane activates the coupling components of the elastic (or compliance) tensor as schematically represented in Figure 11, and then it will conduct to different soil responses. The coaxiality assumption between the structural and intrinsic tensor leads to obtaining, as usual, a constant Biot’s coefficient tensor $\mathit{a}=B\mathit{I}=0.998\mathit{I}$, so the fluid pressure interacts only with the volumetric stresses. The permeability tensor is also anisotropic and coaxial with the direction of material anisotropy; the permeability constants are shown in Table 3, the lowest permeability assumed along the z-direction.

A summary of the main results considering three different rotations is shown in Figure 12. By increasing $\alpha$, the soil stiffness along with z decreases, and it increases along with y. This leads to an increase in peak values of the vertical displacement at node A (Figure 12a) and in the solid velocity field (Figure 12b), together with a decrease in the peak values of fluid velocity at node A (Figure 12c). Evidently, by considering no rotation of the isotropic plane, $\alpha =0$, the motion of nodes B and D is the same (black continuous line; Figure 12d,e), while increasing $\alpha$ leads to different motions. This allows for obtaining two different Rayleigh waves along x and y, being Rayleigh waves linear combinations of P- and S- waves at the surface. In Figure 12f, the surface motion at node C is plotted: it can be observed that, by increasing $\alpha$, a horizontal motion orthogonal to the wave propagation appears, so along direction x=y, representing a Love wave (horizontal shifting) that consequently is coupled with the Rayleigh one.

In Figure 13, the pore pressure evolution is reported for two different pairs of nodes (close and far from the impulsive source respectively, see Figure 7) belonging to orthogonal directions: a comparison with the isotropic case leads to observing that the behavior is now strongly different (and no longer superimposed, see Figure 8e), with aligned peak values evidencing a high speed pressure velocity.

By analyzing soil deformation (Figure 14) where $\alpha =90{}^{\circ }$, the different shear waves propagate with different speed and magnitude within the soil medium; hence, the shear wave on the plane YZ reaches the boundary before the one on the XZ plane.

An alternative way to appreciate the shear wave splitting is plotting the effective shear stresses evolution along two planar orthogonal directions $a=\{\forall (x,y,z)\in {\mathcal{R}}^3|x\in \mathcal{R},y=10.5m,z=9.5m\}$ and $b=\{\forall (x,y,z)\in {\mathcal{R}}^3|y\in \mathcal{R},x=10.5m,z=9.5m\}$, Figure 15, considering two different rotations of the material axis. For $\alpha =0.0{}^{\circ }$ (continuous solid line), the propagation of effective shear stresses along a and b is the same, see, in fact, solid lines of Figure 15a,b and those of Figure 15c,d; whereas, considering $\alpha =45{}^{\circ }$ (marked dashed line), a different behavior is visible. By varying the angle, the shear stress ${\tau }_{xy}$ speed decreases along a but increases along b. In the case of mutual stresses ${\tau }_{xz}$ and ${\tau }_{yz}$, both speeds increase in different ways.

Particularly, Figure 16 summarizes the reported behavior, i.e., normalized peak stresses and their relative velocity show and confirm a measure of shear wave splitting.

4.3.2. Effect of Biot’s Effective Stress Coefficient Tensor on Wave Propagation

By adopting the same geometric model as before, a series of analyses have been additionally performed by fixing the material parameters of the porous material and rotating the plane of isotropy of the solid phase (intrinsic anisotropy), from $\alpha =0{}^{\circ }$ to 90${}^{\circ }$ as shown in Figure 10c. This assumption leads to obtaining a full Biot’s coefficient tensor:

$$\mathit{A}\left(\alpha \right)=\mathit{I}-\mathit{C}·\left[{\mathit{M}}^{-T}\left(\alpha \right){\mathit{S}}_S\left(0{}^{\circ }\right){\mathit{M}}^{-1}\left(\alpha \right)\right]=\left[\begin{array}{ccc}{a}_{11}& 0& 0\\ 0& a_{22}& a_{23}\\ 0& a_{23}& a_{33}\end{array}\right]$$

(44)

with $\mathit{M}\left(\alpha \right)$ the fourth-order transformation tensor and ${\mathit{S}}_S\left(0{}^{\circ }\right)$ the fourth-order compliance constitutive tensor of the solid phase only. In this case, two non zero extra-terms are obtained and, together with the values on the diagonal, the fluid part interacts with the deviatoric stress of the soil. For the initial case of $\alpha =0{}^{\circ }$, the material parameters of the intrinsic transversely isotropic model are listed in

Table 4. Intrinsic transversely isotropic model: material parameters.

Parameter	Values	S.I. Unit
$E_x,E_y$	1.8 × 10${}^7$	Pa
$E_z$	3.0 × 10${}^7$	Pa
${\nu }_{xy},{\nu }_{yx}$	0.25
${\nu }_{yz},{\nu }_{xz}$	0.21
${\nu }_{zx},{\nu }_{zy}$	0.35
$G_{xy}$	7.2 × 10${}^6$	Pa
$G_{yz},G_{zx}$	1.2 × 10${}^7$	Pa

. The anisotropic permeability remains the same as in the previous analyses, and it is coaxial with the symmetry axis of the structural anisotropy.

By varying the Biot’s tensorial coefficients through $\alpha$, no differences are appreciable in terms of vertical displacements (Figure 17a) and vertical fluid velocity (Figure 17b) at node A, while two different plane motions at nodes B and D are clearly visible. These latter graphs allow us to catch the Rayleigh waves spreading horizontally and to estimate the splitting of the shear waves. In case of high compressibility, different fluid pressure velocities are shown (see Figure 17e,f), the values being no longer aligned as previously evidenced.

By considering a rotation of 90${}^{\circ }$ (this one only for sake of brevity), Figure 18, the splitting of the shear wave is again visible, although less evident than in the previous section. The deformed mesh and the Euclidean norm of the displacements in different time steps of the analysis are plotted considering the solid phase rotated by 90${}^{\circ }$ with respect to the symmetry axis of transversely isotropic model of the mixture (structural anisotropy). A higher solid stiffness phase along Y leads to a decrease in the effective stress along the same direction, provoking a slower shear waves along Y and faster along X.

The results in terms of effective shear stresses (Figure 19a) confirm the behavior already appreciable when rotating the structural anisotropy; when choosing e.g., the second invariant of the Biot’s coefficient tensor as reference (Figure 19b), it can be noticed that higher angles correspond to lower values, with an anomalous splitting mechanism.

5. Discussion

The results indicate that the numerical model is able to replicate the following physical phenomena, i.e.,:

(i)

the p-waves produce polarized vibrations along the direction of propagation (particles move along the wave’s direction of propagation) and subsequent compression and extension deformations along the same direction: they are visible along the vertical direction under the impulsive load (Figure 6 and Figure 9), even considering the anisotropic models (in this case they are coupled with the shear contribution, Figure 14 and Figure 18);

(ii)

p-waves are faster than s-waves: in all the models in fact the domain borders are reached in different times;

(iii)

s-waves generate polarized vibrations on a plane containing the direction of propagation and shear deformation (Figure 6, Figure 9, Figure 14 and Figure 18);

(iv)

the s-wave decouples into a wave polarized on the horizontal plane and into another one on the vertical plane: visible in the curves of effective shear stresses, Figure 15;

Particularly, for s-waves, it has been evidenced that, when propagating vertically, they are always polarized on a vertical plane, whereas, in case of horizontal propagation, one component of the wave belongs to an horizontal plane, the other on a vertical plane. This particularly happens when isotropy or transversely isotropy with vertical isotropy axis is assumed; otherwise, they show inclined polarizations.

With regard to the surface waves, the Rayleigh waves have been properly reproduced: such waves propagate according to cylindrical wavefronts; the resulting motion on the vertical plane is retrograde elliptical (Figure 12d,e and Figure 17c,d) (compare, e.g., with Yang and Li [44]). If considering Love waves (Figure 12f, the generated horizontal vibrations clearly appear polarized along a direction orthogonal to the propagation one (shear deformations).

Even the geometric attenuation of the waves seems to be properly reconstructed: their energetic content being reduced at a far distance from the source, the amplitude of the medium displacement correspondingly decreases (geometric damping): this is evidenced in the effective shear stress curves themselves, Figure 15 and Figure 19a,b. More importantly, when numerically reproducing an anisotropic soil domain, the physical phenomenon of waves splitting appears to be reproduced as well (Figure 16 and Figure 19c,d) with generation of waves with different intensity and speed.

6. Conclusions

The propagation of waves in soils, developing from a point source of a dynamic load, have been analyzed with attention focused on polarization and shear wave splitting due to anisotropy of the permeability tensor, anisotropy of the solid skeleton, as well as to a novel Biot’s tensor. The mathematical-numerical model adopts a u–v–p formulation enhanced by the introduction of Taylor–Hood mixed finite elements and comparisons with different integration strategies have revealed to prefer a semi-explicit/implicit scheme with equal order interpolation due to its satisfied stability requirements. The numerical implementation of both an anisotropic permeability tensor together with a Biot’s tensor has allowed for boosting the contribution of anisotropy when reproducing waves polarization and splitting: the conducted analyses have correctly reproduced, to recall a few details, polarized vibrations along the direction of propagation produced by P-waves, polarized vibrations on a plane containing the direction of propagation and shear deformation generated by S-waves, the generation of surface waves, and, more importantly, waves splitting due to intrinsic or structural anisotropy, with enhanced effects when considering the coupling between the volumetric fluid pressure and shear stresses.