Fully-automated adaptive mesh refinement for media embedding complex heterogeneities: application to poroelastic fluid pressure diffusion

Abstract

Relating the attenuation and velocity dispersion of seismic waves to fluid pressure diffusion (FDP) by means of numerical simulations is essential for constraining the mechanical and hydraulic properties of heterogeneous porous rocks. This, in turn, is of significant importance for a wide range of prominent applications throughout the Earth, environmental, and engineering sciences, such as, for example, geothermal energy production, hydrocarbon exploration, nuclear waste disposal, and CO₂ storage. In order to assess the effects of wave-induced FDP in heterogeneous porous rocks, we simulate time-harmonic oscillatory tests based on a finite element (FE) discretization of Biot’s equations in the time-frequency domain for representative elementary volumes (REVs) of the considered rock masses. The major challenge for these types of simulations is the creation of adequate computational meshes, which resolve the numerous and complex interfaces between the heterogeneities and the embedding background. To this end, we have developed a novel method based on adaptive mesh refinement (AMR), which allows for the fully automatic creation of meshes for strongly heterogenous media. The key concept of the proposed method is to start from an initially uniform coarse mesh and then to gradually refine elements which have non-empty overlaps with the embedded heterogeneities. This results in a hierarchy of non-uniform meshes with a large number of elements close to the interfaces, which do, however, not need to be explicitly resolved. This dramatically simplifies and accelerates the laborious and time-consuming process of meshing strongly heterogeneous poroelastic media, thus enabling the efficient simulation of REVs containing heterogeneities of quasi-arbitrary complexity. After a detailed description of the methodological foundations, we proceed to demonstrate that the FE discretization with low-order FE has a unique solution and hence does not present spurious modes. We assess the practical effectiveness and accuracy of the proposed method by means of four case studies of increasing complexity.

Appendices

Appendix A: weak formulation of Biot’s equations

1.1 A.1 Mathematical preliminaries

We denote by $\mathbb C$ the set of complex numbers and by $\mathrm {j}~=~\sqrt {-1}$ the imaginary unit. For a number $a~\in ~\mathbb C$, the symbol a^∗ will denote the complex conjugate, $\Re {\{a\}}\in \mathbb R$ the real part, and $\Im {\{a\}}\in \mathbb R$ the imaginary part. Given two elements $a,b\in \mathbb C$, we introduce the product $(a,b)_{\mathbb C} : =a b^{*}$ so that the modulus of a can be computed as $|a|~=~\sqrt {(a,a)_{\mathbb C}}$. $\mathbb V$ is the complex Euclidean space $\mathbb C^{d}$ provided with the inner product

$$ (\underline{u} , \underline{v})_{\mathbb V} := \underline{u} \cdot \underline{v}^{*} = \sum\limits_{i=1}^{d} u_{i} {v}_{i}^{*}, \quad \underline{u}, \underline v\in\mathbb V. $$

(15)

The space of the second-order tensors over $\mathbb V$ is denoted by $\mathbb M$. Once a basis is introduced over $\mathbb V$, the space $\mathbb M$ can be identified with $\mathbb C^{d\times d}$ and we provide it with the scalar product

$$ (\underline{\underline{\sigma}}, \underline{\underline{\tau}})_{\mathbb M} := \underline{\underline{\sigma}} : \underline{\underline{\tau}}^{*} = \sum\limits_{i=1}^{d} \sum\limits_{k=1}^{d} \sigma_{ik} \tau_{ik}^{*}, \quad \underline{\underline{\sigma}},\underline{\underline{\tau}}\in\mathbb M. $$

(16)

For a tensor σ̲̲ and the symbol σ̲ ̲^H denotes the adjoint whose components are $\sigma _{ik}^{H} = (\sigma _{ki})^{*}$ and the symbol σ^T denotes the transpose whose components are $\sigma _{ik}^{T}~=~\sigma _{ki}$. The symmetric part of a tensor is denoted by

