Impacts of saturation-dependent anisotropy on the shrinkage behavior of clay rocks

Abstract

Geomaterials such as soils and rocks can exhibit inherent anisotropy due to the preferred orientation of mineral grains and/or cracks. They can also be partially saturated with multiple types of fluids occupying the pore space. The anisotropic and unsaturated behaviors of geomaterials can be highly interdependent. Experimental studies have shown that the elastic parameters of rocks evolve with saturation. The effect of saturation has also been shown to differ between directions in transversely isotropic clay rock. This gives rise to saturation-dependent stiffness anisotropy. Similarly, permeability anisotropy can also be saturation-dependent. In this study, constitutive equations accommodating saturation-dependent stiffness and hydraulic anisotropy are presented. A linear function is used to describe the relationship between the elastic parameters and saturation, while the relative permeability–saturation relationship is characterized with a log-linear function. These equations are implemented into a hydromechanical framework to investigate the effects of saturation-dependent properties on the shrinkage behavior of clay rocks. Numerical simulations are presented to demonstrate the role of saturation-dependent stiffness and hydraulic anisotropy in shrinkage behavior. The results highlight that strain anisotropy and time evolution of pore pressures are substantially influenced by saturation-dependent stiffness and hydraulic anisotropy.

Appendices

Appendix 1: Linearized system

Equations (31) and (32) are solved simultaneously by first writing them in residual form:

