An unfitted finite element method for the Darcy problem in a fracture network

TL;DR

This paper introduces an unfitted finite element method for simulating Darcy flow in fracture networks, allowing arbitrary fracture cuts through the mesh and simplifying junction coupling while maintaining optimal convergence.

Contribution

It extends the Hughes--Masud stabilized formulation with penalty terms for interface conditions, enabling accurate and efficient modeling of complex fracture networks without mesh fitting.

Findings

Achieves optimal convergence order in numerical experiments.

Handles both planar and curvilinear fractures seamlessly.

Simplifies junction coupling by over-penalizing pressure continuity.

Abstract

The paper develops an unfitted finite element method for solving the Darcy system of equations posed in a network of fractures embedded in a porous matrix. The approach builds on the Hughes--Masud stabilized formulation of the Darcy problem and the trace finite element method. The system of fractures is allowed to cut through the background mesh in an arbitrary way. Moreover, the fractures are not triangulated in the common sense and the junctions of fractures are not fitted by the mesh. To couple the flow variables at multiple fracture junctions, we extend the Hughes--Masud formulation by including penalty terms to handle interface conditions. One observation made here is that by over-penalizing the pressure continuity interface condition one can avoid including additional jump terms along the fracture junctions. This simplifies the formulation while ensuring the optimal convergence order of the method. The application of the trace finite element allows to treat both planar and curvilinear fractures with the same ease. The paper presents convergence analysis and assesses the performance of the method in a series of numerical experiments. For the background mesh we use an octree grid with cubic cells. The flow in the fracture can be easily coupled with the flow in matrix, but we do not pursue the topic of discretizing such coupled system here.

Figures

Tables

Table 1. Table 1: Errors norms and convergence rates for Example  1 with α = 20 o 𝛼 superscript 20 𝑜 \alpha=20^{o} , β = 0 o 𝛽 superscript 0 𝑜 \beta=0^{o} .

#d.o.f.	${\Vert 𝐮-{𝐮}_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^{\infty }}$	rate
384	7.606e-2		4.602e-3		2.612e-2
1728	3.779e-2	1.01	1.371e-3	1.75	1.250e-2	1.06
7904	2.081e-2	0.86	3.925e-4	1.80	6.118e-3	1.03
33072	1.095e-2	0.93	1.097e-4	1.84	3.006e-3	1.02

Table 2. Table 2: Errors norms and convergence rates for Example  1 with α = 24 o 𝛼 superscript 24 𝑜 \alpha=24^{o} , β = 4 o 𝛽 superscript 4 𝑜 \beta=4^{o} .

#d.o.f.	${\Vert 𝐮-{𝐮}_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^{\infty }}$	rate
354	9.451e-2		6.023e-3		3.738e-2
1766	4.446e-2	1.09	1.518e-3	1.99	2.310e-2	0.70
7320	1.926e-2	1.20	2.879e-4	2.39	1.026e-2	1.17
30388	8.879e-3	1.13	8.149e-5	1.82	4.890e-3	1.07

Table 3. Table 3: Errors norms and convergence rates for Example  1 with α = 20 o 𝛼 superscript 20 𝑜 \alpha=20^{o} , β = 0 o 𝛽 superscript 0 𝑜 \beta=0^{o} and immersed part of the boundary.

#d.o.f.	${\Vert 𝐮-{𝐮}_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^{\infty }}$	rate
280	7.016e-2		4.275e-3		2.459e-2
1422	3.331e-2	0.93	1.182e-3	1.60	1.104e-2	0.99
6612	1.821e-2	0.88	3.419e-4	1.79	5.805e-3	0.92
27924	9.575e-3	0.96	9.388e-5	1.86	3.015e-3	0.95

Table 4. Table 4: Errors norms and convergence rates for the example  2

#d.o.f.	${\Vert 𝐮-{𝐮}_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^{\infty }}$	rate
556	2.250e-0		2.273e-1		6.005e-1
2332	5.978e-1	1.91	5.392e-2	2.07	1.593e-1	1.91
9532	1.559e-1	1.94	1.372e-2	1.97	4.121e-2	1.95
38212	4.907e-2	1.66	3.192e-3	2.10	9.613e-3	2.10

Table 5. Table 5: Errors norms and convergence rates for the example  3

#d.o.f.	${\Vert 𝐮-{𝐮}_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^2}$	rate	${\Vert p-p_h\Vert }_{L^{\infty }}$	rate
560	6.979e-2		1.749e-2		3.762e-2
2300	2.042e-2	1.77	3.775e-3	2.21	1.016e-2	1.89
9096	6.321e-3	2.13	8.759e-4	2.65	2.328e-3	2.68
35528	2.626e-3	1.60	2.154e-4	2.55	5.933e-4	2.48

Equations

