Finite Element Solvers for Biot’s Poroelasticity Equations in Porous Media

Abstract

We study and compare five different combinations of finite element spaces for approximating the coupled flow and solid deformation system, so-called Biot’s equations. The permeability and porosity fields are heterogeneous and depend on solid displacement and fluid pressure. We provide detailed comparisons among the continuous Galerkin, discontinuous Galerkin, enriched Galerkin, and two types of mixed finite element methods. Several advantages and disadvantages for each of the above techniques are investigated by comparing local mass conservation properties, the accuracy of the flux approximation, number of degrees of freedom (DOF), and wall and CPU times. Three-field formulation methods with fluid velocity as an additional primary variable generally require a larger number of DOF, longer wall and CPU times, and a greater number of iterations in the linear solver in order to converge. The two-field formulation, a combination of continuous and enriched Galerkin function space, requires the fewest DOF among the methods that conserve local mass. Moreover, our results illustrate that three out of the five methods conserve local mass and produce similar flux approximations when conductivity alteration is included. These comparisons of the key performance indicators of different combinations of finite element methods can be utilized to choose the preferred method based on the required accuracy and the available computational resources.

Appendices

Appendices

In these appendices, we describe the fully discretized formulation of each method listed in Table 1.

Continuous Galerkin (CG) Method

We begin by defining the finite element space for the continuous Galerkin (CG) method for a vector-valued function

$$\begin{aligned} {\mathscr {U}}_{h}^{{\mathrm {CG}}_{k}}\left( {\mathscr {T}}_{h}\right) :=\left\{ \varvec{\psi _u} \in {\mathbb {C}}^{0}({\varOmega }{; {\mathbb {R}}^d}) :\left. \varvec{\psi _u}\right| _{T} \in {\mathbb {Q}}_{k}(T{; {\mathbb {R}}^d}), \forall T \in {\mathscr {T}}_{h}\right\} , \end{aligned}$$

(37)

where k indicates the order of polynomials, ${\mathbb {C}}^0({\varOmega }{; {\mathbb {R}}^d})$ denotes the space of vector-valued piecewise continuous polynomials, and ${\mathbb {Q}}_{k}(T{; {\mathbb {R}}^d})$ is the space of polynomials of degree at most k over each element T. For our problem, the vector-valued displacement field is approximated using CG space. Next, the CG space for scalar-valued functions is defined as

$$\begin{aligned} {\mathscr {P}}_{h}^{{\mathrm {CG}}_{k}}\left( {\mathscr {T}}_{h}\right) :=\left\{ \psi _p \in {\mathbb {C}}^{0}({\varOmega }) :\left. \psi _p\right| _{T} \in {\mathbb {Q}}_{k}(T), \forall T \in {\mathscr {T}}_{h}\right\} , \end{aligned}$$

(38)

for the pressure solution.

Thus, we seek the approximated displacement solution ($\varvec{u_h}$) by discretizing the linear momentum balance Eq. (6) employing the above CG finite element spaces as done in Choo and Lee (2018), Kadeethum et al. (2019a) and Vik et al. (2018). The fully discretized linear momentum balance Eq. (6) can be defined using the following forms

$$\begin{aligned} {\mathscr {A}}_u\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \varvec{\psi _u} \right) = {\mathscr {L}}_u\left( \varvec{\psi _u} \right) , \quad \forall \varvec{\psi _u} \in {\mathscr {U}}_{h}^{{\mathrm {CG}}_{2}}\left( {\mathscr {T}}_{h}\right) , \end{aligned}$$

(39)

at each time step $t^n$, where

$$\begin{aligned} {\mathscr {A}}_u\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \varvec{\psi _u} \right) = \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \varvec{\sigma }^{\prime }\left( \varvec{u}_{h}\right) : \nabla ^{s} \varvec{\psi _u} \, d V - \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \alpha p_{h} \varvec{I} : \nabla ^{s} \varvec{\psi _u} \, d V,\nonumber \\ \end{aligned}$$

(40)

and

$$\begin{aligned} {\mathscr {L}}_u\left( \varvec{\psi _u} \right) =\sum _{T \in {\mathscr {T}}_{h}} \int _{T} \varvec{f} \varvec{\psi _u} \, d V+\sum _{e \in {\mathscr {E}}_{h}^{N}} \int _{e} {{\varvec{\sigma }}_{D}} \varvec{\psi _u} \, d S, \quad \forall \varvec{\psi _u} \in {\mathscr {U}}_{h}^{{\mathrm {CG}}_{2}}\left( {\mathscr {T}}_{h}\right) . \end{aligned}$$

(41)

Here, $\nabla ^{s}$ is a symmetric gradient operator. We can split

$$\begin{aligned} {{\mathscr {A}}_u\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \varvec{\psi _u} \right) } := a\left( \varvec{u}_{h}^{n}, \varvec{\psi _u} \right) +b\left( p_{h}^{n}, \varvec{\psi _u} \right) , \end{aligned}$$

(42)

by using

$$\begin{aligned} a\left( \varvec{u}_{h}^{n}, \varvec{\psi _u} \right) := \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \varvec{\sigma }^{\prime }\left( \varvec{u}_{h}^n\right) : \nabla ^{s} \varvec{\psi _u} \, d V, \end{aligned}$$

(43)

