A multi-layer reduced model for flow in porous media with a fault and surrounding damage zones

TL;DR

This paper introduces a multi-layer reduced model for simulating fluid flow in porous media with complex fault structures, including damage zones, improving accuracy in heterogeneous and high-contrast permeability scenarios.

Contribution

The work presents a novel multi-layer reduced model that explicitly incorporates faults and surrounding damage zones as lower-dimensional features within porous media.

Findings

Model effectively captures high permeability contrast.

Simulation results demonstrate robustness across configurations.

Accurately represents fault and damage zone interactions.

Abstract

In this work we present a new conceptual model to describe fluid flow in a porous media system in presence of a large fault. Geological faults are often modeled simply as interfaces in the rock matrix, but they are complex structure where the high strain core is surrounded by the so called damage zones, characterized by the presence of smaller fractures which enhance the permeability of the medium. To obtain reliable simulation outcomes these damage zone, as well as the fault, have to be accurately described. The new model proposed in this work considers both these two regions as lower dimensional and embedded in the rock matrix. The model is presented, analyzed, and tested in several configurations to prove its robustness and ability to capture many important features, such as hight contrast and heterogeneity of permeability.

Tables

Table 1. Table 1: Errors corresponding to different values of the ϵ italic-ϵ \epsilon for the three different cases of the test in Subsection 5.1 .

$ϵ$	case (i)	case (ii)	case (iii)
${10}^{-2}$	1.12812e-2	9.72846e-3	8.93785e-3
$5\cdot {10}^{-3}$	7.29157e-3	6.83495e-3	6.36322e-3
$2.5\cdot {10}^{-3}$	5.7667e-3	5.24812e-3	5.21005e-3

Equations

$$\begin{array}{rlrl}& \begin{array}{rl}& {\sigma }_{{\rm Y}}^2{\mathbf{u}}_{{\rm Y}}+\nabla p_{{\rm Y}}=0\\ & \nabla \cdot {\mathbf{u}}_{{\rm Y}}+q_{{\rm Y}}=0\end{array}& & \text{in }{\rm Y}\\ & \mathrm{tr}{p}_{{\rm Y}}=\overline{p_{{\rm Y}}}& & \text{on }\partial {{\rm Y}}_p\\ & \mathrm{tr}{\mathbf{u}}_{{\rm Y}}\cdot \mathbf{\nu }=\overline{u_{{\rm Y}}}& & \text{on }\partial {{\rm Y}}_u\end{array},$$

$$\begin{array}{rlrl}& \begin{array}{rl}& \mathrm{tr}{\mathbf{u}}_{\Omega }\cdot \mathbf{n}=\mathrm{tr}{\mathbf{u}}_{\mu }\cdot \mathbf{n}\\ & \mathrm{tr}{p}_{\Omega }=\mathrm{tr}{p}_{\mu }\end{array}& & \text{on }M\\ & \begin{array}{rl}& \mathrm{tr}{\mathbf{u}}_{\mu }\cdot {\mathbf{n}}_{\mu }=\mathrm{tr}{\mathbf{u}}_{\gamma }\cdot {\mathbf{n}}_{\mu }\\ & \mathrm{tr}{p}_{\mu }=\mathrm{tr}{p}_{\gamma }\end{array}& & \text{on }\Gamma \end{array}.$$

$$p=(p_{\Omega },p_{\mu },p_{\gamma })\text{and}\mathbf{u}=({\mathbf{u}}_{\Omega },{\mathbf{u}}_{\mu },{\mathbf{u}}_{\gamma },u_{\Gamma }).$$

$$\begin{array}{rlrl}& \begin{array}{rl}& {\alpha }_{\Omega }^2{\mathbf{u}}_{\Omega }+\nabla p_{\Omega }=0\\ & {\nabla }_{\Omega }\cdot \mathbf{u}+q_{\Omega }=0\end{array}& & \text{in }\Omega \\ & \mathrm{tr}{p}_{\Omega }=\overline{p_{\Omega }}& & \text{on }\partial {\Omega }_p\\ & \mathrm{tr}{\mathbf{u}}_{\Omega }\cdot \mathbf{\nu }=\overline{u_{\Omega }}& & \text{on }\partial {\Omega }_u\end{array},$$

$$\begin{array}{rlrl}& \begin{array}{rl}& {\alpha }_{\mu }^2{\mathbf{u}}_{\mu }+\nabla p_{\mu }=0\\ & {\nabla }_{\mu }\cdot \mathbf{u}+q_{\mu }=0\end{array}& & \text{in }\mu \\ & \mathrm{tr}{p}_{\mu }=\overline{p_{\mu }}& & \text{on }\partial {\mu }_p\\ & \mathrm{tr}{\mathbf{u}}_{\mu }\cdot \mathbf{\nu }=\overline{u_{\mu }}& & \text{on }\partial {\mu }_u\end{array},$$

$$\begin{array}{rlrl}& \begin{array}{rl}& {\alpha }_{\gamma }^2{\mathbf{u}}_{\gamma }+\nabla p_{\gamma }=0\\ & {\nabla }_{\gamma }\cdot \mathbf{u}+q_{\gamma }=0\end{array}& & \text{in }\gamma \\ & \mathrm{tr}{p}_{\gamma }=\overline{p_{\gamma }}& & \text{on }\partial {\gamma }_p\\ & \mathrm{tr}{\mathbf{u}}_{\gamma }\cdot \mathbf{\nu }=\overline{u_{\gamma }}& & \text{on }\partial {\gamma }_u\end{array},$$

$${\nabla }_{\Omega }\cdot \mathbf{u}:=\nabla \cdot {\mathbf{u}}_{\Omega },$$

$${\nabla }_{\mu }\cdot \mathbf{u}:=\nabla \cdot {\mathbf{u}}_{\mu }-\mathrm{tr}{\mathbf{u}}_{\Omega }\cdot \mathbf{n}+u_{\Gamma },$$

$${\nabla }_{\gamma }\cdot \mathbf{u}:=\nabla \cdot {\mathbf{u}}_{\gamma }-u_{\Gamma },$$

$$\begin{array}{rlrl}& {\alpha }_M^2\mathrm{tr}{\mathbf{u}}_{\Omega }\cdot \mathbf{n}+p_{\mu }-\mathrm{tr}{p}_{\Omega }=0& & \text{on }M\\ & {\alpha }_{\Gamma }^2{u}_{\Gamma }+p_{\gamma }-p_{\mu }=0& & \text{on }\Gamma \end{array}.$$