$$ \text{sym} \underline{\underline{\sigma}} = \frac{ \underline{\underline{\sigma}} + \underline{\underline{\sigma}}^{T}}{2}, $$

(17)

while the Hermitian part is denoted by

$$ \text{Herm} \underline{\underline{\sigma}} = \frac{ \underline{\underline{\sigma}}+ \underline{\underline{\sigma}}^{H}}{2}. $$

(18)

The symbol $\text {tr}(\underline {\underline {\sigma }})$ is the usual trace operator. The operators ^∗, R, and, I are intended to be applied component-wise to a vector or a tensor.

Before deriving the weak formulation of Biot’s equations, we first introduce the suitable function spaces. Given a vector space $\mathbb W$, which can be $\mathbb R$, $\mathbb C$, $\mathbb V$, or $\mathbb M$, for $p \in [1,\infty )$, we denote with $L^{p}({\Omega },\mathbb W)$ the space of p-integrable functions defined on Ω with values in $\mathbb W$. The space $L^{2}({\Omega }, \mathbb W)$ is endowed with the scalar product

$$ (p,q)_{L^{2}({\Omega}, \mathbb W)} := {\int}_{\Omega} (p,q)_{\mathbb W} \mathrm{d} {\Omega}, \quad p,q\in L^{2}({\Omega}, \mathbb W), $$

(19)

and the norm

$$ \| p \|_{L^{2}({\Omega}, \mathbb W)} = \sqrt{(p,p)_{L^{2}({\Omega}, \mathbb W)}} . $$

(20)

Instead, $L^{\infty }({\Omega },W)$ is the space of the essentially bounded functions defined on Ω. The operators ∇ and div are the usual gradient and divergence operators, respectively. They are intended to be applied to the real and the imaginary parts separately. For example, given a function $p : {\Omega } \to \mathbb C$, the gradient is defined as

$$ \nabla p = \nabla \Re{\{p\}} + \j \nabla \Im{\{p\}}. $$

(21)

The symmetric gradient of a function $\underline {v} : {\Omega } \to \mathbb V$ is denoted by

$$ \underline{\underline{\varepsilon}} (\underline{v}) = \text{sym} \nabla \underline{v} $$

(22)

and the following identity holds

$$ \text{tr} (\underline{\underline{\varepsilon }}(\underline{v}) ) = \text{tr} (\nabla \underline{v}) = \underline{\underline{I}} : \nabla \underline{v} = \text{div} \underline{v}, $$

(23)

where $\underline {\underline {I}}$ is the second-order identity tensor. We define the following spaces:

$$ H^{1}({\Omega}, \mathbb C) = \{ p\in L^{2}({\Omega} , \mathbb C) | \nabla p\in L^{2}({\Omega} , \mathbb V) \} $$

(24)

and

$$ H^{1}({\Omega}, \mathbb V) = \{ \underline{v}\in L^{2}({\Omega} , \mathbb V) | \nabla \underline{v}\in L^{2}({\Omega} , \mathbb M) \} $$

(25)

endowed with the standard H¹ scalar products. We denote as

$$ H^1_{\Gamma}({\Omega}, \mathbb W) = \{ w \in H^1({\Omega}, \mathbb W) | {w} |_{{\Gamma}^i_{L}} = {w}|_{{\Gamma}^i_{-L}},\quad \text{with } i = 1,\ldots,d \}. $$

(26)

Finally, in order to simplify the notation, we define the following function spaces $V=H^1_{\Gamma }({\Omega }, \mathbb V)$ and $Q=H^1_{\Gamma }({\Omega }, \mathbb C)$.

1.2 A.2 Weak formulation of Biot’s equations with periodic boundary conditions

In order to formally derive the weak formulation of the system of Eqs. 1 and 2, we first consider a test function $\underline {v}\in V$, multiply (1) by $\underline {v}^{*}$, and integrate over Ω

$$ {\int}_{\Omega} (\text{div } \underline{\underline{\sigma}}) \cdot \underline{v}^{*} \mathrm{d}{\Omega}=0. $$

(27)

Using Green’s formula, we obtain