$$\begin{aligned} b\left( p_{h}^{n}, \varvec{\psi _u} \right) := \sum _{T \in {\mathscr {T}}_{h}} \int _{T} -\alpha p_{h}^n \varvec{I} : \nabla ^{s} \varvec{\psi _u} \, d V. \end{aligned}$$

(44)

For the approximated pressure solution ($p_h$), we solve the mass balance Eq. (10) combined with Eqs. (11) and (25) as follows

$$\begin{aligned} {\mathscr {A}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) = {\mathscr {L}}_p\left( \psi _p; \varvec{u}_{h}^{n} \right) , \quad \forall \psi _p \in {\mathscr {P}}_{h}^{{\mathrm {CG}}_{1}}\left( {\mathscr {T}}_{h}\right) , \end{aligned}$$

(45)

for each time step $t^n$, where

$$\begin{aligned}&{\mathscr {A}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) \nonumber \\&\quad := \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \alpha \nabla \cdot {\varvec{u}}_{h}^n\psi _{p} \, d V + \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \left( \phi c_{f}+\frac{\alpha -\phi }{K_{s}}\right) p_{h}^n \psi _{p} \, d V \nonumber \\&\quad \quad + {\varDelta } t^n \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \varvec{\kappa }( \varvec{u}_{h}^{n})\left( \nabla p_{h}^n-\rho {\mathbf {g}}\right) \cdot \nabla \psi _{p} \, d V \nonumber \\&\quad \quad - {\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }( \varvec{u}_{h}^{n})\left( \nabla p_{h}^n-\rho {\mathbf {g}}\right) \right\} _{\delta _{e}} \cdot \llbracket \psi _p \rrbracket \, d S \nonumber \\&\quad \quad - {\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }( \varvec{u}_{h}^{n}) \nabla \psi _{p}\right\} _{\delta _{e}} \cdot \llbracket p_h^n \rrbracket \, d S \nonumber \\&\quad \quad + {\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e} \frac{\beta }{h_{e}} {\varvec{\kappa }( \varvec{u}_{h}^{n})}_{{e}} \llbracket p_h^n \rrbracket \cdot \llbracket \psi _p \rrbracket \, d S, \end{aligned}$$

(46)

$$\begin{aligned}&{\mathscr {L}}_p\left( \psi _p {; \varvec{u}_{h}^{n}} \right) \nonumber \\&\quad := \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \alpha \nabla \cdot {\varvec{u}}_{h}^{n-1}\psi _{p} \, d V+ \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \left( \phi c_{f}+\frac{\alpha -\phi }{K_{s}}\right) p_{h}^{n-1} \psi _{p} \, d V \nonumber \\&\quad \quad + {\varDelta } t^n \sum _{T \in {\mathscr {T}}_{h}} \int _{T} g \psi _{p} \, d V+ {\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{N}} \int _{e} q_{D} \psi _{p} \, d S \nonumber \\&\quad \quad -{\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e} \varvec{\kappa }( \varvec{u}_{h}^{n}) \nabla \psi _{p} \cdot p_{D} {\varvec{n}} \, d S \nonumber \\&\quad \quad + {\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e} \frac{\beta }{h_{e}} {\varvec{\kappa }( \varvec{u}_{h}^{n})}_{{e}} \llbracket \psi _p \rrbracket \cdot p_D \varvec{n} \, d S. \end{aligned}$$

(47)

Here, again, we can split

$$\begin{aligned} {{\mathscr {A}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) } := c\left( \varvec{u}_{h}^{n}, \psi _p \right) +d\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right) , \end{aligned}$$

(48)

with

$$\begin{aligned} c\left( \varvec{u}_{h}^{n}, \psi _p \right)&:= \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \alpha \nabla \cdot {\varvec{u}}_{h}^n\psi _{p} \, d V, \end{aligned}$$

(49)

$$\begin{aligned} d\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right)&:= \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \left( \phi c_{f}+\frac{\alpha -\phi }{K_{s}}\right) p_{h}^n \psi _{p} \, d V \nonumber \\&\qquad + {\varDelta } t^n \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \varvec{\kappa }( \varvec{u}_{h}^{n})\left( \nabla p_{h}^n-\rho {\mathbf {g}}\right) \cdot \nabla \psi _{p} \, d V \nonumber \\&\qquad - {\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }( \varvec{u}_{h}^{n})\left( \nabla p_{h}^n-\rho {\mathbf {g}}\right) \right\} _{\delta _{e}} \cdot \llbracket \psi _p \rrbracket \, d S \nonumber \\&\qquad - {\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }( \varvec{u}_{h}^{n}) \nabla \psi _{p}\right\} _{\delta _{e}} \cdot \llbracket p_h^n \rrbracket \, d S \nonumber \\&\qquad + {\varDelta } t^n \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e} \frac{\beta }{h_{e}} {\varvec{\kappa }( \varvec{u}_{h}^{n})}_{{e}} \llbracket p_h^n \rrbracket \cdot \llbracket \psi _p \rrbracket \, d S. \end{aligned}$$

(50)

We note that the symmetric internal penalty approach is employed to weakly impose the Dirichlet boundary condition (Bazilevs and Hughes 2007; Nitsche 1971; Hansbo 2005). The interior penalty parameter $\beta$ and $h_e$ characteristic length are calculated as

$$\begin{aligned} h_{e} :=\frac{{\text {meas}}\left( T^{+}\right) +{\text {meas}}\left( T^{-}\right) }{2 {\text {meas}}(e)}. \end{aligned}$$