$${\sigma }_{\Omega }={\alpha }_{\Omega },{\sigma }_{\mu }^2={\alpha }_M^2{ϵ}_{\mu }\mathbf{n}\otimes \mathbf{n}+\frac{{\alpha }_{\mu }^2}{{ϵ}_{\mu }}(I-\mathbf{n}\otimes \mathbf{n}),{\sigma }_{\gamma }^2={\alpha }_{\Gamma }^2{ϵ}_{\gamma }{\mathbf{n}}_{\mu }\otimes {\mathbf{n}}_{\mu }+\frac{{\alpha }_{\gamma }^2}{{ϵ}_{\gamma }}(I-{\mathbf{n}}_{\mu }\otimes {\mathbf{n}}_{\mu }).$$

$$\begin{array}{rl}V:=\{& \mathbf{v}=({\mathbf{v}}_{\Omega },{\mathbf{v}}_{\mu },{\mathbf{v}}_{\gamma },v_{\Gamma }):\\ & {\alpha }_{\Omega }{\mathbf{v}}_{\Omega }\in L^2(\Omega )\wedge {\alpha }_M\mathrm{tr}{\mathbf{v}}_{\Omega }\cdot \mathbf{n}\in L^2(M)\wedge {\nabla }_{\Omega }\cdot \mathbf{v}\in L^2(\Omega )\wedge \mathrm{tr}{\mathbf{v}}_{\Omega }\cdot \mathbf{\nu }=0\text{ on }\partial {\Omega }_u,\\ & {\alpha }_{\mu }{\mathbf{v}}_{\mu }\in L^2(\mu )\wedge {\nabla }_{\mu }\cdot \mathbf{v}\in L^2(\mu )\wedge \mathrm{tr}{\mathbf{v}}_{\mu }\cdot \mathbf{\nu }=0\text{ on }\partial {\mu }_u,\\ & {\alpha }_{\gamma }{\mathbf{v}}_{\gamma }\in L^2(\gamma )\wedge {\nabla }_{\gamma }\cdot \mathbf{v}\in L^2(\gamma )\wedge \mathrm{tr}{\mathbf{v}}_{\gamma }\cdot \mathbf{\nu }=0\text{ on }\partial {\gamma }_u,\\ & {\alpha }_{\Gamma }{v}_{\Gamma }\in L^2(\Gamma )\},\end{array}$$

$$\begin{array}{rl}{\Vert \mathbf{v}\Vert }_V^2:=& {\Vert {\alpha }_{\Omega }{\mathbf{v}}_{\Omega }\Vert }_{\Omega }^2+{\Vert {\nabla }_{\Omega }\cdot \mathbf{v}\Vert }_{\Omega }^2+{\Vert {\alpha }_M\mathrm{tr}{\mathbf{v}}_{\Omega }\cdot \mathbf{n}\Vert }_M^2+{\Vert {\alpha }_{\mu }{\mathbf{v}}_{\mu }\Vert }_{\mu }^2+{\Vert {\nabla }_{\mu }\cdot \mathbf{v}\Vert }_{\mu }^2+\\ & {\Vert {\alpha }_{\gamma }{\mathbf{v}}_{\gamma }\Vert }_{\gamma }^2+{\Vert {\nabla }_{\gamma }\cdot \mathbf{v}\Vert }_{\gamma }^2+{\Vert {\alpha }_{\Gamma }{v}_{\Gamma }\Vert }_{\Gamma }^2.\end{array}$$

$$Q:=\{v=(v_{\Omega },v_{\mu },v_{\gamma }):v_{\Omega }\in L^2(\Omega ),v_{\mu }\in L^2(\mu ),v_{\gamma }\in L^2(\gamma )\},$$

$${\Vert v\Vert }_Q^2:={\Vert v_{\Omega }\Vert }_{\Omega }^2+{\Vert v_{\mu }\Vert }_{\mu }^2+{\Vert v_{\gamma }\Vert }_{\gamma }^2.$$

$$⟨\mathrm{tr}{\mathbf{v}}_{\Omega }\cdot \mathbf{\nu },w{⟩}_{\partial {\Omega }_u}=0\text{for all}w\in H_{0,0}^{\frac{1}{2}}(\partial {\Omega }_u)$$

$$⟨\mathrm{tr}{\mathbf{v}}_{\mu }\cdot \mathbf{\nu },w{⟩}_{\partial {\mu }_u}=0\text{for all}w\in H_{0,0}^{\frac{1}{2}}(\partial {\mu }_u).$$

$$⟨\mathrm{tr}{\mathbf{v}}_{\gamma }\cdot \mathbf{\nu },w{⟩}_{\partial {\gamma }_u}=0\text{for all}w\in H_{0,0}^{\frac{1}{2}}(\partial {\gamma }_u).$$

$$\alpha (\mathbf{u},\mathbf{v}):=$$

$$({\alpha }_{\gamma }{\mathbf{u}}_{\gamma },{\alpha }_{\gamma }{\mathbf{v}}_{\gamma }{)}_{\gamma }+({\alpha }_{\Gamma }{u}_{\Gamma },{\alpha }_{\Gamma }{v}_{\Gamma }{)}_{\Gamma }.$$

$$\begin{array}{rl}\beta (p,\mathbf{v}):=& -(p_{\Omega },{\nabla }_{\Omega }\cdot \mathbf{v}{)}_{\Omega }-(p_{\mu },{\nabla }_{\mu }\cdot \mathbf{v}{)}_{\mu }-(p_{\gamma },{\nabla }_{\gamma }\cdot \mathbf{v}{)}_{\gamma }\\ =& -(p_{\Omega },\nabla \cdot {\mathbf{v}}_{\Omega }{)}_{\Omega }+(p_{\mu },\mathrm{tr}{\mathbf{v}}_{\Omega }\cdot \mathbf{n}{)}_M-(p_{\mu },\nabla \cdot {\mathbf{v}}_{\mu }{)}_{\mu }-(p_{\gamma },\nabla \cdot {\mathbf{v}}_{\gamma }{)}_{\gamma }-\\ & (p_{\mu },v_{\Gamma }{)}_{\Gamma }+(p_{\gamma },v_{\Gamma }{)}_{\Gamma }.\end{array}$$

$$G(\mathbf{v}):=-⟨\overline{p_{\Omega }},\mathrm{tr}{\mathbf{v}}_{\Omega }\cdot \mathbf{\nu }{⟩}_{\partial {\Omega }_p}-⟨\overline{p_{\mu }},\mathrm{tr}{\mathbf{v}}_{\mu }\cdot \mathbf{\nu }{⟩}_{\partial {\mu }_p}-⟨\overline{p_{\gamma }},\mathrm{tr}{\mathbf{v}}_{\gamma }\cdot \mathbf{\nu }{⟩}_{\partial {\gamma }_p}$$