$$ {\int}_{\Omega} \underline{\underline{\sigma}} : \nabla \underline{v}^{*} \mathrm{d}{\Omega}= \sum\limits_{{\mathrm s}\in \{ -1,1 \}} \sum\limits_{i=1}^{d} {\int}_{{\Gamma}^{i}_{{\mathrm s} \mathrm L}} \underline{\underline{\sigma}} \underline{n} \cdot \underline{v}^{*} \mathrm{d}{\Gamma}. $$

(28)

The symmetry of $\underline {\underline {\sigma }}$ implies that

$$ \underline{\underline{\sigma}} : \nabla \underline{v}^{*} = \underline{\underline{\sigma}} : \underline{\underline{\varepsilon}}({ \underline{v}^{*}}) $$

(29)

and, finally, using the fact that the test functions belong to V and the boundary conditions (7), we observe that the boundary integral at the right-hand side vanishes giving

$$ {\int}_{\Omega} \underline{\underline{\sigma}} : {\underline{\underline{\varepsilon}}}({ \underline{v}^{*}}) \mathrm{d}{\Omega}= 0. $$

(30)

Similarly, for Eq. 2, we consider a test function q ∈ Q. Multiplying this equation by q^∗ and integrating over Ω, we obtain

$$ - \mathrm{j} {\int}_{\Omega} \alpha \text{div} \underline{u} q^{*} \mathrm{d}{\Omega} - \mathrm{j} {\int}_{\Omega}\frac{p}{M} q^{*} \mathrm{d}{\Omega} + \frac{1}{\omega} {\int}_{\Omega} \text{div} \left( \frac{k}{\eta}\nabla p\right) q^{*} \mathrm{d}{\Omega}=0. $$

(31)

Applying again Green’s formula, the last term becomes

$$ - {\int}_{\Omega}\left( \frac{k}{\eta}\nabla p\right) \cdot \nabla q^{*} \mathrm{d}{\Omega} + \sum\limits_{{\mathrm s}\in\{ -1,1 \}} \sum\limits_{i=1}^{d} {\int}_{{\Gamma}^{i}_{{\mathrm s} \mathrm L}} \frac{k}{\eta}\nabla p \cdot \underline{n} q^{*}|_{{\Gamma}^{i}_{{\mathrm s} \mathrm L}} \mathrm{d}{\Gamma}. $$

(32)

As for the linear momentum equation, the boundary integral vanishes and we obtain

$$ - \mathrm{j} {\int}_{\Omega} \alpha \text{div} \underline{u} q^{*} \mathrm{d}{\Omega} - \mathrm{j} {\int}_{\Omega} \frac{1}{M} p q^{*} \mathrm{d}{\Omega} - \frac{1}{\omega} {\int}_{\Omega} \frac{k}{\eta}\nabla p \cdot \nabla q^{*} \mathrm{d}{\Omega}=0. $$

(33)

In order to include the external relative displacement, we define the set

$$ U = \{ \underline{v}\in H^{1}({\Omega}, \mathbb V) | \underline{v} |_{{\Gamma}^{i}_{\mathrm L}}- \underline{v}|_{{\Gamma}^{i}_{-\mathrm L}} = \underline{\alpha}_{i} \quad \text{with } i=1,\ldots,d \}. $$

(34)

Hence, exploiting the definition of stress (4), the weak formulation of Biot’s equation can be written as Find $(\underline {u}, p) \in U\times Q$ such that

$$ \begin{array}{r c c c c c c} & a(\underline{u} , \underline{v}) & - & b(\underline{v}^{*} , p^{*} ) & = & 0 & \quad \forall \underline{v}\in V, \\ - & \mathrm{j} b(\underline{u} , q) & - & d(p,q ) & = & 0 & \quad \forall q\in Q, \end{array} $$

(35)

where