(51)

The system of equations is nonlinear (the nonlinear variable being reported after a semicolon for the sake of clarity). Thus, Newton’s method is employed in this work. The Jacobian matrix that arises from Eqs. (39) and (45) is

$$\begin{aligned} \left[ \begin{array}{c@{\quad }c}{{\mathscr {J}}_{ {uu }}} &{} {{\mathscr {J}}_{u p}} \\ {{\mathscr {J}}_{ {pu }}} &{} {{\mathscr {J}}_{ {pp }}} \end{array}\right] \left\{ \begin{array}{l}{\delta \varvec{u}_{h}} \\ {\delta p_{h}} \end{array}\right\} =-\left\{ \begin{array}{l}{R_{u}} \\ {R_{p}} \end{array}\right\} , \end{aligned}$$

(52)

where

$$\begin{aligned} \begin{aligned}&{\mathscr {J}}_{u u}:=\frac{\partial a\left( \varvec{u}_{h}^{n}, \varvec{\psi _u} \right) }{\partial \varvec{u}_{h}^{n}}, \quad {\mathscr {J}}_{u p}:=\frac{\partial b\left( p_{h}^{n}, \varvec{\psi _u} \right) }{\partial p_{h}^{n}}, \\&{\mathscr {J}}_{p u}:=\frac{\partial c\left( \varvec{u}_{h}^{n}, \psi _p \right) }{\partial \varvec{u}_{h}^{n}}, \quad {\mathscr {J}}_{p p}:=\frac{\partial d\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right) }{\partial p_{h}^{n}}, \end{aligned} \end{aligned}$$

(53)

and the residual vector is defined as

$$\begin{aligned} \begin{array}{l} {R_{u}:=a\left( \varvec{u}_{h}^{n}, \varvec{\psi _u} \right) +b\left( p_{h}^{n}, \varvec{\psi _u} \right) } - {\mathscr {L}}_u\left( \varvec{\psi _u} \right) , \\ {R_{p}:=c\left( \varvec{u}_{h}^{n}, \psi _p \right) +d\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right) - {\mathscr {L}}_p\left( \psi _p {; \varvec{u}_{h}^{n}} \right) . } \end{array} \end{aligned}$$

(54)

Here, $\delta \varvec{u}_{h}$ and $\delta p_{h}$ are Newton increments of the $\varvec{u}_h$ and $p_h$, respectively. The block structure for the CG method arises as follows

$$\begin{aligned} \left[ \begin{array}{l@{\quad }l}{{\mathscr {J}}_{u u}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{up}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_1}} \\ {{\mathscr {J}}_{pu}^{{\mathrm {CG}}_1 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{p p}^{{\mathrm {CG}}_1 \times {\mathrm {CG}}_1}} \end{array}\right] \left\{ \begin{array}{l}{\left( \delta \varvec{u}_{h}^{n}\right) ^{{\mathrm {CG}}_2}} \\ {\left( \delta p_{h}^{n}\right) ^{{\mathrm {CG}}_1}} \end{array}\right\} = - \left\{ \begin{array}{l}{R_{u}^{ {\mathrm {CG}}_2}} \\ {R_{p}^{{\mathrm {CG}}_1}} \end{array}\right\} . \end{aligned}$$

(55)

Discontinuous Galerkin (DG) Method

For the discontinuous Galerkin (DG) method presented in Table 1, we approximate the $\varvec{u}_h$ using the CG finite element space, Eqs. (37) and (39), but we utilize the DG finite element space to approximate the $p_h$. The DG finite element space with the polynomial order k is given as

$$\begin{aligned} {\mathscr {P}}_{h}^{{\mathrm {DG}}_{k}}\left( {\mathscr {T}}_{h}\right) :=\left\{ \psi _p \in L^{2}({\varOmega }) :\left. \psi _p\right| _{T} \in {\mathbb {Q}}_{k}(T), \forall T \in {\mathscr {T}}_{h}\right\} , \end{aligned}$$

(56)

where $L^{2}({\varOmega })$ is the space of the square integrable functions. This nonconforming finite element space allows us to consider discontinuous coefficients and preserves the local mass conservation. We then discretize Eqs. (10), (11), and (25) as

$$\begin{aligned} {\mathscr {B}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) = {\mathscr {L}}_p\left( \psi _p; \varvec{u}_{h}^{n} \right) , \quad \forall \psi _p \in {\mathscr {P}}_{h}^{{\mathrm {DG}}_{1}}\left( {\mathscr {T}}_{h}\right) , \end{aligned}$$

(57)

for each time step $t^n$, where