$$F(v):=(q_{\Omega },v_{\Omega }{)}_{\Omega }+(q_{\mu },v_{\mu }{)}_{\mu }+(q_{\gamma },v_{\gamma }{)}_{\gamma }.$$

$$\begin{array}{rlrl}& \alpha (\mathbf{u},\mathbf{v})+\beta (p,\mathbf{v})=G(\mathbf{v})& & \forall \mathbf{v}\in V\\ & \beta (v,\mathbf{u})=F(v)& & \forall v\in Q\end{array}.$$

$$\left|\alpha (\mathbf{u},\mathbf{v})\right|\le$$

$${\Vert {\alpha }_{\gamma }{\mathbf{u}}_{\gamma }\Vert }_{\gamma }{\Vert {\alpha }_{\gamma }{\mathbf{v}}_{\gamma }\Vert }_{\gamma }+{\Vert {\alpha }_{\Gamma }{u}_{\Gamma }\Vert }_{\Gamma }{\Vert {\alpha }_{\Gamma }{v}_{\Gamma }\Vert }_{\Gamma }\le {\Vert \mathbf{u}\Vert }_V{\Vert \mathbf{v}\Vert }_V$$

$$\left|\beta (p,\mathbf{v})\right|\le$$

$$\left|G(\mathbf{v})\right|\le$$

$${\Vert \overline{p_{\gamma }}\Vert }_{H^{\frac{1}{2}}(\partial {\gamma }_p)}{\Vert \mathrm{tr}{\mathbf{v}}_{\gamma }\cdot \mathbf{\nu }\Vert }_{H^{-\frac{1}{2}}(\partial {\gamma }_p)}\le c_G{\Vert \mathbf{v}\Vert }_V$$

$$\left|F(v)\right|\le$$

$$c_G=\max \left\{{\Vert \overline{p_{\Omega }}\Vert }_{H^{\frac{1}{2}}(\partial {\Omega }_p)},{\Vert \overline{p_{\mu }}\Vert }_{H^{\frac{1}{2}}(\partial {\mu }_p)},{\Vert \overline{p_{\gamma }}\Vert }_{H^{\frac{1}{2}}(\partial {\gamma }_p)}\right\},$$

$$c_F=\max \left\{{\Vert q_{\Omega }\Vert }_{\Omega },{\Vert q_{\mu }\Vert }_{\mu },{\Vert q_{\gamma }\Vert }_{\gamma }\right\}.$$

$$v=(v_{\Omega },0,0)$$

$$v=(0,v_{\mu },0)$$

$$v=(0,0,v_{\gamma })$$

$${\Vert \mathbf{w}\Vert }_V^2={\Vert {\alpha }_{\Omega }{\mathbf{v}}_{\Omega }\Vert }_{\Omega }^2+{\Vert {\alpha }_M\mathrm{tr}{\mathbf{v}}_{\Omega }\cdot \mathbf{n}\Vert }_M^2+{\Vert {\alpha }_{\mu }{\mathbf{v}}_{\mu }\Vert }_{\mu }^2+{\Vert {\alpha }_{\gamma }{\mathbf{v}}_{\gamma }\Vert }_{\gamma }^2+{\Vert {\alpha }_{\Gamma }{v}_{\Gamma }\Vert }_{\Gamma }^2.$$

$$\begin{array}{rlrl}& -\nabla \cdot \nabla {\phi }_{\Omega }=v_{\Omega }& & \text{in }\Omega \\ & -\mathrm{tr}\nabla {\phi }_{\Omega }\cdot \mathbf{n}=v_{\mu }& & \text{on }M\\ & \mathrm{tr}{\phi }_{\Omega }=0& & \text{on }M\\ & \mathrm{tr}{\phi }_{\Omega }=0& & \text{on }\partial {\Omega }_p\end{array}\begin{array}{rlrl}& -\nabla \cdot \nabla {\phi }_{\mu }=v_{\mu }+v_{\gamma }& & \text{in }\mu \\ & {\phi }_{\mu }=0& & \text{on }\partial {\mu }_p\\ & -\nabla \cdot \nabla {\phi }_{\gamma }=-v_{\mu }& & \text{in }\gamma \\ & {\phi }_{\gamma }=0& & \text{on }\partial {\gamma }_p\end{array}$$

$${\phi }_{\Omega }\in H_{\ast }^2(\Omega )\text{with}{\Vert {\phi }_{\Omega }\Vert }_{H_{\ast }^2(\Omega )}\le {\Vert v_{\Omega }\Vert }_{\Omega }+{\Vert v_{\mu }\Vert }_{\mu },$$

Full text

1 Introduction

The accurate description and simulation of fluid flow in geological porous media are of paramount importance for several industrial applications, such as: CO2 injection and sequestration, geothermal exploitation, oil migration and recovery, and prevention of groundwater contamination due to nuclear waste disposal, just to name a few. See [5, 28, 24]. Despite the relevance of the topic, a number of challenges are not yet fully solved, in particular related to the presence of faults and fractures in the domain of interest. Faults and fractures are the portion of the porous media where the rock has been broken due to geological movements of the upper crust. In this study we consider only large faults.

A tectonic fault is the result of a relative displacement of two parts of the upper crust happened over geological eras. This displacement is accommodate by a high strain region, the fault core, surrounded on both sides by a highly fractured region, the damage zone. These layer contain several fractures, on a much smaller scale than the fault, which may alter the local properties of the flow path. Faults have a thickness which is several order of magnitude smaller than any other characteristic sizes in the porous domain, however their physical properties may greatly differ from the porous media. Due to infilling processes and chemical reactions these objects may be partially or completely occluded, thus they can behave as preferential conduits or geological barriers for the fluid flow. The surrounding damage zone may or may not had experienced the same processes behaving similarly to, or differently from than the related fault. The development of accurate conceptual models is a key factor to be able to include these objects and their effect in a simulation code and obtain reliable outcomes and predictions. The aim of this work is to devise a new effective conceptual model to account for multiple thin regions in porous media.

One possibility to account for the presence of faults is to characterize this region with a different permeability, but its small thickness makes difficult or even unrealistic its inclusion in the grid representing the rock matrix. A common approach, introduced in [25, 40, 16, 20], is to consider faults as a lower dimensional objects and derive a new conceptual model to describe the flow and pressure behavior inside and across these objects. This approach has been successfully applied to different kind of physics, ranging from advection of a passive scalar, heat transport, multi-phase flow. An incomplete list of references is the following [31, 35, 23, 3, 32, 26, 12, 19, 18, 33, 22, 37, 30, 38, 13, 1, 10, 27, 8, 4, 39, 17, 7, 9, 34].