$$ \begin{array}{r c l} a(\underline{u} , \underline{v}) &= & {\int}_{\Omega} 2\mu \underline{\underline{\varepsilon}}(\underline{u}) : \varepsilon(\underline{v}^{*}) + \lambda \text{div} \underline{u} \text{div} \underline{v}^{*} \mathrm{d}{\Omega},\\ b(\underline{u} , q) &= & {\int}_{\Omega} \alpha \text{div} \underline{u} q^{*} \mathrm{d}{\Omega}, \\ d(p , q) &= & \mathrm{j} m(p,q) + c(p,q),\\ m(p,q) &=& {\int}_{\Omega} \frac{1}{M} p q^{*} \mathrm{d}{\Omega}, \\ c(p,q) &=& {\int}_{\Omega} \frac{k}{ \eta}\nabla p \cdot \nabla q^{*} \mathrm{d}{\Omega}. \end{array} $$

(36)

In Eq. 35, the right-hand side is null since we assume no external force is acting on the REV. The external loads are encoded in the essential boundary conditions in the set U. In order to formally express these boundary conditions, the displacement can be written as

$$ \underline{u}= \underline{u}_{0} + \underline{u}_{\underline{\boldsymbol{\alpha}}}, $$

(37)

where $\underline {u}_{0}$ belongs to V and $\underline {u}_{\underline {\boldsymbol {\alpha }}}$ is a lifting function, which can be any element in U. Hence, for a given $\underline {u}_{\underline {\boldsymbol {\alpha }}}$, the problem (35) can be written as Find (u̲₀,p) ∈ V × Q such that

$$ \begin{array}{r c c c c c c} & a(\underline{u}_{0} , \underline{v}) & - & b(\underline{v}^{*} , p^{*} ) & = & -a(\underline{u}_{\underline{\boldsymbol{\alpha}}} , \underline{v}) & \quad \forall \underline{v}\in V, \\ - & \mathrm{j} b(\underline{u}_{0} , q) & - & d(p,q ) & = & \mathrm{j} b(\underline{u}_{\underline{\boldsymbol{\alpha}}} , q) & \quad \forall q\in Q. \end{array} $$

(38)

The solution can then be computed from Eq. 37. Finally, we also observe that $\underline {u}_{0}$ is defined up to a constant vector, since for any constant $\underline {c}$, $(\underline {u}_{0}+ \underline {c},p)$ is also a solution of Eq. 38. In order to eliminate this ambiguity, the space V can be replaced with

$$ V_{0}=\left\{ \underline{v} \in V : {\int}_{\Omega} \underline{v} \mathrm{d}{\Omega} = \underline{0} \right\}. $$

(39)

1.3 A.3 Properties of bilinear forms

The analysis of problem (38) is not trivial. Similar problems on real function spaces have been considered in [12]. Here, we limit ourselves to prove two properties that are useful for the analysis of the discrete problem.

While for Hermitian complex bilinear forms, the proof of coercivity is just a simple extension of the real case, it needs to be generalized for non-Hermitian bilinear forms [45]. In particular, we say that a bilinear form c(⋅,⋅) defined over a complex Hilbert space Q is $\mathbb C$-coercive, if there exists a constant γ such that $|c(p,p)|>\gamma \| p \|_{Q}^{2}$ for all p ∈ Q.

Property 1

The bilinear form a(⋅,⋅) is coercive over V₀ with coercivity constant

$$ \min(\mu_{b},\mu_{f}), $$

(40)

where μ is the shear modulus. This can be seen observing that $a(\underline{u}, \underline{u})= a(\Re \underline{u}, \Re \underline{u}) + a(\Im \underline{u}, \Im \underline{u})$. Applying Korn’s inequality to both bilinear forms at the right-hand side, coercivity can be easily proven [9].

Property 2

The bilinear form d(⋅,⋅) is $\mathbb C$-coercive over Q, with coercivity constant $\delta = \frac {1}{2} \min \limits (\frac {c_{0}}{\omega } , m_{0} )$, where

$$c_{0}=\min({k_{b}/\eta_{b}, k_{f}/\eta_{f} } ) \textrm{ and } m_{0}=\min(1/M_{b}, 1/M_{f}).$$

We observe that