$$\begin{aligned} \begin{aligned}&{\mathscr {B}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) \\&\quad := \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \alpha \nabla \cdot {\varvec{u}}_{h}^n\psi _{p} \, d V + \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \left( \phi c_{f}+\frac{\alpha -\phi }{K_{s}}\right) p_{h}^n \psi _{p} \, d V \\&\quad \quad + {\varDelta } t^n \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \varvec{\kappa }( \varvec{u}_{h}^{n})\left( \nabla p_{h}^n-\rho {\mathbf {g}}\right) \cdot \nabla \psi _{p} \, d V \\&\quad \quad - {\varDelta } t^n \sum _{e \in {\mathscr {E}}_h^{I} \cup {\mathscr {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }( \varvec{u}_{h}^{n})\left( \nabla p_{h}^n-\rho {\mathbf {g}}\right) \right\} _{\delta _{e}} \cdot \llbracket \psi _p \rrbracket \, d S \\&\quad \quad - {\varDelta } t^n \sum _{e \in {\mathscr {E}}_h^{I} \cup {\mathscr {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }( \varvec{u}_{h}^{n}) \nabla \psi _{p}\right\} _{\delta _{e}} \cdot \llbracket p_h^n \rrbracket \, d S \\&\quad \quad + {\varDelta } t^n \sum _{e \in {\mathscr {E}}_h^{I} \cup {\mathscr {E}}_{h}^{D}} \int _{e} \frac{\beta }{h_{e}} {\varvec{\kappa }( \varvec{u}_{h}^{n})}_{{e}} \llbracket p_h^n \rrbracket \cdot \llbracket \psi _p \rrbracket \, d S. \end{aligned} \end{aligned}$$

(58)

The ${\mathscr {L}}_p\left( \psi _p; \varvec{u}_{h}^{n} \right)$ term is similar to Eq. (47) as presented in the previous section. Next, we split

$$\begin{aligned} {{\mathscr {B}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) } := c\left( \varvec{u}_{h}^{n}, \psi _p \right) +e\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right) . \end{aligned}$$

(59)

We utilize the same definition for $c\left( \varvec{u}_{h}^{n}, \psi _p \right)$ from the CG method, Eq. (49), but we define $e\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right)$ as