$$\begin{aligned}\varvec{{\mathcal {R}}}_1 &=\int _{{\mathcal {B}}}\varvec{B}^{{\mathsf {T}}}\{\varvec{\sigma }'-\psi ^wp^h\varvec{b}\} \,\mathrm{d}V- \int _{{\mathcal {B}}}\varvec{N}^{{\mathsf {T}}}\rho \varvec{g}\,\mathrm{d}V \nonumber \\&\qquad - \int _{\partial {{\mathcal {B}}}_t} \varvec{N}^{{\mathsf {T}}} \mathring{\varvec{t}} \,\mathrm{d}A \end{aligned}$$

(35)

and

$$\begin{aligned}&\varvec{{\mathcal {R}}}_2 \nonumber \\&\quad = -\int _{{\mathcal {B}}}\hat{\varvec{N}}{}^{{\mathsf {T}}}\psi ^{w}\varvec{b}:(\nabla ^s\varvec{u}^h - \nabla ^s\varvec{u}^h_n) \,\mathrm{d}V - \int _{{\mathcal {B}}}\hat{\varvec{N}}{}^{{\mathsf {T}}}c(\psi ^w-\psi ^w_n) \,\mathrm{d}V\nonumber \\&\qquad -\int _{{\mathcal {B}}}\hat{\varvec{N}}{}^{{\mathsf {T}}}\frac{1}{{\mathcal {M}}}(p^h-p^h_n) \,\mathrm{d}V+\varDelta t \int _{{\mathcal {B}}} \varvec{E}^{{\mathsf {T}}}\varvec{q} \,\mathrm{d}V \nonumber \\&\qquad + {\varDelta } t \int _{\partial {{\mathcal {B}}}_q} \hat{\varvec{N}}{}^{{\mathsf {T}}} \mathring{q}\,\mathrm{d}A .\qquad \end{aligned}$$

(36)

From Remark 2, we drop the $\nabla p^h\cdot \varvec{q}/K_w$-term and simply set $\varvec{E} = \nabla \hat{\varvec{N}}$. Furthermore, since the variation of Biot tensor $\varvec{b}$ is very small, we shall assume it to be constant in what follows.

The residuals defined above are highly nonlinear with respect to the unknown vectors $\varvec{x}_1=\varvec{d}$ and $\varvec{x}_2=\varvec{p}$, so a Newton–Raphson iteration scheme is chosen to solve this problem. The linearized system is defined by

$$\begin{aligned} \begin{bmatrix}\varvec{K}_{11}&{}\varvec{K}_{12}\\ \varvec{K}_{21}&{}\varvec{K}_{22}\\ \end{bmatrix} \begin{Bmatrix}\delta \varvec{x}_1\\ \delta \varvec{x}_2\end{Bmatrix} = \begin{Bmatrix}\varvec{{\mathcal {R}}}_1\\ \varvec{{\mathcal {R}}}_2\end{Bmatrix} , \end{aligned}$$

(37)

where $\delta \varvec{x}_1=\delta \varvec{d}$ and $\delta \varvec{x}_2=\delta \varvec{p}$ are the search directions [13]. The submatrices of the tangent operator take the following forms:

$$\begin{aligned} \varvec{K}_{11}=\frac{\partial \varvec{{\mathcal {R}}}_1}{\partial \varvec{x}_1} = \int _{{\mathcal {B}}} \varvec{B}^{{\mathsf {T}}}\varvec{C}^e\varvec{B} \,\mathrm{d}V , \end{aligned}$$

(38)

where $\varvec{C}^e$ is the elasticity tensor in matrix form.

The coupling operators are

$$\begin{aligned}\varvec{K}_{12}&=\frac{\partial \varvec{{\mathcal {R}}}_1}{\partial \varvec{x}_2} =-\int _{{\mathcal {B}}}\varvec{B}^{{\mathsf {T}}}\varvec{\xi }\hat{\varvec{N}}\,\mathrm{d}V \nonumber \\&\qquad - \int _{{\mathcal {B}}} \varPsi \varvec{B}^{{\mathsf {T}}}\{\varvec{b}\}\hat{\varvec{N}}\,\mathrm{d}V + \int _{{\mathcal {B}}}r \varvec{N}^{{\mathsf {T}}}\varvec{g}\hat{\varvec{N}} \,\mathrm{d}V , \end{aligned}$$

(39)

where

$$\begin{aligned} \varvec{\xi }=\frac{\partial \{\varvec{\sigma }'\}}{\partial \psi ^w}\frac{\partial \psi ^w}{\partial s},\qquad \varPsi = \psi ^w-\frac{\partial \psi ^w}{\partial s}p^h,\qquad r=\frac{\partial \rho }{\partial \psi ^w}\frac{\partial \psi ^w}{\partial s}, \end{aligned}$$

(40)

and

$$\begin{aligned} \varvec{K}_{21}=\frac{\partial \varvec{{\mathcal {R}}}_2}{\partial \varvec{x}_1} = -\int _{{\mathcal {B}}}\psi ^w\hat{\varvec{N}}{}^{{\mathsf {T}}}\{\varvec{b}\}^{{\mathsf {T}}}\varvec{B} \,\mathrm{d}V . \end{aligned}$$

(41)

For fully saturated condition, $\varvec{\xi }=\varvec{0}$, $\varPsi =\psi ^w=1$, and $r=0$, which gives $\varvec{K}_{12}=\varvec{K}_{21}$. In the unsaturated range, the system is asymmetric.

The structure of (2,2) submatrix block is given by

$$\begin{aligned} \varvec{K}_{22}&=\frac{\partial \varvec{{\mathcal {R}}}_2}{\partial \varvec{x}_2} =-\varDelta t\int _{{\mathcal {B}}}\varvec{E}^{{\mathsf {T}}}\bar{\varvec{K}}\varvec{E}\,\mathrm{d}V \nonumber \\&\qquad +\varDelta t\int _{{\mathcal {B}}}\varvec{E}^{{\mathsf {T}}}\varvec{\varPhi }\hat{\varvec{N}}\,\mathrm{d}V +\int _{{\mathcal {B}}}\hat{\varvec{N}}{}^{{\mathsf {T}}}\varTheta \hat{\varvec{N}}\,\mathrm{d}V, \end{aligned}$$

(42)

where $\bar{\varvec{K}}=\varvec{k}_\mathrm{rel} \cdot \varvec{\kappa }/\mu _w$ is the matrix of hydraulic conductivity,

$$\begin{aligned} \varvec{\varPhi }=\frac{\partial \bar{\varvec{K}}}{\partial \psi ^w}\frac{\partial \psi ^w}{\partial s}\big (\nabla p^h - \rho _w\varvec{g}\big ), \end{aligned}$$

(43)

and

$$\begin{aligned} \varTheta&= {} \varvec{b}:\nabla ^s\big (\varvec{u}^h-\varvec{u}^h_n\big )\frac{\partial \psi ^w}{\partial s} + c\frac{\partial \psi ^w}{\partial s}+\big (\psi ^w-\psi ^w_n\big )\frac{\partial c}{\partial s}\nonumber \\&\quad - {} \frac{1}{{\mathcal {M}}}+\big (p^h-p^h_n\big )\frac{\partial (1/M)}{\partial \psi ^w}\frac{\partial \psi ^w}{\partial s} . \end{aligned}$$

(44)

For fully saturated condition, $\varvec{\varPhi }=\varvec{0}$ and $\varTheta =-1/{{\mathcal {M}}}$, and

$$\begin{aligned} \frac{1}{{\mathcal {M}}} =\frac{\phi }{K_w}+\frac{\beta }{K_s} \end{aligned}$$

(45)

is the Biot modulus for fully saturated media.

Appendix 2: Stability of transversely isotropic elastic material

Ting and Chen [89] showed that for a transversely isotropic material with positive elastic moduli to be stable, the following inequality must be satisfied:

$$\begin{aligned} s_{55}\big [2\varDelta _2-s_{11}(s_{22}-s_{23})\big ]>0 , \end{aligned}$$

(46)

where $s_{11}$, $s_{22}$, $s_{23}$ and $s_{55}$ are components of the compliance tensor in Voigt notation and $\varDelta _2=s_{11}s_{22}-s_{12}^2$. This inequality can be expressed with the elastic constants

$$\begin{aligned} \frac{1}{G_{12}}\Big [2\varDelta _2-\frac{1}{E_1}\Big (\frac{1}{E_2}+\frac{\nu _{23}}{E_2}\Big )\Big ]>0 , \end{aligned}$$

(47)

where

$$\begin{aligned} \varDelta _2=\frac{1}{E_1 E_2}-\Big (\frac{\nu _{12}}{E_1}\Big )^2 . \end{aligned}$$

(48)

If all elastic moduli are positive, the equations can be simplified to

$$\begin{aligned} 2(E_1-E_2{\nu ^2_{12}})-E_1(1+\nu _{23})>0 , \end{aligned}$$

(49)

which can be rearranged to give

$$\begin{aligned} 1-\nu _{23}>2 {\nu ^2_{12}} \frac{E_2}{E_1} . \end{aligned}$$

(50)