In this work we extend the previously introduced models to consider also the damage zone as a lower dimensional objects which are connected on one side to the rock matrix and, on the other side, to the fault. The aim is to be able to simulate different scenarios where the rock matrix, damage zone, and fault may have different permeability values without resorting to extreme mesh refinement to capture the thickness of the damage and core zones. Moreover this multilayer approach allows for different apertures and asymmetries across the fault, unlike the previous models.

The numerical discretization is based on the classical Raviart-Thomas-Nédélec approximation for the flux field and a constant piecewise representation of the pressure. The resulting scheme is locally mass conservative and is thus suitable for coupling with a transport problem. For clarity in the exposition, we consider only one single fault, since the case of intersecting faults requires additional model complexities.

For the implementation of the numerical examples, we have used the library PorePy [29], which is a simulation tool written in Python for fractured and deformable porous media. The numerical tests presented in this paper are available in the GitHub repository of the library. The main contribution of this work is the introduction of the multi-layers interface law, valid for any dimension. Even if we present an approximation and analysis based on the lowest order Raviart-Thomas-Nédélec, the implementation is agnostic with respect to the numerical scheme. It is indeed possible to use any other scheme present in the library, like two and multi point flux approximation or the mixed virtual element method.

The paper is organised as follow: in Section 2 both the equi and mixed-dimensional mathematical models are introduced and discussed. Section 3 deals with the weak formulation of the mixed-dimensional problem: functional spaces, weak problem, and it well posedness. In Section 4 we briefly describe the numerical scheme and how to handle the mixed-dimensional nature of the problem. Section 5 contains two numerical test to validate the proposed model. Finally, Section 6 is devoted to conclusions.

2 Mathematical problem

In this section we present the mathematical models considered to describe the pressure and Darcy velocity governed by the single-phase Darcy flow equations. In Subsection 2.1, we recall the classical formulation where the fault and the damage zone are considered equi-dimensional with respect to the rock matrix, i.e. are represented as two dimensional in a two dimensional porous matrix, and three dimensional in 3D. Their thickness is thus explicitly represented in the domain and captured by the computational grid. By considering a model reduction approach, Subsection 2.2 presents the new formulation where the fault and damage zone are now represented as lower-dimensional objects and new equations are introduced.

2.1 Equi-dimensional model

Let $\Omega$ denote the rock host matrix, $\mu$ denote the damage zone and $\gamma$ the fault core. The interfaces $M$ and $\Gamma$ between these objects can be defined as $\overline{M}=\overline{\partial\Omega}\cap\overline{\mu}$ and $\overline{\Gamma}=\overline{\mu}\cap\overline{\gamma}$, respectively. We define $\partial\Omega_{\rm int}$ as the part of the boundary of $\Omega$ which is in contact with $M$, and with $\partial\Omega$ the part of the boundary of $\Omega$ which is not in contact with $M$. We further subdivide $\partial\Omega$ into two parts, $\partial\Omega_{p}$ and $\partial\Omega_{u}$, such that $\overline{\partial\Omega}=\overline{\partial\Omega_{p}}\cup\overline{\partial\Omega_{u}}$ with $\partial\Omega_{p}\neq\emptyset$ and $\emptyset=\mathring{\partial\Omega_{p}}\cup\mathring{\partial\Omega_{u}}$. In $\partial\Omega_{p}$ we will impose natural boundary conditions, while in $\partial\Omega_{u}$ essential boundary conditions. In the same way, we define the boundary of the damage zone $\mu$ as $\partial\mu_{\rm int}$ and $\partial\mu$, being the internal and external portion of the boundary of $\mu$, respectively. The external part $\partial\mu$ can be divided into two parts $\partial\mu_{p}$ and $\partial\mu_{u}$ such that $\overline{\partial\mu}=\overline{\partial\mu_{p}}\cup\overline{\partial\mu_{u}}$ with $\partial\mu_{p}\neq\emptyset$ and $\emptyset=\mathring{\partial\mu_{p}}\cup\mathring{\partial\mu_{u}}$. In $\partial\mu_{p}$ we will impose natural boundary conditions, while in $\partial\mu_{u}$ essential boundary conditions. The same nomenclature is introduced for the fault $\gamma$. We define the unit normal $\bm{n}$, associated with $\partial\Omega_{\rm int}$, which points from $\Omega$ to $M$. Similarly, we introduce the unit normal $\bm{n}_{\mu}$, associated with $\partial\mu_{\rm int}$, which points from $\mu$ to $\Gamma$. Finally, $\bm{\nu}$ is the external outward unit normal of the domain. See Figure 1 as an example.

We assume that the rock matrix $\Omega$, damage zone $\mu$, and fault $\gamma$ have the same spatial dimension. The damage zone and fault are characterized by one of the dimension to be much (orders of magnitude) smaller than the others.

We are interested in the mathematical description of the Darcy velocity and pore pressure described by the Darcy problem. We indicate with $p_{\Upsilon}$ the pressure and with $\bm{u}_{\Upsilon}$ the Darcy velocity, on the portion (rock matrix, damage zone, or fault) $\Upsilon$ of the problem, being $\Upsilon$ equal to $\Omega$, $\mu$, and $\gamma$. The system of equations reads: find $(\bm{u}_{\Upsilon},p_{\Upsilon})$ such that

[TABLE]

Even if not explicitly used in the previous set of equations, we indicate with $\epsilon_{\mu}$ and $\epsilon_{\gamma}$ the thickness of the damage zone and fault, respectively. To summarize the equi-dimensional system of equations, we introduce the following problem.

Problem 1** (Equi-dimensional model problem).**

Find $(\bm{u}_{\Upsilon},p_{\Upsilon})$ such that (1) is satisfied, for $\Upsilon$ equal to $\Omega$, $\mu$, and $\gamma$.

2.2 Mixed-dimensional model

The geometrical reduction of the model approximates the thin regions, in our case the damage zone and the fault, by their center line (a lower dimensional object) and derives new equations and coupling conditions to describe the Darcy velocity and pressure field in the new setting. See Figure 2 as an example. In this work we follow the reduction procedure described in the literature in [2, 15, 31]. To keep the notation simple, we preserve the same notation for the rock matrix $\Omega$, the damage zone $\mu$ and fault $\gamma$ even if the domains are geometrically different from the equi-dimensional case: in particular $\Omega$ is extended up to the center line of the fault core, while the fault and the damage zone shrink and become overlapped lower dimensional interfaces.