$$\begin{aligned} \begin{aligned} e\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right)&:= \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \left( \phi c_{f}+\frac{\alpha -\phi }{K_{s}}\right) p_{h}^n \psi _{p} \, d V \\&\quad + {\varDelta } t^n \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \varvec{\kappa }( \varvec{u}_{h}^{n})\left( \nabla p_{h}^n-\rho {\mathbf {g}}\right) \cdot \nabla \psi _{p} \, d V \\&\quad - {\varDelta } t^n \sum _{e \in {\mathscr {E}}_h^{I} \cup {\mathscr {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }( \varvec{u}_{h}^{n})\left( \nabla p_{h}^n-\rho {\mathbf {g}}\right) \right\} _{\delta _{e}} \cdot \llbracket \psi _p \rrbracket \, d S \\&\quad - {\varDelta } t^n \sum _{e \in {\mathscr {E}}_h^{I} \cup {\mathscr {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }( \varvec{u}_{h}^{n}) \nabla \psi _{p}\right\} _{\delta _{e}} \cdot \llbracket p_h^n \rrbracket \, d S \\&\quad + {\varDelta } t^n \sum _{e \in {\mathscr {E}}_h^{I} \cup {\mathscr {E}}_{h}^{D}} \int _{e} \frac{\beta }{h_{e}} {\varvec{\kappa }( \varvec{u}_{h}^{n})}_{{e}} \llbracket p_h^n \rrbracket \cdot \llbracket \psi _p \rrbracket \, d S. \end{aligned} \end{aligned}$$

(60)

Compared to ${{\mathscr {A}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) }$ from the CG method, we note that there are additional surface integral terms $\int _{e}\cdot \, d S$ on the interior faces ($e \in {\mathscr {E}}_h^{I}$) in ${{\mathscr {B}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) }$. This results in more demanding computational time when the matrix is assembled.

Again, we note that this system of equations are nonlinear forms in the displacement $\varvec{u}_{h}^{n}$. We form the Jacobian matrix that arises from Eqs. (39) and (57) as follows

$$\begin{aligned} \left[ \begin{array}{c@{\quad }c}{{\mathscr {J}}_{ {uu }}} &{} {{\mathscr {J}}_{u p}} \\ {{\mathscr {J}}_{ {pu }}} &{} {{\mathscr {J}}_{ {pp }}} \end{array}\right] \left\{ \begin{array}{l}{\delta \varvec{u}_{h}} \\ {\delta p_{h}} \end{array}\right\} =-\left\{ \begin{array}{l}{R_{u}} \\ {R_{p}} \end{array}\right\} , \end{aligned}$$

(61)

where

$$\begin{aligned} \begin{aligned}&{\mathscr {J}}_{u u}:=\frac{\partial a\left( \varvec{u}_{h}^{n}, \varvec{\psi _u} \right) }{\partial \varvec{u}_{h}^{n}}, \quad {\mathscr {J}}_{u p}:=\frac{\partial b\left( p_{h}^{n}, \varvec{\psi _u} \right) }{\partial p_{h}^{n}}, \\&{\mathscr {J}}_{p u}:=\frac{\partial c\left( \varvec{u}_{h}^{n}, \psi _p \right) }{\partial \varvec{u}_{h}^{n}}, \quad {\mathscr {J}}_{p p}:=\frac{\partial e\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right) }{\partial p_{h}^{n}}, \end{aligned} \end{aligned}$$

(62)

and the residual vector is defined as

$$\begin{aligned} \begin{array}{l} {R_{u}:=a\left( \varvec{u}_{h}^{n}, \varvec{\psi _u} \right) +b\left( p_{h}^{n}, \varvec{\psi _u} \right) } - {\mathscr {L}}_u\left( \varvec{\psi _u} \right) , \\ {R_{p}:=c\left( \varvec{u}_{h}^{n}, \psi _p \right) +e\left( p_{h}^{n}, \psi _p{; \varvec{u}_{h}^{n}} \right) - {\mathscr {L}}_p\left( \psi _p {; \varvec{u}_{h}^{n}} \right) , } \end{array} \end{aligned}$$

(63)

at each time step $t^n$. The block structure for the DG method arises as follows

$$\begin{aligned} \left[ \begin{array}{l@{\quad }l} {{\mathscr {J}}_{u u}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{up}^{{\mathrm {CG}}_2 \times {\mathrm {DG}}_1}} \\ {{\mathscr {J}}_{pu}^{{\mathrm {DG}}_1 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{p p}^{{\mathrm {DG}}_1 \times {\mathrm {DG}}_1}} \end{array}\right] \left\{ \begin{array}{l} {\left( \delta \varvec{u}_{h}^{n}\right) ^{{\mathrm {CG}}_2}} \\ {\left( \delta p_{h}^{n}\right) ^{{\mathrm {DG}}_1}} \end{array}\right\} = - \left\{ \begin{array}{l}{R_{u}^{ {\mathrm {CG}}_2}} \\ {R_{p}^{{\mathrm {DG}}_1}} \end{array}\right\} . \end{aligned}$$

(64)

Enriched Galerkin (EG) Method

For the EG method presented in Table 1, we define the EG finite element space with the polynomial order k as

$$\begin{aligned} {\mathscr {P}}_{h}^{{\mathrm {EG}}_{k}}\left( {\mathscr {T}}_{h}\right) :={\mathscr {P}}_{h}^{{\mathrm {CG}}_{k}}\left( {\mathscr {T}}_{h}\right) +{\mathscr {P}}_{h}^{{\mathrm {DG}}_{0}}\left( {\mathscr {T}}_{h}\right) , \end{aligned}$$

(65)

where the CG finite element space is enriched by ${\mathscr {P}}_{h}^{{\mathrm {DG}}_{0}}\left( {\mathscr {T}}_{h}\right)$, the piecewise constant functions. As in the previous section, we approximate the displacement ($\varvec{u}_h$) using the CG finite element space Eq. (39), but here the $p_h$ is approximated by the EG finite element space Eq. (65). Thus, we seek $p_h$ by solving

$$\begin{aligned} {\mathscr {B}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right) = {\mathscr {L}}_p\left( \psi _p; \varvec{u}_{h}^{n} \right) , \quad \forall \psi _p \in {\mathscr {P}}_{h}^{{\mathrm {EG}}_{1}}\left( {\mathscr {T}}_{h}\right) . \end{aligned}$$

(66)

We note that the EG method employs the same bilinear ${\mathscr {B}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right)$ and linear ${\mathscr {L}}_p\left( \psi _p; \varvec{u}_{h}^{n} \right)$ forms as the DG method. Since the EG space Eq. (65) only takes piecewise constants (DG$_0$) to enforce the discontinuity, several terms involving the derivatives in ${\mathscr {B}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right)$ and ${\mathscr {L}}_p\left( \psi _p; \varvec{u}_{h}^{n} \right)$ become zero. Thus, not only does the EG method have fewer degrees of freedom than DG, but it requires less computational time for assembling the matrix.

We again note that this system of equations are nonlinear forms in the displacement $\varvec{u}_{h}^{n}$. We employ the same Jacobian matrix, Eqs. (61) and (62), and residual vector, Eq. (63), built for the DG method. The block structure used for the DG method Eq. (64), however, is decomposed into

$$\begin{aligned} \left[ \begin{array}{l@{\quad }l@{\quad }l} {{\mathscr {J}}_{u u}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{u p}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_1}} &{} {{\mathscr {J}}_{u p}^{{\mathrm {CG}}_2 \times {\mathrm {DG}}_0}} \\ {{\mathscr {J}}_{p u}^{{\mathrm {CG}}_1 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{p p}^{{\mathrm {CG}}_1 \times {\mathrm {CG}}_1}} &{} {{\mathscr {J}}_{p p}^{{\mathrm {CG}}_1 \times {\mathrm {DG}}_0}} \\ {{\mathscr {J}}_{p u}^{{\mathrm {DG}}_0 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{p p}^{{\mathrm {DG}}_0 \times {\mathrm {CG}}_1}} &{} {{\mathscr {J}}_{p p}^{{\mathrm {DG}}_0 \times {\mathrm {DG}}_0}} \end{array} \right] \left\{ \begin{array}{l} {\left( \delta \varvec{u}_{h}^{n}\right) ^{{\mathrm {CG}}_2}} \\ {\left( \delta {p_{h}}^{n}\right) ^{{\mathrm {CG}}_1}}\\ {\left( \delta {p_{h}}^{n}\right) ^{{\mathrm {DG}}_0}} \end{array}\right\} = - \left\{ \begin{array}{l} {{R}_u^{{\mathrm {CG}}_2}} \\ {{R}_p^{{\mathrm {CG}}_1}} \\ {{R}_p^{{\mathrm {DG}}_0}} \end{array}\right\} ,\nonumber \\ \end{aligned}$$

(67)

to satisfy Eq. (66).

Raviart–Thomas Mixed Finite Element (MFE-RT) Method

For the MFE-RT method presented in Table 1, we also approximate the $\varvec{u}_h$ using the CG finite element space [Eqs. (37) and (39)]. The mass balance Eq. (10) for the MFE-RT method is discretized by using the piecewise constant (i.e., DG finite element space Eq. (56) with $k=0$).

Thus, we seek $p_h$ by solving

$$\begin{aligned} {{\mathscr {C}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n},\varvec{v}_{h}^{n}), \psi _p \right) } = {\mathscr {M}}_p\left( \psi _p \right) , \quad \forall \psi _p \in {\mathscr {P}}_{h}^{{\mathrm {DG}}_{0}}\left( {\mathscr {T}}_{h}\right) , \end{aligned}$$

(68)

for each time step $t^n$, where

$$\begin{aligned} \begin{aligned} {\mathscr {C}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n},\varvec{v}_{h}^{n}), \psi _p; \varvec{u}_{h}^{n} \right)&:= \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \alpha \nabla \cdot {\varvec{u}}_{h}^n\psi _{p} \, d V \\&\quad + \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \left( \phi c_{f}+\frac{\alpha -\phi }{K_{s}}\right) p_{h}^n \psi _{p} \, d V \\&\quad + {\varDelta } t^n \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \nabla \cdot \left( \rho \varvec{v}_{h}^{n} \right) \psi _{p} \, d V, \end{aligned} \end{aligned}$$