$$ | d(u,u)| \geq \Re d(p,p) = \frac{1}{\omega} {\int}_{\Omega} \frac{k}{\eta} \nabla p \cdot \nabla p^{*} \mathrm{d}{\Omega} \geq \frac{1}{\omega} \min(k_{b}/\eta_{b},k_{f}/\eta_{f}) \| \nabla p\|^{2}_{L^{2}({\Omega}, \mathbb V)} $$

(41)

and

$$ | d(u,u)| \geq \Im d(u,u) = {\int}_{\Omega} \frac{1}{M} p^{2} \mathrm{d}{\Omega} \geq \min(1/M_{b}, 1/M_{f}) \| p\|^{2}_{L^{2}({\Omega}, \mathbb C)}. $$

(42)

Summing up (41) and (42), we obtain

$$ | d(p,p)| \geq \frac{1}{2} \min(m_{0}, \frac{c_{0}}{\omega} ) \| p \|^{2}_{V}. $$

(43)

B Finite element discretization of Biot’s equations

Let $\mathcal {T}$ be a 1-irregular triangulation of Ω. We denote by $\mathbb Q_{1}$ the space of real bilinear functions for d = 2 and of real trilinear functions for d = 3. We define the conforming interpolation spaces over $\mathcal {T}$

$$ X_{\mathcal{T}}= \{p_{h} \in C^{0}({\Omega},\mathbb R) : p_{h}|_{K}\in\mathbb Q_{1} \forall K \in \mathcal{T} \} $$

(44)

and

$$ Y_{\mathcal{T}}= \{ \underline{u}_{h}\in C^{0}({\Omega},\mathbb R^{d}) : (\underline{u}_{h})_{i} |_{K}\in\mathbb Q_{1} \forall K\in\mathcal{T} \text{ and } i=1,\ldots, d \}. $$

(45)

The dimension of the space $X_{\mathcal {T}}$ is N^r and the dimension of the space $Y_{\mathcal {T}}$ is dN^r. In order to include boundary conditions in the interpolation spaces above, we define $V_{\mathcal {T}}= Y_{\mathcal {T}} \cap V$, and $Q_{\mathcal {T}}=X_{\mathcal {T}} \cap Q$. For any $\underline {u}_{\underline {\boldsymbol {\alpha }}}$, the FE approximation of problem (38) has the following form:

Find $(\underline {u}_{0h}, p_{h}) \in V_{\mathcal {T}} \times Q_{\mathcal {T}}$ such that

$$ \begin{array}{r c c c c c c} & a(\underline{u}_{0h} , \underline{v}_{h}) & - & b(\underline{v}^{*}_{h} , p^{*}_{h} ) & = & -a(\underline{u}_{\underline{\boldsymbol{\alpha}}} , \underline{v}_{h}) & \quad \forall \underline{v}_{h} \in V_{\mathcal{T}}, \\ - & \mathrm{j} b(\underline{u}_{0h} , q_{h}) & - & d(p_{h},q_{h} ) & = & \mathrm{j} b(\underline{u}_{\underline{\boldsymbol{\alpha}}} , q_{h}) & \quad \forall q_{h}\in Q_{\mathcal{T}}. \end{array} $$

(46)

2.1 B.1 Uniqueness of the solution of the discretized of Biot’s equations

In order to prove the uniqueness of solution of problem (46), we first rewrite it in an algebraic form by introducing the basis functions $\{{\Phi }_{j} \in V_{\mathcal {T}}\}$ and $\{\phi _{l} \in Q_{\mathcal {T}}\}$ of $V_{\mathcal {T}}$ and $Q_{\mathcal {T}}$, respectively. Expanding the two components of the solution with respect to the basis functions, we obtain a linear system with the following block structure:

$$ \left| \begin{array}{c c} \mathbf{A} & -\mathbf{B}^{T}\\ - \mathrm{j} \mathbf{B} & -\mathbf{D} \end{array} \right| \left| \begin{array}{c} \mathbf{u} \\ \mathbf{p} \end{array} \right| = \left| \begin{array}{c} \mathbf{f} \\ \mathbf{g} \end{array} \right|, $$