The mixed-dimensional problem describes the Darcy velocity and pressure field defined on the tangent space of each domain. At the interfaces, to couple the problems, we consider traces of the variables and the additional unknown $u_{\Gamma}$, which can be interpreted as a normal Darcy velocity from the damage zone $\mu$ to the fault $\gamma$. To simplify the notation, we introduce the pressure and Darcy velocity compounds $p$ and $\bm{u}$, respectively, as

[TABLE]

Note that $\bm{u}_{\mu}$ and $\bm{u}_{\gamma}$ are fluxes in the lower dimensional layers and $u_{\Gamma}$ is a scalar variable representing the normal flux exchanged by lower dimensional objects. Following the approaches presented in [31, 40, 16, 20], we consider the Darcy model in each domain separately. First we consider the rock matrix,

[TABLE]

In the previous equations the data are: $\alpha_{M}^{2}$ which is the inverse of the effective normal matrix permeability between the rock matrix and the damage zone and $\alpha_{\Gamma}^{2}$ which is the inverse of the effective normal matrix permeability between the damage zone and the fault.

Problem 2** (Mixed-dimensional model problem).**

Find $(\bm{u},p)$ such that (2) is satisfied.

Remark 1**.**

In the mixed-dimensional Problem 2 we have used directly the effective permeabilities, i.e. the permeabilities already scaled by the thickness of the related object. They can be related with the parameters of the equi-dimensional Problem 1 as follows

[TABLE]

3 Weak problem

In this section we introduce the functional setting and the weak formulation of Problem 2. For simplicity, we consider null flux essential boundary conditions which, through a lifting technique, can be generalized. In the forthcoming parts we indicate with ${\left\|{\cdot}\right\|}_{A}$ the $L^{2}$-norm on the set $A$ and with $(\cdot,\cdot)_{A}$ the $L^{2}$-scalar product on the set A.

3.1 Functional space setting

Since the previous problem couples unknowns in all the domains, it is natural to introduce a global functional space to respect this structure. The key motivation is related to the mixed-dimensional divergence operators, (2d) that are related to a space that generalise the classical (weighted) $H_{\nabla\cdot}-$space in this framework. For the vector fields, i.e. $\bm{u}=(\bm{u}_{\Omega},\bm{u}_{\mu},\bm{u}_{\gamma},u_{\Gamma})$, we introduce the following

[TABLE]

with the associated norm that make the space $V$ complete, defined as

[TABLE]

We assume that exists $c_{\Omega}>0$ such that $c_{\Omega}\leq\alpha_{\Omega}$ a.e., and $\alpha_{\Omega}\in L^{\infty}(\Omega)$. Similarly, we assume that there exists $c_{M}>0$ such that $c_{M}\leq\alpha_{M}$ a.e., and $\alpha_{M}\in L^{\infty}(M)$. For the damage zone, we require that exists $c_{\mu}>0$ such that $c_{\mu}\leq\alpha_{\mu}$ a.e., and also $\alpha_{\mu}\in L^{\infty}(\mu)$. Again, we assume that exists $c_{\gamma}>0$ such that $c_{\gamma}\leq\alpha_{\gamma}$ a.e., and also $\alpha_{\gamma}\in L^{\infty}(\gamma)$. Finally, we require that exists $c_{\Gamma}>0$ such that $c_{\Gamma}\leq\alpha_{\Gamma}$ a.e., and also $\alpha_{\Gamma}\in L^{\infty}(\Gamma)$. The extra regularity required for the trace of $\bm{v}_{\Omega}$ on $M$ is related to the interface conditions on $M$ which can be seen as a Robin-type boundary conditions for a problem in mixed form. For more details see [36, 31].

The functional space for the pressure is defined again on the compound, namely

[TABLE]

with associated norm that make the space $Q$ complete, given by

[TABLE]

Remark 2** (Boundary conditions).**

The conditions related to the boundaries in the definition of (3) have to be interpreted in a proper way. The first condition $\operatorname*{\operatorname{tr}}\bm{v}_{\Omega}\cdot\bm{\nu}=0$ on $\partial\Omega_{u}$ corresponds to

[TABLE]

with $\langle\cdot,\cdot\rangle_{A}$ the duality pairing defined as $\langle\cdot,\cdot\rangle_{A}:H^{-\frac{1}{2}}(A)\times H^{\frac{1}{2}}_{0,0}(A)\rightarrow\mathbb{R}$, with $A$ a generic set. From standard results we have $\operatorname*{\operatorname{tr}}\bm{v}_{\Omega}\cdot\bm{\nu}\in H^{-\frac{1}{2}}(\partial\Omega_{u})$. In a similar way, the condition $\operatorname*{\operatorname{tr}}\bm{v}_{\mu}\cdot\bm{\nu}=0$ on $\partial\mu_{u}$ means

[TABLE]

In this case, we have $\operatorname*{\operatorname{tr}}\bm{v}_{\mu}\cdot\bm{\nu}\in H^{-\frac{1}{2}}(\partial\mu_{u})$. Finally, the condition $\operatorname*{\operatorname{tr}}\bm{v}_{\gamma}\cdot\bm{\nu}=0$ on $\partial\gamma_{u}$ is interpreted similarly as

[TABLE]

In this case we have $\operatorname*{\operatorname{tr}}\bm{v}_{\gamma}\cdot\bm{\nu}\in H^{-\frac{1}{2}}(\partial\gamma_{u})$.

Remark 3**.**

A more rigorous approach requires the introduction of interpolation operators between $\partial\Omega_{\rm int}$ and $M$, as well as between $\mu$ and $\Gamma$ and also between $\gamma$ and $\Gamma$. However, to simplify the presentation we consider them implicit.

3.2 Weak formulation

We propose now the weak formulation of the previous problem. Again, we respect the mixed-dimensional nature of Problem 2 by introducing bilinear forms that consider variables defined on different domains. We first introduce the weighted $H_{\nabla\cdot}-$mass bilinear form, indicated by $\alpha$ and defined as $\alpha:V\times V\rightarrow\mathbb{R}$ such that

[TABLE]

The bilinear form associated to the conservation statement is given by: $\beta:Q\times V\rightarrow\mathbb{R}$ such that

[TABLE]

The functionals, which contain the boundary data and source terms, are defined as: $G:V\rightarrow\mathbb{R}$ and $F:Q\rightarrow\mathbb{R}$ such that

[TABLE]

Where we assume that the pressure boundary data are such that $\overline{p_{\Omega}}\in H^{\frac{1}{2}}(\partial\Omega_{p})$, $\overline{p_{\mu}}\in H^{\frac{1}{2}}(\partial\mu_{p})$, and $\overline{p_{\gamma}}\in H^{\frac{1}{2}}(\partial\gamma_{p})$. We require also that the source terms belong to $(q_{\Omega},q_{\mu},q_{\gamma})\in Q$.