(69)

and

$$\begin{aligned} \begin{aligned} {\mathscr {M}}_p\left( \psi _p {; \varvec{u}_{h}^{n}} \right)&:= \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \alpha \nabla \cdot {\varvec{u}}_{h}^{n-1}\psi _{p} \, d V \\&\quad + \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \left( \phi c_{f}+\frac{\alpha -\phi }{K_{s}}\right) p_{h}^{n-1} \psi _{p} \, d V \\&\quad + {\varDelta } t^n \sum _{T \in {\mathscr {T}}_{h}} \int _{T} g \psi _{p} \, d V. \end{aligned} \end{aligned}$$

(70)

Here, we can split Eq. (69) as

$$\begin{aligned} {{\mathscr {C}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n},\varvec{v}_{h}^{n}), \psi _p \right) } := f\left( \varvec{u}_{h}^{n}, \psi _p \right) +g\left( p_{h}^{n}, \psi _p \right) +h\left( \varvec{v}_{h}^{n}, \psi _p \right) , \end{aligned}$$

(71)

where

$$\begin{aligned} f\left( \varvec{u}_{h}^{n}, \psi _p \right):= & {} \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \alpha \nabla \cdot {\varvec{u}}_{h}^n\psi _{p} \, d V, \end{aligned}$$

(72)

$$\begin{aligned} g\left( p_{h}^{n}, \psi _p \right):= & {} \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \left( \phi c_{f}+\frac{\alpha -\phi }{K_{s}}\right) p_{h}^n \psi _{p} \, d V, \end{aligned}$$

(73)

and

$$\begin{aligned} \begin{aligned} h\left( \varvec{v}_{h}^{n}, \psi _p \right)&:= {\varDelta } t^n \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \nabla \cdot \left( \rho \varvec{v}_{h}^{n} \right) \psi _{p} \, d V. \end{aligned} \end{aligned}$$

(74)

For the velocity approximation ($\varvec{v}_h$), we use the RT finite element space, which is defined as

$$\begin{aligned} {\mathscr {V}}_{h}^{{\mathrm {RT}}_{k}}\left( {\mathscr {T}}_{h}\right) :=\left\{ \varvec{\psi _v} \in {\mathbb {C}}^{0}({\varOmega }) :\left. \varvec{\psi _v}\right| _{T} \in {\mathbb {Q}}_{k}(T)^m + x{\mathbb {Q}}_{k}(T), \forall T \in {\mathscr {T}}_{h}\right\} , \end{aligned}$$

(75)

${\mathbf {q}}(x) =\left( \begin{array}{c}{q_{1}(x)} \\ {q_{2}(x)} \\ {\vdots } \\ {q_{m}(x)} \end{array}\right) +q_{0}(x) \left( \begin{array}{c}{x_{1}} \\ {x_{2}} \\ {\vdots } \\ {x_{m}} \end{array}\right) \quad \forall x:=\left( x_{1}, x_{2}, \ldots , x_{m}\right) ^{{\mathrm {T}}} \in {\varOmega }$, ${\mathbf {q}}(x) \subset {\mathbb {Q}}_{k}(T)^m$, and $q_{0}(x) \subset {\mathbb {Q}}_{k}(T)$. We then solve the following Darcy velocity Eq. (11) to obtain $\varvec{v}_{h}^{n}$ by

$$\begin{aligned} {\mathscr {C}}_v\left( (\varvec{v}_{h}^{n}, p_{h}^{n}),\varvec{\psi }_{v}; \varvec{u}_{h}^{n} \right) = {\mathscr {M}}_v\left( \varvec{\psi _v} \right) , \quad \forall \varvec{\psi }_{v} \in {\mathscr {V}}_{h}^{{\mathrm {RT}}_{1}}\left( {\mathscr {T}}_{h}\right) , \end{aligned}$$

(76)

where

$$\begin{aligned} {\mathscr {C}}_v\left( (\varvec{v}_{h}^{n}, p_{h}^{n}),\varvec{\psi }_{v}; \varvec{u}_{h}^{n} \right) := \sum _{T \in {\mathscr {T}}_{h}} \int _{T} p_{h}^{n} \nabla \cdot \varvec{\psi _v} \, d V+\sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \varvec{\kappa }( \varvec{u}_{h}^{n})^{-1} \varvec{v}_{h}^{n} \varvec{\psi _v} \, d V, \nonumber \\ \end{aligned}$$