$$\left\{\begin{aligned} {\color[rgb]{0,0,0}K_{i}^{-1}}\mathbf{u}_{i}+\nabla_{\Gamma}p_{i}&=\mathbf{f}_{i}\\ \mbox{\rm div}_{\Gamma}\mathbf{u}_{i}&=g\\ \mathbf{u}_{i}\cdot\mathbf{n}_{i}&=0\end{aligned}\right.\quad\text{in}~{}\Gamma_{i},\quad i=1,\dots,N,$$

$$\sum _{k=1}^{M_e}{\mathbf{u}}_{i_k}\cdot {\mathbf{m}}_{i_k}=0\text{on}e.$$

$$p_{i_1}=\cdots =p_{i_{M_e}}\text{on}e.$$

$$\{\begin{array}{rl}{\mathbf{m}}_i\cdot {\mathbf{u}}_i& ={\varphi }_i\text{on}{\partial \Gamma }_N\cap {\partial \Gamma }_i,i=1,\dots ,N\\ p& =p_D\text{on}{\partial \Gamma }_D.\end{array}$$

$$V_h:=\{v\in C(\Omega )|v{|}_S\in Q_1\forall S\in {\mathcal{T}}_h\},\text{with}{Q}_1=\text{mbox}span\{1,x_1,x_2,x_3,x_1{x}_2,x_1{x}_3,x_2{x}_3,x_1{x}_2{x}_3\}.$$

$${\Omega }_h^i=\bigcup \{T\in {\mathcal{T}}_h:T\cap {\Gamma }_i\ne \varnothing \},$$

$$V_h^i:=\{u\in C({\Omega }_h^i)|\exists v\in V_h\text{such that }u=v{|}_{{\Omega }_h^i}\}.$$

$$Q_h=\underset{i=1}{\overset{N}{⨂}}{V}_h^i\text{and}{\mathbf{U}}_h=\underset{i=1}{\overset{N}{⨂}}[V_h^i{]}^3.$$

$$(\text{mbox}{\mathrm{div}}_{\Gamma }\mathbf{v},q{)}_{{\Gamma }_i}=-(\mathbf{v},{\nabla }_{\Gamma }q{)}_{{\Gamma }_i}+({\mathbf{m}}_i\cdot \mathbf{v},q{)}_{{\partial \Gamma }_i}$$

$$({\color[rgb]{0,0,0}K_{i}^{-1}}\mathbf{u}+\nabla_{\Gamma}p,\mathbf{v})_{\Gamma_{i}}+(\mbox{\rm div}_{\Gamma}\mathbf{u},q)_{\Gamma_{i}}+\frac{1}{2}({\color[rgb]{0,0,0}K_{i}^{-1}}\mathbf{u}+\nabla_{\Gamma}p,-\mathbf{v}+{\color[rgb]{0,0,0}K_{i}}\nabla_{\Gamma}q)_{\Gamma_{i}}=(g,q)_{\Gamma_{i}}+\frac{1}{2}(\mathbf{f},-\mathbf{v}+{\color[rgb]{0,0,0}K_{i}}\nabla_{\Gamma}q)_{\Gamma_{i}}$$

$$({\color[rgb]{0,0,0}K_{i}^{-1}}\mathbf{u},\mathbf{v})_{\Gamma_{i}}+(\nabla_{\Gamma}p,\mathbf{v})_{\Gamma_{i}}-(\nabla_{\Gamma}q,\mathbf{u})_{\Gamma_{i}}+({\color[rgb]{0,0,0}K_{i}}\nabla_{\Gamma}p,\nabla_{\Gamma}q)_{\Gamma_{i}}+2(\mathbf{m}_{i}\cdot\mathbf{u},q)_{{\partial\Gamma}_{i}}=2(g,q)_{\Gamma_{i}}+(\mathbf{f},-\mathbf{v}+{\color[rgb]{0,0,0}K_{i}}\nabla_{\Gamma}q)_{\Gamma_{i}}.$$

$$({\color[rgb]{0,0,0}K_{i}^{-1}}\mathbf{u},\mathbf{v})_{\Gamma_{i}}+(\nabla p,\mathbf{v})_{\Gamma_{i}}-(\nabla q,\mathbf{u})_{\Gamma_{i}}+({\color[rgb]{0,0,0}K_{i}}\nabla p,\nabla q)_{\Gamma_{i}}+2(\mathbf{m}_{i}\cdot\mathbf{u},q)_{{\partial\Gamma}_{i}}=2(g,q)_{\Gamma_{i}}+(\mathbf{f},-\mathbf{v}+{\color[rgb]{0,0,0}K_{i}}\nabla q)_{\Gamma_{i}}.$$

$$({\color[rgb]{0,0,0}K^{-1}}\mathbf{u},\mathbf{v})_{\Gamma}+(\nabla p,\mathbf{v})_{\Gamma}-(\nabla q,\mathbf{u})_{\Gamma}+({\color[rgb]{0,0,0}K_{i}}\nabla p,\nabla q)_{\Gamma}+\sum_{e\in E}\frac{2}{M_{e}}\sum_{k=1}^{M_{e}-1}\sum_{\ell={k+1}}^{M_{e}}(\mathbf{m}_{i_{k}}\cdot\mathbf{u}_{i_{k}}-\mathbf{m}_{i_{\ell}}\cdot\mathbf{u}_{i_{\ell}},q_{i_{k}}-q_{i_{\ell}})_{e}\\ =2(g,q)_{\Gamma}+(\mathbf{f},-\mathbf{v}+{\color[rgb]{0,0,0}K_{i}}\nabla q)_{\Gamma}-2(\psi,q)_{{\partial\Gamma}_{N}}$$

$$\sum_{i=1}^{M}a_{i}b_{i}=\frac{1}{M}\left(\big{(}\sum_{i=1}^{M}a_{i}\big{)}\big{(}\sum_{i=1}^{M}b_{i}\big{)}+\sum_{i=1}^{M-1}\sum_{j={i+1}}^{M}(a_{i}-a_{j})(b_{i}-b_{j})\right)$$

$$a({\mathbf{u}}_h,p_h;{\mathbf{v}}_h,q_h)=f({\mathbf{v}}_h,q_h)$$

$$a(\mathbf{u},p;\mathbf{v},q)$$

$$+\underset{\text{normal volume stabilization}}{\underbrace{\sum _{i=1}^N{\rho }_{u}h({\mathbf{n}}_i\cdot \nabla {\mathbf{u}}_i,{\mathbf{n}}_i\cdot \nabla {\mathbf{v}}_i{)}_{{\Omega }_h^i}+\sum _{i=1}^N{\rho }_{p}h({\mathbf{n}}_i\cdot \nabla p_i,{\mathbf{n}}_i\cdot \nabla q_i{)}_{{\Omega }_h^i}}},$$

$$f(\mathbf{v},q)$$

$$\sum _{i=1}^N\frac{{\rho }_i}{h^2}(p_i-p_D,q_i{)}_{{\Gamma }_i\cap {\Gamma }_D}$$

$$\text{mbox}dist({\Gamma }_i,{\Gamma }_{h,i})\le c{h}_{\mathrm{loc}}^2,|\mathbf{n}(\mathbf{x})-{\mathbf{n}}_h(\tilde{\mathbf{x}})|\le c{h}_{\mathrm{loc}},$$

$$Q=\underset{i=1}{\overset{N}{⨂}}{H}^1({\Omega }_h^i)\cap H^1({\Gamma }_i)\text{and}\mathbf{U}=\underset{i=1}{\overset{N}{⨂}}[H^1({\Omega }_h^i){]}^3,Q_h\subset Q,{\mathbf{U}}_h\subset \mathbf{U}.$$

$$\|\mathbf{v}\|_{\ast}^{2}=({\color[rgb]{0,0,0}K^{-1}}\mathbf{v},\mathbf{v})_{\Gamma}+\sum_{i=1}^{N}\rho_{u}h\|\mathbf{n}_{i}\cdot\nabla\mathbf{v}_{i}\|^{2}_{\Omega_{h}^{i}},\quad\|q\|_{\ast}^{2}=\|{\color[rgb]{0,0,0}K}\nabla q\|^{2}_{\Gamma}+\sum_{e\in E}\frac{\rho_{e}}{h^{2}}\sum_{k=1}^{M_{e}-1}\sum_{\ell=k}^{M_{e}}\|q_{i_{k}}-q_{i_{\ell}}\|^{2}_{e}+\sum_{i=1}^{N}\rho_{p}h\|\mathbf{n}_{i}\cdot\nabla q_{i}\|^{2}_{\Omega_{h}^{i}},$$

$$\Vert \mathbf{v},q\Vert =(\Vert \mathbf{v}{\Vert }_{\ast }^2+\Vert q{\Vert }_{\ast }^2{)}^{\frac{1}{2}}.$$

$$\Vert \mathbf{v},q{\Vert }^2\le a(\mathbf{v},q;\mathbf{v},q)\forall \mathbf{v}\in \mathbf{U},q\in Q,$$

$$a(\mathbf{u},p;\mathbf{v},q)\le 2\Vert \mathbf{u},p\Vert \Vert \mathbf{v},q\Vert \forall \mathbf{u},\mathbf{v}\in \mathbf{U},p,q\in Q.$$

$$\Vert q{\Vert }_{{\Gamma }_i}\le C({\Gamma }_i,e)(\Vert {\nabla }_{\Gamma }q{\Vert }_{{\Gamma }_i}+\Vert q{\Vert }_e)q\in H^1({\Gamma }_i),e\in {\partial \Gamma }_i,\text{mbox}{meas}_1(e)>0,$$

$$\Vert q{\Vert }_e\le C({\Gamma }_i,e)(\Vert {\nabla }_{\Gamma }q{\Vert }_{{\Gamma }_i}+\Vert q{\Vert }_{{\Gamma }_i})q\in H^1({\Gamma }_i),e\in {\partial \Gamma }_i,\text{mbox}{meas}_1(e)>0,$$

$$\Vert q{\Vert }_{\Gamma }\le C{\left(\Vert \nabla q{\Vert }_{\Gamma }^2+\sum _{e\in E}\frac{{\rho }_e}{h^2}\sum _{k=1}^{M_e-1}\sum _{\ell =k}^{M_e}\Vert q_{i_k}-q_{i_{\ell }}{\Vert }_e^2\right)}^{\frac{1}{2}}\le C\Vert q{\Vert }_{\ast }\forall q\in Q^0,$$

$$\begin{array}{rl}f(\mathbf{v},q)& \le 2\Vert g{\Vert }_{\Gamma }\Vert q{\Vert }_{\Gamma }+(\Vert K^{\frac{1}{2}}\mathbf{f}{\Vert }_{\Gamma }\Vert K^{-\frac{1}{2}}\mathbf{v}{\Vert }_{\Gamma }+\Vert K^{\frac{1}{2}}\mathbf{f}{\Vert }_{\Gamma }\Vert K^{\frac{1}{2}}\nabla q{\Vert }_{\Gamma })-2\Vert \psi {\Vert }_{{\partial \Gamma }_N}\Vert q{\Vert }_{\partial \Gamma }\\ & \le C\left((\Vert g{\Vert }_{\Gamma }+\Vert \psi {\Vert }_{{\partial \Gamma }_N}+\Vert K^{\frac{1}{2}}\mathbf{f}{\Vert }_{\Gamma })\Vert q{\Vert }_{\ast }+\Vert K^{\frac{1}{2}}\mathbf{f}{\Vert }_{\Gamma }\Vert \mathbf{v}{\Vert }_{\ast }\right).\end{array}$$

$$\Vert {\mathbf{u}}_h{\Vert }_{\ast }+\Vert p_h{\Vert }_{\ast }\le C\left(\Vert g{\Vert }_{\Gamma }+\Vert \psi {\Vert }_{{\partial \Gamma }_N}+\Vert \mathbf{f}{\Vert }_{\Gamma }\right),$$

$$E{\mathbf{u}}_i={\mathbf{u}}_i\text{on}{\Gamma }_i,E{p}_i=p_i\text{on}{\Gamma }_i\text{and}\Vert E{\mathbf{u}}_i{\Vert }_{H^1({\hat{\Gamma }}_i)}\le C\Vert {\mathbf{u}}_i{\Vert }_{H^1({\Gamma }_i)},\Vert E{p}_i{\Vert }_{H^3({\hat{\Gamma }}_i)}\le C\Vert p_i{\Vert }_{H^3({\Gamma }_i)}$$

Full text

Keywords: fractured porous media, Darcy, Trace FEM, unfitted meshes, octree grid

1 Introduction

Numerical modelling of a flow in a fractured porous medium is a standard problem in geosciences and reservoir simulation [14, 37]. While the literature on this topic is overwhelming (see, e.g., [3, 6, 11, 16, 28] for a snapshot of recent research), the problem of developing an accurate and effective numerical method for a complex network of fractures still constitutes a challenge. The present paper contributes to the topic by introducing a finite element method for the Darcy problem posed in a system of intersecting fractures represented by a set of 2D surfaces embedded in a bulk domain. The enabling feature of the method is that it solely uses the background triangulation of the bulk domain (i.e., a tessellation in simplexes or more general polytopes) which is completely independent of the fracture network. Moreover, it does not require any 2D mesh fitted to the fracture surfaces or their intersections.

Application of geometrically unfitted finite element methods for the modelling of flow and transport in fractured porous medium have been addressed recently in a number of publications; see, e.g., [5, 17, 27]. Developments most closely related to the approach taken in the present paper are those found in [12, 18, 25]. Thus in [18] the authors consider a low order Raviart-Thomas finite element method for the Darcy flow on a 1D network of fractures. Although a triangulation of each fracture surface were build, these triangulations do not match the fracture intersection points, and the authors applied XFEM methodology to handle discontinuities in the solution over the junctions. The recent paper [20] reviews this and other numerical approaches, where a different degree of the conformity of fracture meshes at junction interfaces is assumed. However, triangulating each fracture branch can be itself a demanding task for large networks or complex geometry. Hence, the next level of nonconformity is to abandon the triangulation of the fracture in the usual sense and to discretize the flow problem along the fracture network only with the help of degrees of freedom tailored to the ambient mesh in the matrix. This ambient (background) mesh should be independent of the embedded fracture network. This approach was first taken in [31] to discretize scalar elliptic PDEs posed on surfaces and later it evolves to become the Trace FEM methodology [33] and a part of the Cut FEM [9]. Trace FEM for the transport and diffusion of a contaminant in a fractured porous media was recently developed in [12]. For the Darcy problem, the Trace FEM was first studied in [25]. In that paper, the authors considered the Darcy problem posed on a surface embedded in a bulk tetrahedra grid. They use a variant of Hughes–Masud week formulation to solve for the pressure and tangential velocity.

Following [25], we apply Trace FEM in combination with a variant of Hughes–Masud week formulation. The novelty of the present work is two-fold. First, for the ambient mesh we consider octree Cartesian grids, which can be easily adapted. Second and more importantly, we assume intersecting piecewise smooth surfaces (representing branching fractures), while in the previous work only closed smooth manifolds were considered. The branching leads to discontinuous fluxes and only piecewise smooth pressure field. Handling those without mesh fitting, but preserving optimal convergence order, is not straightforward. In the paper, this is achieved by allowing discontinuous velocity and pressure fields in background cells intersected by the fracture junctions and by including a penalty term. This treatment of fracture junctions draw an analogy with the Nitsche-XFEM method of Hansbo and Hansbo [24] for interface problems and more general with CutFEM [9]. One interesting difference, however, is that we use another scaling for the penalty term and skip certain consistency terms (typical for discontinuous Galerkin FEM and Nitsche’s method) along the junctions without sacrificing optimal asymptotic accuracy. The paper includes both numerical analysis and computational assessment of the method.

The remainder of the paper is organized as follows. Section 2 defines the mathematical model. Section 3 introduces the finite element method. Section 4 presents the convergence analysis of the method. Section 5 collects the results of several numerical experiments that illustrate the analysis and the performance of the method.

2 Mathematical model

Assume a piecewise smooth surface $\Gamma\subset\Omega$ embedded in the given bulk domain $\Omega\subset\mathbb{R}^{3}$. The surface $\Gamma$ represents a 2D fracture network and consists of several connected components $\overline{\Gamma}=\cup_{i=1}^{N}\overline{\Gamma}_{i}$, where each $\Gamma_{i}$ is smooth orientable surface without self-intersections. For the purpose of analysis, we shall assume that each $\Gamma_{i}$ is a subdomain of a larger $C^{2}$–smooth surface $\widehat{\Gamma}_{i}$, such that $\partial\widehat{\Gamma}_{i}\cap\Omega=\emptyset$ and $\partial\Gamma_{i}$ is piecewise smooth and Lipschitz as a curve in $\widehat{\Gamma}_{i}$. The individual components $\Gamma_{i}$ may intersect only by a curve, i.e. $\mbox{meas}_{2}(\overline{\Gamma}_{i}\cap\overline{\Gamma}_{j})=0$ for $i\neq j$, and also ${\Gamma}_{i}\cap{\Gamma}_{j}=\emptyset$, for $i\neq j$ (this condition means that parts of a fracture separated by a junction are treated as different components). Further $\mathbf{n}$ is a unit normal vector defined everywhere on $\Gamma$ except junction interfaces. We shall write $\mathbf{n}_{i}$ for $\mathbf{n}$ on $\Gamma_{i}$ and similar for other vector and scala fields defined on $\cup_{i=1}^{N}\Gamma_{i}$.

Modeling fractures as 2D interfaces for flow in porous media has been considered in many places in the literature; see, e.g., [1, 2, 19, 29]. In this framework, the flow along the fracture component ${\Gamma}_{i}$ is described in terms of tangential velocity field $\mathbf{u}_{i}(\mathbf{x})$, having the physical meaning of the flow rate through the cross-section of the fracture, and pressure field $p_{i}(\mathbf{x})$, $\mathbf{x}\in{\Gamma}_{i}$. The steady state flow in $\Gamma$, is governed by the Darcy systems

[TABLE]

together with interface and boundary conditions specified below. In (1) and further in the text, $\nabla_{\Gamma}$ and $\mbox{\rm div}_{\Gamma}$ denote the surface tangential gradient and divergence operators; $g$ stands for the source term, which is typically due to the fluid exchange with the porous matrix (not treated in this paper); $\mathbf{f}_{i}$ is an exterior force per unit area, $\mathbf{f}_{i}$ is tangential to $\Gamma_{i}$; $K_{i}$ denotes the permeability tensor along the fracture; all $K_{i}$ are symmetric and such that for any tangential vector field $\mathbf{v}$, i.e. $\mathbf{v}\cdot\mathbf{n}_{i}=0$, it holds $\mathbf{v}^{T}K_{i}\mathbf{v}\geq\zeta_{i}|\mathbf{v}|^{2}$ with some $\zeta_{i}>0$, and $\mathbf{n}^{T}_{i}K_{i}\mathbf{v}=0$. Hence, $K_{i}^{-1}\mathbf{v}$ is well defined for a tangential field $\mathbf{v}$. Note that fracture aperture can be included in $K_{i}$ by scaling; see, e.g. [1].

When $\Gamma$ is piecewise smooth, we need further conditions on the edges (fracture junctions). Consider an edge $e$ shared by $M_{e}$ smooth components $\Gamma_{i_{k}}$, $k=1,\dots,M_{e}$. Here and in the rest of the paper, $\{i_{k}\}_{k=1,\dots,M_{e}}$ denote the subset of indexes from $\{1,\dots,N\}$, which is specific for each given $e$. Denote by $\mathbf{m}_{i}$ the normal vector on ${\partial\Gamma}_{i}$ in the plane tangential to $\Gamma_{i}$ and pointing outward. The conservation of fluid mass yields

[TABLE]

The second interface condition is the continuity of pressure over $e$,

[TABLE]

Denote by $E$ the collections of all fracture junctions. It is reasonable to assume that $E$ is a finite set and $0<\mbox{meas}_{1}(e)<+\infty$ for any $e\in E$.

Finally, we prescribe the pressure boundary condition on ${\partial\Gamma}_{D}$ and the flux boundary condition on ${\partial\Gamma}_{N}$, respectively, with $\overline{{\partial\Gamma}}=\overline{{\partial\Gamma}_{D}}\cup\overline{{\partial\Gamma}_{N}}$,

[TABLE]

3 Finite element method

First we assume a tessellation $\mathcal{T}_{h}$ of the bulk domain $\Omega$ (matrix). $\mathcal{T}_{h}$ can be a consistent subdivision into shape-regular tetrahedra. In this paper, we consider Cartesian background mesh with cubic cells. We allow local refinement of the mesh by sequential division of any cubic cell into 8 cubic subcells. This leads to a grid with an octree hierarchical structure. This mesh gives the tessellation $\mathcal{T}_{h}$ of the bulk domain $\Omega$, $\overline{\Omega}=\cup_{T\in\mathcal{T}_{h}}\overline{T}$. We allow the fracture network $\Gamma\subset\Omega$ to cut through this mesh in an arbitrary way. For the purpose of analysis, we shall assume that the cells cut by $\Gamma$ have a quasi-uniform size with the characteristic size $h$.

Consider now the ambient finite element space of all piecewise trilinear continuous functions with respect to the bulk octree mesh $\mathcal{T}_{h}$:

[TABLE]

For every fracture $\Gamma_{i}$ in the network $\Gamma$ we define the subdomain of $\Omega$ consisting of all cells cut by $\Gamma_{i}$,

[TABLE]

and define the restriction of $V_{h}$ on $\Omega^{i}_{h}$, i.e. the space of piecewise trilinear continuous functions on $\Omega^{i}_{h}$,

[TABLE]

Our trial and test finite element spaces are built from $V_{h}^{i}$: We define the pressure space and velocity spaces

[TABLE]

According to the Trace FEM approach, the finite element solutions of (1)–(4) will be given by traces of functions from $Q_{h}$ and $\mathbf{U}_{h}$ on $\Gamma$, but the finite element formulation will be written in terms of function defined on $\bigcup_{i=1}^{N}\Omega^{i}_{h}$. Hence, the method leads to a system of algebraic equations for standard nodal degrees of freedom in the ambient mesh $\mathcal{T}_{h}$.

Further we use the notation $(\cdot,\cdot)_{Q}$ for the $L^{2}$ scalar product over a domain $Q$, which can be a 3D, 2D or 1D manifold on different occasions. For example, with this notation, the Green formula on $\Gamma_{i}$ reads:

[TABLE]

for any smooth tangential vector field $\mathbf{v}$ and scalar function $q$ on $\Gamma_{i}$.

The proposed finite element formulation extends the stabilized mixed formulation for the Darcy problem originally introduced in [30] for the planar domains. The key observation here is that the smooth solution to (1)–(4) satisfies the identity

[TABLE]

for all $q\in H^{1}(\Gamma_{i})$, $\mathbf{v}\in L^{2}(\Gamma_{i})$ and $i=1,\dots,N$. We now set $q=0$ on ${\partial\Gamma}_{D}$ and apply (7). After simple calculations this gives

[TABLE]

One further helpful observation is that $p$ and $q$ can be identified with their normal extensions to a neighborhood of $\Gamma_{i}$ (for each $i$). This identification (which is assumed further in the paper) implies the equality $\nabla_{\Gamma}p=\nabla p$, which can be further used in (8) to yield

[TABLE]

This corresponds to so-called full gradient formulation of the surface PDES; see [13, 35]. The full gradient formulation exploits the embedding of $\Gamma$ in the ambient Euclidian space and, in general, provides extra stability for a finite element method based on external elements. The formulation is consistent for any ambient finite element method, which aims to approximate the surface solution together with its normal extension. In the context of the surface Darcy problem, the full gradient formulation was used in [25]. We finally sum up equalities (9) for all $i=1,\dots,N$ and use the interface condition (2) and the boundary condition for fluxes from (4) to conclude that any smooth solution of (1)–(4) satisfies

[TABLE]

for any $q\in\bigotimes_{i=1}^{N}H^{1}(\Gamma_{i})$ such that $q=0$ on ${\partial\Gamma}_{D}$ and $\mathbf{v}\in L^{2}(\Gamma)^{3}$. To handle the sum of the edge terms, we used (2) and the identity

[TABLE]

for any $a_{i},b_{i}\in\mathbb{R}$.

Our finite element method is based on the equality (10). Note that we may assume that $K_{i}$ is extended to be symmetric positive definite in $\mathbb{R}^{3}$, rather than only on the tangential space, since this does not affect any quantities in (10), but would be helpful, when we proceed with the finite element formulation. To approximate pressure, we use finite element functions from $Q_{h}$, which are discontinuous across $e\in E$. Therefore, we add a penalty term to our formulation to weekly enforce the pressure continuity condition from (4). Furthermore, we over-penalize this condition, by choosing a different scaling of the penalty parameter compared to the standard Nitsche’s [24] or discontinuous Galerkin methods [4]. It turns out that the over-penalization allows one to skip other edge terms in the finite formulation. This greatly simplifies the method while keeping the consistency order optimal. Summarizing, the finite element method reads: Find $\mathbf{u}_{h}\in\mathbf{U}_{h}$ and $p_{h}\in Q_{h}$ such that p_{h}|_{\partial\Gamma_{D}}={\color[rgb]{0,0,0}I_{h}^{b}}(p_{D}) and

[TABLE]

for all $\mathbf{v}_{h}\in\mathbf{U}_{h}$ and $q_{h}\in Q_{h}$ such that $q_{h}|_{\partial\Gamma_{D}}=0$, with

[TABLE]

Here $\rho$’s are tunable parameters, which we set (in both analysis and experiments) to be equal to $1$; $I_{h}^{b}(p_{D})$ is the interpolation of the boundary condition, which we define by extending pressure values from $\partial\Gamma_{D}$ along normal directions in ${\partial\Omega}$ to the corresponding nodal values from $\Omega_{h}^{i}\cap{\partial\Omega}$.

Remark 3.1** (Normal volume stabilization).**

We briefly discuss the “normal volume stabilization” terms in (12). The term involve the extension of the normal vector to $\Omega_{h}^{i}$, which can be defined as $\mathbf{n}_{i}(\mathbf{x})=\nabla\mbox{dist}(\mathbf{x},\widehat{\Gamma}_{i})$. Assuming that the mesh is fine enough to resolve the (curvilinear) geometry, this definition gives the meaning to the normal vector in all $\Omega_{h}^{i}$, including mesh cells cut by ${\partial\Gamma}_{i}$. Next, we note that the normal volume stabilization terms vanish for the solution $\mathbf{u}$, $p$ of the Darcy equations (1), because we assume the normal extension of the solution off the fracture components. Finally, these terms are included, following [10, 21], to ensure that algebraic properties of the resulting linear systems are insensitive to the position of $\Gamma$ against the background mesh. Indeed, if $\rho_{u}=0$ or $\rho_{p}=0$, then for a natural nodal basis in $\mathbf{U}_{h}$ and $Q_{h}$, small cuts of the background elements by the surface may lead to arbitrarily small diagonal entries in the resulting matrix. The stabilization terms in (12) eliminate this problem since for the choice of $\rho_{u}=O(1)$ and $\rho_{p}=O(1)$ they allow to get control over the $L^{2}(\Omega_{h}^{i})$-norms of $\mathbf{v}_{h}\in(V_{h}^{i})^{3}$, $q_{h}\in V_{h}^{i}$ by the problem induced norms. The analysis of this acquired algebraic stability can be found at several places in the literature, e.g. [10, 21] for the Laplace–Beltrami problem or [32] for the surface Stokes problem, so we omit repeating it here.

Remark 3.2** (Overpenalty).**

In the framework of discontinuous Galerkin (DG) methods, a technique similar to the overpenalty used to enforce pressure continuity here is known as a superpenalty; see [4] and references there. Compared to the superpenalty technique in DG FEM, we have a weaker dependence of the penalty parameter on the negative power of $h$, which is beneficial for the condition number of the resulting matrices. On the analysis size, the superpenalty DG method exploits the availability of a continuous finite element interpolant across the element edges, which is not the case here. As a consequence, to show a suitable consistency bound, we have to apply a different argument comparing to the analysis of the superpenalty DG method. This results in the extra smoothness assumption for pressure solution, e.g. $p$ is from $H^{3}$ on every $\Gamma_{i}$ rather than from $H^{2}$.

Remark 3.3** (Internal parts of $\partial\Gamma$).**

For the finite element formulation in (12) we assumed that ${\partial\Gamma}\subset{\partial\Omega}$. If $\partial\Gamma_{D}$ has a part strictly inside $\Omega$, then the pressure boundary condition can be enforced by including additional penalty term of the form

[TABLE]

to the finite element formulation (12). Internal boundaries with prescribed fluxes, i.e. $\partial\Gamma_{N}\subset\Omega$, do not affect the formulation in (12).

3.1 Numerical integration

The finite element formulation (12) requires computing surface integrals. If $\Gamma_{i}$ is a planar component, then numerical integration is straightforward. For a curvilinear $\Gamma$, in general, we need to know a (local) parametrization of the surface to compute integrals in (12). For implicitly given surfaces (for example, for surfaces defined as the zero of a distance function), the numerical integration is a more subtle issue; see, e.g., [34]. In the present paper, for numerical tests with curvilinear surfaces we compute surface integrals by using a polygonal second order approximation of $\Gamma_{i}$, denoted by $\Gamma_{h,i}$. We construct $\Gamma_{h,i}$ as follows. For $\Gamma_{i}$ let $\phi$ be a Lipshitz-continuous level set function, such that $\phi(\mathbf{x})=0$ on $\Gamma_{i}$. We set $\phi_{h}=I(\phi)$, a nodal interpolant of $\phi$ by a piecewise trilinear continuous function with respect to the octree grid $\mathcal{T}_{h}$. Further, consider the zero level set of $\phi_{h}$, $\widetilde{\Gamma}_{h,i}:=\{\mathbf{x}\in\Omega\,:\,\phi_{h}(\mathbf{x})=0\}.$ If $\Gamma_{i}$ is smooth, then $\widetilde{\Gamma}_{h,i}$ is an approximation to $\Gamma_{i}$ in the following sense:

[TABLE]

where $\mathbf{x}$ is the closest point on $\Gamma_{i}$ for $\tilde{\mathbf{x}}\in\widetilde{\Gamma}_{h,i}$ and $h_{\rm loc}$ is the local mesh size. We note that in some applications, $\phi_{h}$ is computed from a solution of a discrete indicator function equation, without a direct knowledge of $\Gamma$.

Once $\phi_{h}$ is computed, we recover $\Gamma_{h,i}$ by the cubical marching squares method from [26] (a variant of the very well-known marching cubes method). The method provides a triangulation of $\widetilde{\Gamma}_{h}$ within each cube such that the global triangulation is continuous, the number of triangles within each cube is finite and bounded by a constant independent of $\widetilde{\Gamma}_{h,i}$ and a number of refinement levels. Moreover, the vertices of triangles from $\mathcal{F}_{h}$ are lying on $\widetilde{\Gamma}_{h,i}$. This final discrete surface $\Gamma_{h,i}$ is still an approximation of $\Gamma_{i}$ in the sense of (13). An example of bulk domain with embedded surface and background mesh is illustrated in Figure 1. Note that the resulting “triangulation” of $\Gamma_{h,i}$ is not shape regular. This is not a problem, since this triangulation is used only to define quadratures in the finite element method, while approximation properties of the method depend on the volumetric octree mesh.

4 Error analysis

The error analysis fits the standard finite element framework. Certain care and less standard arguments will be needed to show optimal order consistency for the formulation in (12). Interpolation results rely on approximation properties of the polynomial traces on smooth surfaces. We start with the definition of the norms used further and the proof of numerical stability.

4.1 Stability

For analysis, we need the broken Sobolev spaces

[TABLE]

The subspaces of functions from $Q$ and $Q_{h}$ vanishing on ${\partial\Gamma}_{D}$ are denoted by $Q^{0}$ and $Q_{h}^{0}$, respectively. The stability estimate of the method involves the following problem-dependent velocity and pressure norms:

[TABLE]

where $\mathbf{v}\in\mathbf{U}$, $q\in Q$. For the later expression to define a norm on $Q^{0}$ we shall assume that $\mbox{meas}_{1}({\partial\Gamma}_{D})>0$ and $\Gamma$ is connected. Otherwise, if $\mbox{meas}_{1}({\partial\Gamma}_{D})=0$, one uses a factor norm, since the pressure solution to (1)–(4) is defined only up to the addition of the hydrostatic constant mode. On the product space $\mathbf{U}\times Q$, we define

[TABLE]

Given the definitions above, one immediately checks the coercivity and continuity bounds for the finite element method bilinear form,

[TABLE]

and

[TABLE]

The repeating application of the Sobolev inequality

[TABLE]

and the trace inequality

[TABLE]

leads to the Poincare inequality on $\Gamma$:

[TABLE]

where $C$ depends on $\Gamma$, ${\partial\Gamma}_{D}$ and the permeability tensor $K$.

With the help of the Poincare, Cauchy-Schwartz, trace and triangle inequalities, one shows

[TABLE]

Finally, from (11), (14) and the estimate above we get the stability bound for the finite element solution

[TABLE]

where $C$ depends on $\Gamma$, ${\partial\Gamma}_{D}$ and $K$, but not on the position of $\Gamma$ over the background mesh. Penalty and volumetric stabilization terms (cf. Remark 3.1) in the definition of the norms on the left-hand side depend on $h$.

4.2 Error bound

For the error analysis, we shall assume that all fractures are planar, so that $\widehat{\Gamma}_{i}\simeq\mathbb{R}^{2}$ is just a plane. We believe that the error estimate below still holds for fractures with non-zero curvatures, but the analysis needs trace and extensions results for functions defined on submanifolds as in (19) and (21), which we do not find in the literature and including their proof would made this paper excessively technical.

As usual, the error analysis needs some extra regularity of the solution, namely $\mathbf{u}_{i}\in H^{1}(\Gamma_{i})^{3}$ and $p_{i}\in H^{3}(\Gamma_{i})$ for each $i=1,\dots,N$. We first observe that the finite element formulation (11) is consistent up to missing interface integrals. To see this, we give sense to the solution components $\mathbf{u}_{i}$ and $p_{i}$ not just on $\Gamma_{i}$ but in the neighborhood $\Omega^{i}_{h}$. To this end, we first consider Stein’s [36] extensions of $\mathbf{u}_{i}$, $p_{i}$ to some $E\mathbf{u}_{i}\in H^{1}(\widehat{\Gamma}_{i})^{3}$ and $Ep_{i}\in H^{3}(\widehat{\Gamma}_{i})$ such that

[TABLE]

with some finite $C$ depending only on $\Gamma_{i}$. Now we define normal extensions in $\Omega^{i}_{h}$, $p^{e}_{i}=Ep_{i}\circ\mathbf{p}$, where $\mathbf{p}$ is the closest point projection on $\widehat{\Gamma}_{i}$, and similar we define $\mathbf{u}_{i}^{e}$. From the properties of the normal extension we get $p^{e}_{i}\in H^{3}(\Omega_{h}^{i})$, $\mathbf{u}_{i}^{e}\in H^{1}(\Omega_{h}^{i})^{3}$ and for the norms it holds $\|p_{i}^{e}\|^{2}_{H^{3}(\Omega_{h}^{i})}\leq C\,h\,\|Ep_{i}\|^{2}_{H^{3}(\widehat{\Gamma}_{i})}$, $\|\mathbf{u}_{i}^{e}\|^{2}_{H^{1}(\Omega_{h}^{i})}\leq C\,h\,\|E\mathbf{u}_{i}\|^{2}_{H^{1}(\widehat{\Gamma}_{i})}$; see, e.g., [35]. Combining this with (19) gives

[TABLE]

with some finite $C$ depending only on $\Gamma_{i}$. If no confusion arises, we further identify $\mathbf{u}$ and $p$ with their extensions defined above.

For handling fracture junctions terms, we also consider the $h$-neighborhood of $\partial\Gamma_{i}$ in $\widehat{\Gamma}_{i}$, denoted by $\mathcal{O}(\partial\Gamma_{i})$. Then it holds (see, Lemma 4.10 in [15])

[TABLE]

We can always assume that $\mathcal{O}(\partial\Gamma_{i})$ is wide enough so that $T\cap\widehat{\Gamma}_{i}\subset\mathcal{O}(\partial\Gamma_{i})$ for all $T\in\mathcal{T}_{h}$ such that $T\cap\partial\Gamma_{i}\neq\emptyset$.

Now the normal volume stabilization terms make sense (and equals zero) for $\mathbf{u}$ and $p$, and we see that the piecewise smooth solution $\mathbf{u},p$ to the network Darcy problem (1)–(4) satisfies the equality

[TABLE]

with

[TABLE]

Using the Cauchy–Schwartz and triangle inequalities, the definition of the $\|\cdot\|_{\ast}$ norm on $Q$ and $|\mathbf{m}_{i}|=1$, one readily checks the upper bound,

[TABLE]

with a constant $C$ depending only on $\Gamma$.

We proceed with the interpolation bounds.

Lemma 4.1**.**

Let $\mathbf{u}\in\bigotimes_{i=1}^{N}H^{1}(\Gamma_{i})^{3}$, $p\in\bigotimes_{i=1}^{N}H^{3}(\Gamma_{i})$. Assume $h\leq h_{0}$, where $h_{0}$ may depend on $\Gamma$, then it holds

[TABLE]

with a constant $C$ independent on how $\Gamma$ intersects with $\mathcal{T}_{h}$.

Proof.

The analysis of the interpolation properties of the Trace FEM is commonly based on the local trace inequality, see e.g. [35, 10], which in our case takes the form:

[TABLE]

where $h_{T}=\text{diam}(T)$ and $C$ is independent of $T$, $v$ and how $\widehat{\Gamma}_{i}$ cuts through $T$. A quick proof of (25) consists in dividing the cubic cell $T$ into a finite number of regular tetrahedra and further applying Lemma 4.2 from [24] on each of these tetrahedra. For handling edge terms, we need the extension of (25) for curves cutting through the mesh. More precisely, we need the following inequality:

[TABLE]

where $C$ is independent of $T$, $v$ and how $\partial\Gamma_{i}$ cuts through $T$. We provide the proof of (26) in Appendix.

We recall that $\mathbf{u}_{i}\in H^{1}(\Gamma_{i})^{3}$ and $p_{i}\in H^{3}(\Gamma_{i})$ are identified with their extensions to $\Omega_{h}^{i}$ such that (20) holds. Let \mathbf{w}_{h}=I_{h}\big{(}\mathbf{u}\big{)}\in\mathbf{U}_{h}, $q_{h}=I_{h}(p)\in Q_{h}$ be the finite element (Clement) interpolants. Let us first treat the edge term in the definition of the $\|p-q_{h}\|_{\ast}$: Let $e\subset\partial\Gamma_{i}$ for $e\in E$. Using interpolation properties of bilinear polynomials and (26), we have for any $T\in\Omega_{h}^{i},~{}T\cap e\neq\emptyset$:

[TABLE]

Summing up the above inequality over all $T$ intersecting $e$ (the domain formed by all such cells is denoted by $\widetilde{\Omega}_{h}^{i}$), we get

[TABLE]

Here we used (20) (which remains true with $\widetilde{\Omega}_{h}^{i}$ and $\mathcal{O}(\partial\Gamma_{i})$ instead of ${\Omega}_{h}^{i}$ and $\Gamma_{i}$) and (21). For the rest of $\|p-q_{h}\|_{\ast}$ we use interpolation properties of bilinear polynomials, (25) and (20) to obtain

[TABLE]

and

[TABLE]

Estimates (27)–(29) lead to the desired bound on $\|p-q_{h}\|_{\ast}$. Similar to (28)–(29), we use interpolation properties of bilinear polynomials, (25) and (20), to obtain

[TABLE]

with a constant $C$ independent of $h$ and how $\Gamma$ intersects the background mesh.

∎

Since for each component $\Gamma_{i}$ we associate with $\mathbf{u}_{i}$ and $p_{i}$ there extensions to $\Omega_{h}^{i}$, the error functions $\mathbf{u}-\mathbf{u}_{h}$ and $p-p_{h}$ are well-defined as functions of $\mathbf{U}$ and $Q$. We are ready to prove the following convergence result.

Theorem 4.1**.**

Let $(\mathbf{u},p)$ be the solution of (1)–(4) and assume that $\mathbf{u}\in\bigotimes_{i=1}^{N}H^{1}(\Gamma_{i})^{3}$, $p\in\bigotimes_{i=1}^{N}H^{3}(\Gamma_{i})$. Let $(\mathbf{u}_{h},p_{h})\in\mathbf{U}_{h}\times Q_{h}$ be the solution of (11). The following discretization error bound holds:

[TABLE]

Here $\|\cdot\|_{k}$, $k=1,3$ denotes the broken Sobolev spaces norms for $\bigotimes_{i=1}^{N}H^{k}(\Gamma_{i})$; the constant $C$ depends on $\Gamma$, but not on how $\Gamma$ intersects the background mesh.

Proof.

Using the coercivity and consistency properties in (14) and (22) as well as continuity estimates for the $a$ and E forms (15), (23), we obtain, for arbitrary $(\mathbf{w}_{h},\xi_{h})\in\mathbf{U}_{h}\times Q_{h}$:

[TABLE]

Therefore, after cancellation and using the trace inequality (16) we get

[TABLE]

For $(\mathbf{w}_{h},\xi_{h})\in\mathbf{U}_{h}\times Q_{h}$ we take optimal finite element interpolants for the (normal extensions of the) solution \mathbf{w}_{h}=I_{h}\big{(}\mathbf{u}\big{)}, $q_{h}=I_{h}(p)$. Now, the triangle inequality, (31) and (24) leads to (30):

[TABLE]

∎

Remark 4.1** ($O(h^{2})$ convergence).**

In [25] the trace finite element method as in (12) applied on a smooth closed surface $\Gamma$ was proved to enjoy higher convergence in weaker norms for the pressure and fluxes. In our setting, this would mean the estimate

[TABLE]

for the pressure error. The convergence leverage argument, as usual, is based on the $H^{2}$- regularity estimate for the solution of the dual problem, which is the same system of Darcy equations in the fracture network in our case. However, we are not aware of a suitable regularity results for the case of intersecting fractures. Namely, we would need the estimate of the norm $\|p\|_{2}+\|\mathbf{u}\|_{1}\leq C\|\mathbf{f}\|_{1}+\|g\|_{\Gamma}$, again $\|\cdot\|_{k}$ are norms on the broken Sobolev spaces. Moreover, the studies in [8] of the Poisson problem posed in a domain with intersecting interfaces suggest that this higher regularity results might not hold in our case.

5 Numerical results and discussion

This section collects several numerical examples, which demonstrate the accuracy and capability of our unfitted finite element method. To verify the convergence rates of the method, we start with a few examples where exact solution is known. This includes the case of planar intersecting fractures and Darcy flow along curvilinear surfaces. Further we include an example of a pressure drop driven flow in a more complex network of fractures.

5.1 A multiple fracture problem with a synthetic solution

To test the convergence of the method, we first consider the example of an analytically prescribed solution on the two intersecting fractures; the test is built on an example from [7]. The setup is given below.

Example 1**.**

Consider $\Omega=(0,1)^{3}$. The initial configuration of the fracture network $\Gamma$ is given by the union of the two rectangles $\{(x,y,z)\in\Omega\mid x=0.5\}$ and $\{(x,y,z)\in\Omega\mid y=0.5\}$ so that $\Gamma=\bigcup\limits_{i=1}^{4}\Gamma_{i}$.

To define the exact solution $(\mathbf{u}_{i},p_{i})$ for the initial configuration, we first introduce the functions $t_{1}(\mathbf{x})=y+z-0.5$, $t_{2}(\mathbf{x})=x+z-0.5$, $t_{3}(\mathbf{x})=-y+z+0.5$ and $t_{4}(\mathbf{x})=-x+z+0.5$. The pressure and the Darcy velocity in each fracture component are given by

[TABLE]

where $i=1,..,4$ and $d_{x}=(0,1,0,-1),d_{y}=(1,0,-1,0)$. This $p$ and $\mathbf{u}$ satisfy (1)–(4) with $\mathbf{f}=0$ and

[TABLE]

To generate less regular intersections of the fractures and the junction line with the bulk mesh, we next perform the deformation of the fracture system by applying counterclockwise rotations by the angle $\alpha$ about the axis $x=z=0.5$ and by the angle $\beta$ about the axis $x=y=0.5$. The resulting fracture network is denoted by $\Gamma(\alpha,\beta)$. The corresponding change of variables is applied to prescribe the exact solution on $\Gamma(\alpha,\beta)$ using $(\mathbf{u}_{i},p_{i})$ defined above. For numerical experiments, we take $\alpha=20^{o}$, $\beta=0^{o}$ and $\alpha=24^{o}$, $\beta=4^{o}$.

We next consider a sequence of uniform tessellations of $\Omega$ into cubes with $h\in\{1/9,\,1/19,\,1/39,\,1/79\}$.

The trace of the second level ($h=1/19$) volumetric grid on $\Gamma(24^{o},4^{o})$ and the part of the volumetric grid intersected by the surface are illustrated in Fig. 2. The computed pressure and velocity for $\Gamma(20^{o},0^{o})$ are shown in Figure 3. Tables 1–2 present the error norms for the computed finite element solutions. We measure the error in the $L^{2}(\Gamma)$ and $L^{\infty}(\Gamma)$ for the pressure and $L^{2}(\Gamma)^{3}$ for the velocity. The results show close to the second order convergence for the pressure $L^{2}$ norm (although we are unable to prove it) and the first order convergence for the velocity. The $L^{\infty}$ norm of the error for the pressure also goes to zero as $O(h)$.

Finally, we consider the case when a part of the fracture’s boundary is immersed in the bulk as illustrated in Figure 3 (right). This corresponds to $\Gamma(20^{o},0^{o})$, with $\Gamma_{2}$ cut so that the immersed part of the boundary is vertical and the width of $\Gamma_{2}$ (i.e. the distance between the immersed boundary in the junction) is $0.25$. The background mesh does not fit the immersed boundary. Hence, we impose pressure Dirichlet boundary condition using the penalty term as described in Remark 3.3. The finite element error, reported in Table 3, appears to be almost unaffected by the presence of the immersed boundary.

5.2 Darcy flow over curvilinear surfaces

We now check if the fracture curvature influences convergence rates of the unfitted finite element method. To this end, we consider Darcy problem (1) defined on surface and on the torus. In both examples, $\Gamma$ is given by closed smooth surfaces (no junctions).

Example 2**.**

We consider $\Gamma=\{\mathbf{x}\in\mathbb{R}^{3}\mid\|\mathbf{x}\|_{2}=1\}$ embedded in $\Omega=(-2,\,2)^{3}$. The solution $(\mathbf{u},p)$ is given by

[TABLE]

with $a=12$. One verifies that $\mathbf{u}$ and $p$ satisfy (1) with $\mathbf{f}=0$, and $g=12p$, so that $g$ satisfies the compatibility condition $\int_{\Gamma}g\,\mathrm{d}\mathbf{s}=0$.

Example 3**.**

This example can be found in [25] for the Darcy flow along the torus surface. We consider $\Gamma=\{\mathbf{x}\in\Omega\mid r^{2}=x_{3}^{2}+(\sqrt{x_{1}^{2}+x_{2}^{2}}-R)^{2}\}$ embedded in $\Omega=(-1.6,1.6)\times(-1.6,1.6)\times(-0.8,0.8)$. We set $R=1$ and $r=0.5$. The solution $(\mathbf{u},p)$ to (1) with right-hand sides $g=0$ and

[TABLE]

is given by

[TABLE]

Again, we use a sequence of octree bulk grids. We start with the initial uniform grids with $h=1/4$ for Example 2 and $h=8/25$ for Example 3, which were further gradely refined towards the surfaces. Tables 4 and 5 show finite element errors and convergence rates for the computed solutions over several levels of refinement. We see that the convergence rates overall improve compared to the case with junctions, which is expected from the analysis. At the same time, the fracture bending does not affect the efficiency of the method, which is also well known property of the Trace FEM. The computed solutions and induced surface meshes are illustrated in Figure 4. Convergence rates for this test well agree with those reported in [25] for tetrahedra bulk mesh.

5.3 Pressure drop driven flow in a fracture network

The last example demonstrates the flexibility in applying the method for the case of more complex fracture networks.

Example 4**.**

We consider a fracture network consisting of $5$ components, both curvilinear and planar, and embedded in the bulk domain $\Omega=(-1,1)^{3}$. The fracture network is illustrated in Figure 5 (left), where each component $\Gamma_{i}$, $i=1,\dots,5$ has a distinct color. On the parts of $\partial\Gamma$ crossing the left and the right sides of the cube, we prescribe the Dirichlet pressure boundary conditions: $p=2$ for $\{\mathbf{x}\in\Gamma_{1}\,:\,x=-1\}$ and $p=0$ for $\{\mathbf{x}\in\Gamma_{2}\cup\Gamma_{4}\,:\,x=1\}$ and for $\{\mathbf{x}\in\Gamma_{5}:z=1\}$. On the rest of $\partial\Gamma$ we prescribe zero-flux conditions. Thus the boundary conditions define the pressure drop that drives the flow from left to right through the network. The computed solution is shown in Figure 5 (right).

Acknowledgments

The work of the first author (the method development, setup and analysis of numerical experiments) was supported by the Russian Science Foundation Grant 17-71-10173; the second author was supported by the NSF grant 1717516.

Appendix A Proof of the FE trace inequality (26)

The proof largely follows the arguments given in [23] to prove (25) and makes use of the following result found for example in [22]:

[TABLE]

for a bounded domain $\omega\subset\mathbb{R}^{n}$ with Lipschitz boundary.

Lemma A.1**.**

Let $T\in\mathcal{T}_{h}$. There exists an extension operator $R_{T}:H^{2}(T)\rightarrow H(\mathbb{R}^{3})$ such that $R_{T}v=v$ on $T$ and

[TABLE]

where the constant $C$ is independent of $T$ and $v$.

Proof.

We denote by $\hat{T}$ the reference cube of unit size. The center of $\hat{T}$ is placed at the origin. Then, we know ([36]) there exists an extension operator from $R:H^{1}(\hat{T})\rightarrow H_{0}^{1}(B_{2})$ such that $R\hat{v}=\hat{v}$ on $\hat{T}$ and

[TABLE]

where $B_{2}$ is the ball with radius 2 centered at the origin.

Let $F_{T}:\hat{T}\rightarrow T$ be the onto affine mapping and has the form $F_{T}(\hat{x})=B\hat{x}+b$. For any $v\in H^{2}(T)$ we define $\hat{v}\in H^{2}(\hat{T})$ in the following way: $\hat{v}(\hat{x})=v(F_{T}(\hat{x}))$. Then, the desired extension is given by $(R_{T}v)(x)=(R\hat{v})(F^{-1}_{T}(x)).$ For notational convenience we use $w=R_{T}v$. Then, we see that $\hat{w}=R\hat{v}$. Using a change of variables formula we get

[TABLE]

We have $|B_{ij}|\leq C\,h_{T}$, $|B_{ij}^{-1}|\leq C\,h_{T}^{-1}.$ Therefore, we obtain using (36)

[TABLE]

It is standard to show, again using a change of variable formula, and the bounds for $B$ and $B^{-1}$ above that

[TABLE]

Therefore, we have shown

[TABLE]

The required bounds for $\|R_{T}v\|_{L^{2}(\mathbb{R}^{3})}$ and $\|\nabla R_{T}v\|_{L^{2}(\mathbb{R}^{3})}$ follow a similar argument; see Lemma 5 in [23]. ∎

We are now ready to prove the FE trace inequality (26). Let $T\in\mathcal{T}_{h}$ and let $v\in H^{1}(T)$. Then, we apply (34) first for $\omega=\Gamma_{i}$ and next for $\omega=\Omega$, to get (we again use $w=R_{T}v$ for notation convenience):

[TABLE]

We apply Young’s inequality and use $h_{T}\leq h_{0}$ to handle terms on the right hand side. For example, we estimate

[TABLE]

Other terms are treated in the same way to get

[TABLE]

The desired result in (26) now follows after applying (35).