Problem 3** (Weak formulation).**

The weak formulation of Problem 2 reads: find $(\bm{u},p)\in V\times Q$ such that

[TABLE]

3.3 Well posedness

In the proof of the well posedness of the problem, some parts are inspired by [18, 21].

Theorem 1**.**

Problem 3 is well posed.

Proof.

The bilinear forms and functionals in Problem 3 are linear, so we first prove their continuity. For simplicity, we assume that the portions of the boundary $\partial\Omega_{u}$, $\partial\mu_{u}$, and $\partial\gamma_{u}$ are empty. A lifting technique can be used in the general case. By using Cauchy-Schwarz and triangular inequalities we get

[TABLE]

as well as for the $\beta$ bilinear form

[TABLE]

For the functionals, we consider in addition the trace inequality associated to the $H_{\nabla\cdot}-$spaces. We obtain

[TABLE]

with the constants

[TABLE]

Now, we move on to proving the coercivity of $\alpha$ on the kernel of $\beta$. Considering a function $\bm{w}\in V$ such that $\beta(v,\bm{w})=0$ for all $v\in Q$, we have by a particular choice of $v$ the following results

[TABLE]

thus the norm of $\bm{w}$ simplifies to

[TABLE]

We can prove the coercivity of $\alpha$ on the kernel of $\beta$, by simply notice that $\alpha(\bm{w},\bm{w})={\left\|{\bm{w}}\right\|}_{V}^{2}$.

To prove the inf-sup condition, given a function $v\in Q$, we introduce the following auxiliary problems

[TABLE]

while for the intersection $\Gamma$ we have $\varphi_{\Gamma}=v_{\mu}+v_{\gamma}$. By assuming that the domains are regular enough, from [14] an elliptic regularity result can be used giving

[TABLE]

The space $H^{2}_{*}(\Omega)$ is the broken $H^{2}-$space defined on each connected component of $\Omega$ (e.g. the left and right part in Figure 2), its norm is defined coherently. We clearly have ${\left\|{\varphi_{\Gamma}}\right\|}_{\Gamma}\leq{\left\|{v_{\mu}}\right\|}_{\mu}+{\left\|{v_{\gamma}}\right\|}_{\gamma}$. By considering the function $\bm{w}\in V$ such that

[TABLE]

we obtain that $-\nabla\cdot\bm{w}_{\Omega}=v_{\Omega}$ and $\operatorname*{\operatorname{tr}}\bm{w}_{\Omega}\cdot\bm{n}=v_{\mu}$. For $\bm{w}_{\mu}$ and $\bm{w}_{\gamma}$ we get that $-\nabla\cdot\bm{w}_{\mu}=v_{\mu}+v_{\gamma}$ and $-\nabla\cdot\bm{w}_{\gamma}=-v_{\mu}$, respectively. Finally, we have $w_{\Gamma}=v_{\gamma}+v_{\mu}$. This gives the following expressions for the mixed-dimensional divergences $\nabla_{\Omega}\cdot\bm{w}=-v_{\Omega}$, $\nabla_{\mu}\cdot\bm{w}=-v_{\mu}$, and $\nabla_{\gamma}\cdot\bm{w}=-v_{\gamma}$. The $V-$norm of $\bm{w}$ can be bound by the $Q-$norm of $v$, in fact

[TABLE]

Finally we obtain the boundedness of $\beta$ from below with this choice of $\bm{w}$, namely

[TABLE]

Thus the inf-sup condition is fulfilled, and following [11] we conclude that Problem 3 is well posed. ∎

4 Numerical discretization

To keep the mixed nature of Problem 3 in the numerical approximation, we consider the lowest order Raviart-Thomas-Nédélec $\mathbb{RT}_{0}$ finite element for the Darcy velocity and piecewise constant $\mathbb{P}_{0}$ for the pressure in the domains $\Omega$, $\mu$, and $\gamma$. The choice of the pair $\mathbb{RT}_{0}-\mathbb{P}_{0}$ is also motivated by their local mass conservation property and consistency with the functional spaces considered in the weak formulation.

We introduce a family of simplicial meshes approximations of the domains $\Omega$, $\mu$, $\gamma$, $\Gamma$ respectively, which will be indicated with the same symbol. By assuming a matching discretization of the meshes at the interfaces, the approximation of the functional spaces are defined as

[TABLE]

we introduce a set of base functions for the discrete spaces, such that

[TABLE]

We assume that quadrature is performed exactly. We indicate with $h$ the global mesh size of the discretization.

4.1 Matrix formulation

Given the choice of the discrete spaces, we can recast the weak formulation of the problem in term of a block matrix system. We introduce the block matrices related to the mass matrices

[TABLE]

the matrices associated with the conservation statement in each domain are defined as

[TABLE]

while between the domains the coupling conditions are associated with the matrices

[TABLE]

We indicate with $\mathfrak{u}$ and $\mathfrak{p}$ the values of degrees of freedom associated with the Darcy velocity and pressure. We make use of a similar notation to indicate the source and boundary terms. The linear system associated with Problem 3 is

[TABLE]

where we have avoided explicitly writing empty matrices. Since the bilinear forms associated to the matrices $A_{\Omega}$, $A_{\mu}$, $A_{\gamma}$, $A_{\Gamma}$ are symmetric then also these matrices are symmetric, resulting in a symmetric global system with a saddle-point structure. We can introduce the following problem.

Problem 4** (Matrix formulation).**

The matrix formulation of Problem 3 is: find $(\mathfrak{u},\mathfrak{p})$ such that (4) is satisfied.

To recast the previous system in terms of the pressures alone, we introduce the matrices

[TABLE]

where the matrices $A_{\Omega}$, $A_{\mu}$, $A_{\gamma}$, and $A_{\Gamma}$ are invertible since they are a Raviart-Thomas-Nédélec approximation of the $H_{\nabla\cdot}-$mass bilinear forms and $\partial\Omega_{p}\neq\emptyset$, $\partial\mu_{p}\neq\emptyset$, and $\partial\gamma_{p}\neq\emptyset$, while $A_{\Gamma}$ is invertible by construction. The previous matrices are thus well defined and $S_{\Omega}$, $S_{\mu}$, $S_{\gamma}$ are also symmetric and positive definite. We introduce the vectors

[TABLE]

which are again well defined. The system in terms of pressure can be written as

[TABLE]

Problem 5** (Pressure matrix formulation).**