(47)

where $\mathbf {A}\in \mathbb R^{d{N_{h}^{r}} \times d{N_{h}^{r}}}$, $\mathbf {B}\in \mathbb R^{{N_{h}^{r}} \times d{N_{h}^{r}}}$, and $\mathbf {D}\in \mathbb C^{{N_{h}^{r}} \times {N_{h}^{r}}}$ are the matrices related to the bilinear forms a, b, and d, respectively. The elements of such matrices are given by A_ij = a(Φ_j,Φ_i), B_kj = b(Φ_j,ϕ_k), and D_kl = d(ϕ_l,ϕ_k). The vectors $\textbf {u}\in \mathbb C^{d{N_{h}^{r}}}$ and $\textbf {p}\in \mathbb C^{{N_{h}^{r}}}$ are the vectors collecting the unknown Lagrange coefficients of the discrete displacement and pressure, respectively. The two vectors at the right-hand side are defined as f_i = −a(u̲_α,Φ_i) and g_k = jb(u̲_α,ϕ_k). The matrix D can be written as $\textbf {D}=\frac {1}{\omega } \textbf {C} + \mathrm {j} \textbf {M}$.

The stiffness matrix in the linear system (47) is indefinite and we show that it is invertible thanks to from properties 1 and 2. First, owing to property 1, A is Hermitian and semi-positive definite, having in its kernel the displacements associated with rigid body motions. Hence, u can be formally computed from

$$ \mathbf{u}= \mathbf{A}^{-1}(\mathbf F + \mathbf{B}^{T} \mathbf{p}), $$

(48)

where u is defined up to a vector in the kernel of A.

Replacing this equation in the second line of Eq. 47, we obtain

$$ \mathbf{S} \mathbf{p}= - \mathrm{j} \mathbf{B} \mathbf{A}^{-1} \mathbf{f} - \mathbf{g}, $$

(49)

where S = D + jBA^− 1B^T. Before proving that S is invertible, we observe that BA^− 1B^T is symmetric and semi-positive definite and M, being a scaled mass matrix, is symmetric and positive definite.

Property 3

The matrix S is regular.

Let us suppose by contradiction that S is not regular. This means there exists a non-trivial vector r associated with a function $r_{h}~\varepsilon ~Q_{\mathcal T}$, such that Sr = 0. Hence, by definition of S, we obtain

$$ \begin{aligned} 0=| \mathbf{r}^{T} \mathbf{S} \mathbf{r} | = | \frac{1}{\omega} \mathbf{r}^{T} \mathbf{C} \mathbf{r} + j \mathbf{r}^{T} (\mathbf{M} + \mathbf{B} \mathbf{A}^{-1} \mathbf{B}^{T} )\textbf{r} | > \\ \frac{1}{2} | \frac{1}{\omega} \mathbf{r}^{T} \mathbf{C} \mathbf{r}| + \frac{1}{2} | \mathbf{r}^{T} (\mathbf{M} + \mathbf{B} \mathbf{A}^{-1} \mathbf{B}^{T} ) \mathbf{r} | >\\ \frac{1}{2} | \frac{1}{\omega} \mathbf{r}^{T} \mathbf{C} \mathbf{r}| + \frac{1}{2} | \mathbf{r}^{T} \mathbf{M} \mathbf{r} | > \\ \frac{1}{2} | \mathbf{r}^{T} \mathbf{D} \mathbf{r} | = \frac{1}{2} | d(r_{h}, r_{h}) |. \quad \end{aligned} $$

(50)

Therefore, property (2) leads to a contradiction

$$ 0=| \mathbf{r}^{T} \mathbf{S} \mathbf{r} | > \frac{1}{2} \delta \| r_{h} \|^{2}_{Q} >0,$$

which proves the invertibility of S. This ensures the uniqueness of p, while (48) ensures that of u. The uniqueness of p ensures that the FE problem (46) does not present spurious modes.