(77)

and

$$\begin{aligned} \begin{aligned} {\mathscr {M}}_v\left( \varvec{\psi _v} \right) :=&- \sum _{e \in {\mathscr {E}}_{h}^{D}} \int _{e} p_{D} \varvec{\psi }_{v} \cdot \varvec{n} \, d S. \end{aligned} \end{aligned}$$

(78)

Here, we can split

$$\begin{aligned} {{\mathscr {C}}_v\left( (\varvec{v}_{h}^{n}, p_{h}^{n},), \varvec{\psi }_v {; \varvec{u}_{h}^{n}} \right) } := i\left( p_{h}^{n}, \varvec{\psi }_v \right) + j\left( \varvec{v}_{h}^{n}, \varvec{\psi }_v{; \varvec{u}_{h}^{n}} \right) , \end{aligned}$$

(79)

where

$$\begin{aligned} i\left( p_{h}^{n}, \varvec{\psi }_v \right) := \sum _{T \in {\mathscr {T}}_{h}} \int _{T} p_{h}^{n} \nabla \cdot \varvec{\psi _v} \, d V, \end{aligned}$$

(80)

and

$$\begin{aligned} j\left( \varvec{v}_{h}^{n}, \varvec{\psi }_v{; \varvec{u}_{h}^{n}} \right) := \sum _{T \in {\mathscr {T}}_{h}} \int _{T} \rho \varvec{\kappa }( \varvec{u}_{h}^{n})^{-1} \varvec{v}_{h}^{n} \varvec{\psi _v} \, d V. \end{aligned}$$

(81)

Similar to the CG, DG, and EG methods, we note that this system of equations are nonlinear forms for the displacement $\varvec{u}_{h}^{n}$. We then form the Jacobian matrix that arises from Eqs. (39), (68), and (76) as follows

$$\begin{aligned} \left[ \begin{array}{c@{\quad }c@{\quad }c}{{\mathscr {J}}_{ {uu }}} &{} {{\mathscr {J}}_{up}} &{} {{\mathscr {J}}_{u v}} \\ {{\mathscr {J}}_{ {pu }}} &{} {{\mathscr {J}}_{ {pp }}} &{} {{\mathscr {J}}_{p v}} \\ {{\mathscr {J}}_{ {vu }}} &{} {{\mathscr {J}}_{vp}} &{} {{\mathscr {J}}_{v v}} \end{array}\right] \left\{ \begin{array}{l}{\delta \varvec{u}_{h}} \\ {\delta p_{h}} \\ {\delta \varvec{v}_{h}} \end{array}\right\} =-\left\{ \begin{array}{l}{R_{u}} \\ {R_{p}} \\ {R_{v}} \end{array}\right\} , \end{aligned}$$

(82)

where $\delta \varvec{v}_{h}$ is the Newton increments of the $\varvec{v}_h$ and

$$\begin{aligned} \begin{aligned}&{\mathscr {J}}_{u u}:=\frac{\partial a\left( \varvec{u}_{h}^{n}, \varvec{\psi _u} \right) }{\partial \varvec{u}_{h}^{n}}, \quad {\mathscr {J}}_{u p}:=\frac{\partial b\left( p_{h}^{n}, \varvec{\psi _u} \right) }{\partial p_{h}^{n}}, \quad {\mathscr {J}}_{u v}:=0, \\&{\mathscr {J}}_{p u}:=\frac{\partial f\left( \varvec{u}_{h}^{n}, \psi _p \right) }{\partial \varvec{u}_{h}^{n}}, \quad {\mathscr {J}}_{p p}:=\frac{\partial g\left( p_{h}^{n}, \psi _p \right) }{\partial p_{h}^{n}}, \quad {\mathscr {J}}_{p v}:=\frac{\partial h\left( \varvec{v}_{h}^{n}, \psi _p \right) }{\partial \varvec{v}_{h}^{n}}, \\&{\mathscr {J}}_{v u}:=0, \quad {\mathscr {J}}_{v p}:=\frac{\partial i\left( p_{h}^{n}, \varvec{\psi }_v \right) }{\partial p_{h}^{n}}, \quad {\mathscr {J}}_{v v}:=\frac{\partial j\left( \varvec{v}_{h}^{n}, \varvec{\psi }_v{; \varvec{u}_{h}^{n}} \right) }{\partial \varvec{v}_{h}^{n}}, \end{aligned} \end{aligned}$$

(83)

and the residual vector is defined as

$$\begin{aligned} \begin{array}{l} {R_{u}:=a\left( \varvec{u}_{h}^{n}, \varvec{\psi _u} \right) +b\left( p_{h}^{n}, \varvec{\psi _u} \right) } - {\mathscr {L}}_u\left( \varvec{\psi _u} \right) ,\\ {R_{p}:=f\left( \varvec{u}_{h}^{n}, \psi _p \right) + g\left( p_{h}^{n}, \psi _p \right) +h\left( \varvec{v}_{h}^{n}, \psi _p \right) - {\mathscr {M}}_p\left( \psi _p {; \varvec{u}_{h}^{n}} \right) , } \\ {R_{v}:= i\left( p_{h}^{n}, \varvec{\psi }_v \right) +j\left( \varvec{v}_{h}^{n}, \varvec{\psi }_v{; \varvec{u}_{h}^{n}} \right) - {\mathscr {M}}_v\left( \varvec{\psi _v} \right) .} \end{array} \end{aligned}$$