The pressure matrix formulation of Problem 3 is: find $\mathfrak{p}$ such that (5) is satisfied.

To numerically solve the problem, it is possible to use the equivalent formulations of Problem 4 or Problem 5.

5 Numerical examples

In this part we present some numerical examples to show different aspects of the mathematical model introduced previously. In particular, we focus on its behavior in the presence of high contrasts in permeability among the rock matrix, damage zone and fault. We present also the model error, i.e. the error introduced by the geometrical reduction of the three layers and discuss the obtained results. Finally, we show the applicability in a three-dimensional domain.

All the test cases are implemented in the library PorePy [29]. The current implementation in the code considers a mortar variable on each interface and it is thus capable to handle non-matching discretization. Also, the code is agnostic with respect to the numerical scheme adopted in each domain. However, to keep the presentation simple and coherent with the previous sections, we will consider matching grids and Raviart-Thomas-Nédélec finite element of the lowest order for the numerical approximation.

5.1 Example 1

This first set of numerical tests is divided in two parts. In the first we present the effect of permeability heterogeneity in the fault and damage zone on the solution. This test is inspired by the examples presented in [31, 19]. We study the case of high permeability, low permeability, and a mixed case. The aim is to present the potentialities of the introduced model in the Sub-subsection 5.1.1. In the second part, Sub-subsection 5.1.2, the model error is studied comparing the mixed-dimensional solution with the equi-dimensional one where a full Darcy problem is solved on a computational grid refined enough to resolve the aperture of both the damage zone and the fault.

In all the cases we consider a fixed geometry and boundary conditions on the rock matrix. The rock matrix occupies the domain $\Omega=(0,1)\times(0,1)\cup(1,2)\times(0,1)$, while the damage zone and fault are identified as $\mu=\{1\}\times(0,1)$ and $\gamma=\{1\}\times(0,1)$, respectively. We set $\alpha_{\Omega}=1$ and pressure boundary condition on the left and right of $\Omega$, with values 0 and 1, respectively. On the remaining portions of the boundary we impose zero flux. The computational domain is represented in Figure 2. The grid is composed of $\sim 8.5k$ triangles for $\Omega$, $80$ segments for $\mu$, $40$ segments for $\gamma$, and $40$ segments for $\Gamma$.

We consider three different sub-cases, depending on the value of $\alpha$ in the damage zone and fault. In all the cases, we set $\epsilon=10^{-2}$. In the case (i) we have $\alpha_{\mu}^{2}=10^{2}\epsilon$ and $\alpha_{M}^{2}=10^{2}/\epsilon$ for the damage zone, and $\alpha_{\gamma}^{2}=10^{2}\epsilon$ and $\alpha_{\Gamma}^{2}=10^{2}/\epsilon$ for the fault. For the case (ii), the parameters are chosen as $\alpha_{\mu}^{2}=k\epsilon$ and $\alpha_{M}^{2}=k/\epsilon$ for the damage zone, and $\alpha_{\gamma}^{2}=k\epsilon$ and $\alpha_{\Gamma}^{2}=k/\epsilon$ for the fault, where $k$ is given by

[TABLE]

Finally, in case (iii) we impose the values of $\alpha_{\gamma}^{2}=k\epsilon$ and $\alpha_{\Gamma}^{2}=k/\epsilon$ for the fault. The damage zone is divided into its left and right parts, on the former we have $\alpha_{\mu}^{2}=k\epsilon$ and $\alpha_{M}^{2}=k/\epsilon$ while on the latter $\alpha_{\mu}^{2}=10^{2}\epsilon$ and $\alpha_{M}^{2}=10^{2}/\epsilon$. In case (i) the lower dimensional objects are highly conductive, in case (ii) are heterogeneous in space and less conductive in the central part, and in case (iii) the damage zone is heterogeneous in space and asymmetric.

In case (i) we impose unitary pressure at damage zone and fault tips $y=1$, while zero pressure on the other side. For case (ii) and case (iii) zero flux is imposed on the boundary of the damage zone and fault.

5.1.1 Permeability contrast

In this first section, we present the numerical results obtained from the reduced model in the three different cases. In Figure 3 the results from case (i) are shown.

Due to the high value of $\alpha_{\mu}$, $\alpha_{M}$, $\alpha_{\gamma}$, and $\alpha_{\Gamma}$ the pressure profile is smooth (continuous) across $\mu$ and $\gamma$. The impact of these lower dimensional objects on the flow is still quite remarkable, due to the type of boundary conditions and the permeability contrast.

The solution for the case (ii) is given in Figure 4.

The impact of the low value of the permeabilities in the central parts of $\mu$ and $\gamma$ is evident. We have a pressure jump between the rock matrix and the damage zone, and between the damage zone and the fault. The pressure in the latter is constant due to the symmetry of the problem. The flow tends to focus around these less permeable regions.

Figure 5 shows the results from the case (iii).

In this test a side of $\mu$ has a portion with low $\alpha_{\mu}$ and $\alpha_{M}$, while on the other side these coefficients have hight values. The pressure solution exhibit thus a jump between the rock matrix and the first side of $\mu$, and with the latter and the fault. However, on the other side the solution behaves similarly to case (i) and we obtain a smooth (continuous) profile among the fault and the second part of $\mu$ and the fault. The Darcy velocity tends to avoid the low permeable part of $\mu$, but not the high permeable one.

The test cases above demonstrate the capability of the model to handle different combinations of the model parameters. The obtained solutions are physically sound and show the capability and potentiality of the model.

5.1.2 Model error

In this section we discuss the error associated with the geometrical reduction of the fault and damage zone, which are modeled as $n-1$ dimensional interfaces even if, physically, they are $n$ dimensional regions as $\Omega$. To estimate this error we will compare the solutions provided by the mixed-dimensional model with the numerical solution of a traditional full Darcy problem set on a domain with heterogeneous permeabilities, discretized with a grid that is able to resolve the actual aperture of the layers. The equi-dimensional solution is computed with the same pair of mixed finite elements as the mixed-dimensional one. Let $p_{\mathrm{equi}}^{\eta}$ be the numerical solution of the equi-dimensional problem on a grid $\mathcal{T}_{\eta}$ of size $\eta$. Let $I^{\eta}(p^{h}_{\Omega})$ be the interpolation of $p^{h}_{\Omega}$ on $\mathcal{T}_{\eta}$: in Figure 6 we show the difference $I^{\eta}(p^{h}_{\Omega})-p_{\mathrm{equi}}^{\eta}$ for the three cases presented in the previous section.

We can observe that, as expected, for the second and third cases the difference is mostly focused in the fault and surrounding layers: there indeed the equi-dimensional solution exhibits strong gradients which are replaced by jumps in the mixed-dimensional approximation due to the shrinking of the layers into interfaces. It can also be observed the asymmetry in the model error of the third case, reflecting the different permeabilities of the two layers of the damage zone. In the first case however, where the fault is conductive, pressure is continuous everywhere and the largest values of the error are due to the strong gradients close to the fault which are not well captured by the coarse mixed-dimensional grid.

For a more quantitative analysis we compute $\tilde{e}_{m}$ as the $L^{2}(\Omega)$ norm of $I^{\eta}(p^{h}_{\Omega})-p_{\mathrm{equi}}^{\eta}$ for different apertures $\epsilon$ of the fault and surrounding layers. If the geometrical reduction is consistent we expect a reduction of this error with smaller values of $\epsilon$. The results for the three cases are reported in Table 1. Note that as $\epsilon$ decreases we observe a saturation of the error. This is due to the fact that we are indeed measuring together the model error and the numerical error, i.e. the error due to a coarse mixed-dimensional grid.

To isolate the two effects we proceed as follows. Assuming that $p_{\mathrm{equi}}^{\eta}$ is fully resolved and can replace the exact equi-dimensional solution, we define the model error as

[TABLE]

where $p_{\Omega}^{\mathrm{exact}}$ is the fully resolved mixed-dimensional solution, which is in general not available since the aim of reduced model is to avoid extreme refinement. Let $p_{\Omega}^{h_{2}}$ be the solution for a second mixed-dimensional grid with $h_{2}<h$, then

[TABLE]

Here $\tilde{e}_{m}$ is the previous - incorrect - estimate of the model error, and $\Delta p$ estimates the effect of grid refinement. If we assume that $e_{h_{2}}$ is small, since $h_{2}<h$, then we state that $e_{m}\leq\tilde{e}_{m}+\Delta p$. With similar arguments we can also obtain $e_{m}\geq\tilde{e}_{m}-\Delta p$. The values of this upper and lower bounds are reported for the three cases in Figure 7, where we can observe that the lower bound decreases linearly with $\epsilon$ as expected.

5.2 A three-dimensional example

This last test case is inspired by Case 1: single fracture of the benchmark initiative [6]. We consider a single fault immersed in a three-dimensional rock matrix, defined as $\Omega=(0,100)^{3}$. The fault has thickness equal to $\epsilon_{\gamma}=10^{-3}$ and permeability equal to $\alpha_{\gamma}^{2}=10^{-7}\epsilon_{\gamma}$ and $\alpha_{\Gamma}^{2}=10^{-7}/\epsilon_{\gamma}$. The fault is identified by the following four corners

[TABLE]

In the rock matrix the permeability is $\alpha^{2}_{\Omega}=10^{-6}$ for $z\geq 10$, otherwise $\alpha^{2}_{\Omega}=10^{-5}$. For the boundary condition, we impose 4 for the pressure in the narrow band $\{0\}\times(0,100)\times(90,100)$ and 0 in the portion $(0,100)\times\{0\}\times(0,100)$ of the boundary. The assume zero flux elsewhere. For the damage zone we have a thickness, greater than the fault, equal to $\epsilon_{\mu}=10^{-1}$ with permeability $\alpha_{\mu}^{2}=10^{-2}\epsilon_{\mu}$ and $\alpha_{M}^{2}=10^{-2}/\epsilon_{\mu}$ on the upper part and $\alpha_{\mu}^{2}=10^{-1}\epsilon_{\mu}$ and $\alpha_{M}^{2}=10^{-1}/\epsilon_{\mu}$ on the lower part. The numerical solution considers a grid composed of $\sim 9.5k$ tetrahedra for $\Omega$, $757$ triangles for the fault grid $\gamma$, and $\sim 1.5k$ triangles for the damage zone $\mu$. This correspond to the coarsest mesh of the benchmark [6].

In Figure 8, we present the obtained numerical solution. On the left the pressure field; since the fault is less permeable than the damage zone and rock matrix the solution exhibits a steep variation across the fault. On the right the velocity is represented with different colours for the rock matrix, fault, and each part (top in blue and bottom in red) of the damage zone. The arrows corresponding to the rock matrix and the bottom part of the damage zone are enlarged by a factor $10^{8}$, while for the upper part of the damage zone we use a factor $10^{9}$. The velocity along the fault is zero because of its low permeability. Since the rock matrix is much less permeable than the damage zone the flow tends to be more concentrated in the latter. The arrows are coherently pointing from the inflow part to the outflow

Also in this three-dimensional case, the obtained solution is physically sound and shows the capability and potentiality of the model.

6 Conclusion

In this work we have introduced a new conceptual model for the simulation of Darcy flows in porous media crossed by large faults, i.e. by complex regions characterized by an inner thin core surrounded by damage zone where, due to the accommodation of strain, a large number of small fractures is present. Following a well established literature, we approximate the thin fault region, and in particular both the core and the damage zone as lower-dimensional and geometrically coincident objects, i.e. lines in $2d$ and surfaces in $3d$ to avoid extreme mesh refinement. However, unlike previous works, the presence of three lower dimensional interfaces instead of a single fault object gives us the freedom to better characterize the fluid-dynamic behavior of the structure, by using different permeability and a different thickness for the core and the damage zone, and even accounting for asymmetries across the fault. Moreover, we highlight the fact that this approach can be extended, in different areas of application, to the efficient simulation of thin layered porous media.

We have proven the well posedness of the weak formulation of the mixed dimensional problem and discussed its numerical approximation and a set of examples. The numerical tests confirm the ability of the resulting numerical scheme to handle high contrasts in permeability among the rock matrix, fault, and damage zone. Moreover, by comparing the mixed dimensional solution with fully resolved equi-dimensional simulations we have shown that the model error associated with geometrical model reduction is only focused where the fault or the damage zone give a jump in the pressure field and, moreover, this error decreases for thinner faults as expected. The obtained solutions are thus physically sound and the model can be regarded as a promising strategy. Possible future developments include the possibility to account for heterogeneities, and changes in time, in the aperture and permeability of the different layers. These variability could be the result of geochemical or mechanical processes. Moreover we plan to couple the flow problem with transport, to fully exploit the freedom and flexibility given by the detailed description of the permeability of the layers. We remark that the coupling with transport which is natural thanks to the mixed conservative approximation of the velocity field.

7 Acknowledgements

We acknowledge the PorePy development team: Runar Berge, Inga Berre, Eirik Keilegavlen, and Ivar Stefansson. Finally, the authors warmly thanks Stefano Scialò for many fruitful discussions.