(84)

The block structure of the MFE-RT method arises as

$$\begin{aligned} \left[ \begin{array}{c@{\quad }c@{\quad }c} {{\mathscr {J}}_{u u}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{u p}^{{\mathrm {CG}}_2 \times {\mathrm {DG}}_0}} &{} 0 \\ {{\mathscr {J}}_{p u}^{{\mathrm {CG}}_2 \times {\mathrm {DG}}_0}} &{} {{\mathscr {J}}_{p p}^{{\mathrm {DG}}_0 \times {\mathrm {DG}}_0}} &{} {{\mathscr {J}}_{p v}^{{\mathrm {DG}}_0 \times {\mathrm {RT}}_1}} \\ 0 &{} {{\mathscr {J}}_{v p}^{{\mathrm {RT}}_1 \times {\mathrm {DG}}_0}} &{} {{\mathscr {J}}_{v v}^{{\mathrm {RT}}_1 \times {\mathrm {RT}}_1}} \end{array} \right] \left\{ \begin{array}{l}{\left( \delta \varvec{u}_{h}^{n}\right) ^{{\mathrm {CG}}_2}} \\ {\left( \delta {p_{h}}^{n}\right) ^{{\mathrm {DG}}_0}} \\ {\left( \delta {\varvec{v}_{h}}^{n}\right) ^{{\mathrm {RT}}_1}} \end{array}\right\} = - \left\{ \begin{array}{l} {{R}_u^{{\mathrm {CG}}_2}} \\ {{R}_p^{{\mathrm {DG}}_0}} \\ {{R}_v^{{\mathrm {RT}}_1}} \end{array}\right\} . \end{aligned}$$

(85)

Continuous Mixed Finite Element (MFE-P2) Method

For the MFE-P2 method presented in Table 1, we utilize a similar system of equations as employed by the MFE-RT method. The displacement $\varvec{u}_h$ is approximated by the CG finite element space, which is the same as the CG method Eq. (39). The $p_h$ in this method is approximated by using the CG finite element space Eq. (38) with $k=1$

$$\begin{aligned} {\mathscr {C}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n},\varvec{v}_{h}^{n}), \psi _p \right) = {\mathscr {M}}_p\left( \psi _p \right) , \quad \forall \psi _p \in {\mathscr {P}}_{h}^{{\mathrm {CG}}_{1}}\left( {\mathscr {T}}_{h}\right) . \end{aligned}$$

(86)

Here, the approximated velocity $\varvec{v}_h$ is also discretized by the vector-valued CG finite element space Eq. (37) but with higher-order $k=2$. Thus, we obtain the following for $\varvec{v}_h$

$$\begin{aligned} {\mathscr {C}}_v\left( (\varvec{v}_{h}^{n}, p_{h}^{n}),\varvec{\psi }_{v}; \varvec{u}_{h}^{n} \right) = {\mathscr {M}}_v\left( \varvec{\psi _v} \right) , \quad \forall \varvec{\psi }_{v} \in {\mathscr {U}}_{h}^{{\mathrm {CG}}_{2}}\left( {\mathscr {T}}_{h}\right) . \end{aligned}$$

(87)

As a result, the block structure of the MFE-RT method Eq. (85) is modified to

$$\begin{aligned} \left[ \begin{array}{c@{\quad }c@{\quad }c} {{\mathscr {J}}_{u u}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{u p}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_1}} &{} 0 \\ {{\mathscr {J}}_{p u}^{{\mathrm {CG}}_1 \times {\mathrm {CG}}_2}} &{} {{\mathscr {J}}_{p p}^{{\mathrm {CG}}_1 \times {\mathrm {CG}}_1}} &{} {{\mathscr {J}}_{p v}^{{\mathrm {CG}}_1 \times {\mathrm {CG}}_2}} \\ 0 &{} {{\mathscr {J}}_{v p}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_1}} &{} {{\mathscr {J}}_{v v}^{{\mathrm {CG}}_2 \times {\mathrm {CG}}_2}} \end{array} \right] \left\{ \begin{array}{l}{\left( \delta \varvec{u}_{h}^{n}\right) ^{{\mathrm {CG}}_2}} \\ {\left( \delta {p_{h}}^{n}\right) ^{{\mathrm {CG}}_1}} \\ {\left( \delta {\varvec{v}_{h}}^{n}\right) ^{{\mathrm {CG}}_2}} \end{array}\right\} = - \left\{ \begin{array}{l} {{R}_u^{{\mathrm {CG}}_2}} \\ {{R}_p^{{\mathrm {CG}}_1}} \\ {{R}_v^{{\mathrm {CG}}_2}} \end{array}\right\} .\nonumber \\ \end{aligned}$$

(88)

Remark 2

For both two- and three-field formulations, the Neumann boundary condition is naturally applied on the boundary faces, $e \in {\mathscr {E}}_{h}^{N}$. For the two-field formulation, the Dirichlet boundary condition is weakly enforced on $e \in {\mathscr {E}}_{h}^{D}$, for all CG, EG, and DG methods. On the other hand, the three-field formulation strongly applies the Dirichlet boundary conditions.