#  An extension of the box method discrete fracture model (Box-DFM) to include
low-permeable barriers with minimal additional degrees of freedom

Ziyao Xu111Department of Applied and Computational Mathematics and Statistics,
University of Notre Dame, Notre Dame, IN 46556, USA. E-mail<EMAIL_ADDRESS>Dennis Gläser222Institute for Modelling Hydraulic and Environmental Systems,
University of Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart, Germany. E-mail:
<EMAIL_ADDRESS>
Abstract. The box method discrete fracture model (Box-DFM) is  an important
finite volume-based discrete fracture model (DFM) to simulate flows in
fractured porous media. In this paper, we investigate a simple but effective
extension of the box method discrete fracture model to include low-permeable
barriers. The method remains identical to the traditional Box-DFM [41, 48] in
the absence of barriers. The inclusion of barriers requires only minimal
additional degrees of freedom to accommodate pressure discontinuities and
necessitates minor modifications to the original coding framework of the Box-
DFM. We use extensive numerical tests on published benchmark problems and
comparison with existing finite volume DFMs to demonstrate the validity and
performance of the method.

Key words. box method, discrete fracture model, fractured porous media,
extension, low-permeable barriers

##  1 Introduction

As a result of geological processes and human activities, cracks are widely
distributed in subsurface rocks. Depending on the degree of precipitation,
cracks can act either as highly permeable, flow-preferential channels or as
low-permeability, flow-blocking barriers; both have a significant impact on
the hydraulic properties of rocks. Understanding fluid flow in fractured rocks
is vital for many engineering fields, such as hydrocarbon reservoir
engineering, hydraulic fracturing, subsurface contaminant transport and
geothermal energy extraction. Hence, this area has attracted great interest
from researchers.

Due to the extreme contrast in permeability and lateral scale between
fractures and their surrounding rocks, accurately and efficiently modeling
fluid flow in fractured porous media is highly nontrivial. In recent decades,
many effective approaches have been developed to address this challenge. These
approaches can be divided into two categories: continuum models and discrete
fracture models. In continuum models, fractures are either homogenized into
the surrounding matrix, as seen in the effective permeability approach [47,
17], or integrated into a separate continuum that exchanges flux with the
matrix, as in the dual-porosity dual-permeability models [9, 54] and multi-
continuum models [56, 37]. Continuum models are preferable when fractures are
small, dense, and fully connected within rocks. However, model error becomes
significant when large, individual fractures dominate the flow fields. As a
compromise, discrete fracture models (DFMs) explicitly account for the effects
of individual fractures as lower-dimensional geometries, thus providing a more
accurate description of flows in fractured media.

Various DFMs have been developed based on different numerical methods. Early
works on DFMs [46, 7] were established based on the Galerkin finite element
methods (FEM-DFM) on  conforming grids where fractures are aligned with the
faces of grid cells. In these models, the  stiffness matrices of flows in
fractures are superposed on that of the porous matrix (see also [34, 32]).
Another approach for simulating subsurface flows is the class of control
volume methods. Based on the cell-centered control volume methods with two-
point flux approximation (TPFA), the TPFA-DFM was developed [33]. This method
treats fracture elements as geometrically low-dimensional but conceptually
equi-dimensional control volumes and computes their flux exchanges with other
control volumes based on transmissibility coefficients. It is further extended
to MPFA-DFM [49, 4, 23], based on the multi-point flux approximation (MPFA)
[1, 2], to allow for full tensor permeabilities on grids that are not
K-orthogonal. The Box methods [8, 26, 28], where the control volumes are
vertex-centered boxes and the function space is the  linear Lagrange finite
element, represent another widely used class of control volume methods. Based
on the Box methods, the Box-DFM has been developed [41, 48], where the grids
are conforming, and pressure continuity across fractures is assumed, as also
studied in [42, 43]. Models based on other numerical methods on conforming
meshes are also actively studied; for example, the mixed finite element
methods (MFEM) [39, 29, 30, 60], discontinuous/enriched Galerkin methods
(DG/EG) [6, 45, 15, 31], vertex approximate gradient methods [13, 14], and
mortar approaches [21, 12]. To alleviate the grid constraints when fracture
networks are complex, the non-conforming discrete fracture models are also
studied in the literature, e.g., the embedded discrete fracture model (EDFM)
[38, 40, 22, 52], Lagrange multiplier methods [36, 50], XFEM-DFM [18, 51, 20],
reinterpreted discrete fracture model [57, 58, 32], and immersed finite
element DFM [59].

Among the aforementioned discrete fracture models, control volume DFMs are
especially popular due to their local mass conservation, flexibility for
complex domains, and ease of implementation. Cell-centered control volume DFMs
have advantages in coupling fracture elements of arbitrary permeability with
matrix elements, as additional degrees of freedom (DoFs) are assigned to
fracture elements. However, cell-centered schemes require more DoFs on simplex
grids than vertex-centered schemes. The number of cell elements is roughly
twice that of vertices on triangular meshes, and the ratio is  approximately
six on tetrahedral meshes. Moreover, MPFA methods typically produce more
nonzero entries in the resulting stiffness matrix than box methods [24]. In
[27], a comparison between MPFA-DFM and Box-DFM was presented. The Box method
is capable of handling full tensor permeability on unstructured meshes as well
as MPFA methods, but with fewer DoFs. Furthermore, the stiffness matrix of
Box-DFM is symmetric positive-definite, which allows for the use of many fast
linear solvers. However, the Box-DFM can only be used for high-permeable
fractures, as discussed in [27]. Extended Box-DFM to account for wider type of
fractures in a uniform way has been developed in [24], where the equations of
fracture flows are established and coupled with the matrix equations based on
the interface model established in [39]. A treatment of assigning individual
equations for fracture variables was also developed in [11]. Both methods are
not identical to the Box-DFM [41, 48] when only high-permeable fractures
present.

In this paper, we propose a simple yet effective extension of the Box-DFM to
include low-permeable barriers. The model ignores the flow in the tangential
direction along barriers and adopts a broken Lagrange finite element space as
the trial functions to accommodate pressure discontinuity. The low-dimensional
Darcy’s law is adopted to compute the flux across cell interfaces aligned with
low-permeable barriers. Since no individual equations and variables are
assigned to fractures and barriers, the inclusion of barriers requires only
minimal additional degrees of freedom, and the model is identical to the
original Box-DFM in the absence of barriers. Moreover, we prove that the
resulting stiffness matrix remains symmetric positive-definite, as in the
original Box-DFM. In the numerical section, we test the model on published
benchmarks to demonstrate that such a simple treatment for low-permeable
barriers achieves high accuracy. The range of validity is also studied in this
section.

The rest of the paper is organized as follows. In Section 2, we describe the
model problem of steady-state single-phase flow through a fractured porous
medium. In Section 3, we first provide a review of the classic Box-DFM.
Subsequently, we introduce the algorithm of the extended Box-DFM to include
barriers and explore its properties. In Section 4, the validity and
performance of the extended Box-DFM are tested by published benchmarks before
a conclusion and outlook in Section 5. The source code to all examples
implemented in DuMux can be found at https://git.iws.uni-stuttgart.de/dumux-
pub/xu2024.

##  2 Problem description

We consider the steady-state single-phase flow through a porous medium
governed by the Darcy’s law.

In  $n$ dimensional space, $n\in\\{2,3\\}$, the equi-dimensional model is
formulated as

$-\nabla\cdot\left(\mathbf{K}\nabla
p\right)=q,\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{in}}~{}\Omega,$
(2.1)

subject to the boundary conditions

$p=g_{D}~{}\text{on}~{}\Gamma_{D},\quad-(\mathbf{K}\nabla
p)\cdot\mathbf{n}=g_{N}~{}\text{on}~{}\Gamma_{N}:=\partial\Omega\setminus\Gamma_{D},$
(2.2)

where $\nabla$ is the gradient operator, $\Omega\subset\mathbb{R}^{n}$ is an
open bounded domain, $\mathbf{K}$ is the permeability tensor of porous media,
$p$ is the pressure of fluid, $q$ is the source term, $\Gamma_{D}$ and
$\Gamma_{N}$ are the Dirichlet and Neumann boundaries, respectively, and
$\mathbf{n}$ is the unit outer normal on $\partial\Omega$.

In the absence of fractures and barriers, the problem (2.1), (2.2) is
generally an inhomogeneous anisotropic elliptic equation, which can be solved
by traditional numerical methods on regular meshes, e.g., the Box-methods [26,
28] or the cell-centered finite volume methods [1, 2]. However, in the
presence of fractures and barriers, numerical methods based on the equi-
dimensional setting may not be well-suited due to the extreme contrast in
geometry and permeability between fractures/barriers compared to the porous
matrix.

$\mathbf{n}^{-}$$\mathbf{n}^{+}$$\gamma$${\Omega}_{-}$${\Omega}_{+}$ Figure
2.1: The geometry of an interface model  in $\mathbb{R}^{2}$. $\gamma$
represents the interface of a fracture or barrier. $\Omega^{\pm}$ and
$\mathbf{n}^{\pm}$ are the bulk matrix regions and unit outer normal from each
region, respectively.

The hybrid-dimensional model, which treats the fractures and barriers as co-
dimension one manifolds, provides an efficient way to explicitly take into
account the effects of fractures and barriers. Consider a fractured porous
medium $\Omega$ subdivided into $\Omega^{\pm}$ by the $(n-1)$-dimensional
fracture/barrier $\gamma$. For an illustration, see Figure 2.1. In this paper,
we consider the governing equation

$-\nabla\cdot\left(\mathbf{K}_{m}^{\pm}\nabla
p^{\pm}\right)=q^{\pm},\quad\text{in}~{}\Omega^{\pm},$ (2.3)

in the porous matrix, and the interface condition

$\begin{cases}p_{f}=p^{\pm},\\\
u_{f}=-k_{f}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla_{\tau}p_{f}},\\\
a{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla_{\tau}\cdot
u_{f}}=-[\mathbf{K}_{m}\nabla
p\cdot\mathbf{n}],\end{cases}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{on}~{}}\gamma,$
(2.4)

if $\gamma$ is a high-permeable fracture, or

$-\mathbf{K}^{-}_{m}\nabla p^{-}\cdot\mathbf{n}^{-}=\mathbf{K}^{+}_{m}\nabla
p^{+}\cdot\mathbf{n}^{+}=-k_{b}\frac{p^{+}-p^{-}}{a},\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{on}~{}}\gamma,$
(2.5)

if $\gamma$ is a low-permeable barrier. Here, we denote by
$\mathbf{K}_{m}^{\pm}$, $p^{\pm}$ and $q^{\pm}$ the permeability, pressure and
source on $\Omega^{\pm}$, respectively. Moreover, $k_{f}$ is the permeability
of high-permeable fractures, $k_{b}$ is the permeability of low-permeable
barriers, $\nabla_{\tau}$ is the gradient operator on the tangent space of
$\gamma$, $[\mathbf{K}_{m}\nabla p\cdot\mathbf{n}]=\mathbf{K}_{m}^{-}\nabla
p^{-}\cdot\mathbf{n}^{-}+\mathbf{K}_{m}^{+}\nabla p^{+}\cdot\mathbf{n}^{+}$ is
the jump of normal flux across $\gamma$, $\mathbf{n}^{\pm}$ is the unit outer
normal on $\gamma$ from $\Omega^{\pm}$, and $a$ is the aperture of $\gamma$.

######  Remark 2.1.

The hybrid-dimensional model (2.3), (2.4), (2.5) considered in this paper is
an important limit case of the general interface model established in [39].
The pressure continuity in the fractures is assumed, and the tangential flow
in the barriers is neglected here. See discussions about the equations (4.5),
(4.6) therein.

We would like to emphasize that we do not necessarily assume the
fracture/barrier $\gamma$ split the whole domain $\Omega$ into two separate
pieces. The geometric setting above is only for the ease of illustration. We
shall drop the superscripts $\pm$ in the porous matrix equation (2.3), and use
$\pm$ only for the interface conditions in the next section.

##  3 Algorithm

In this section, we review the algorithm of the traditional Box-DFM and
propose the inclusion of low-permeable barriers. The illustration is based on
two-dimensional space with triangular grids for convenience, but the extension
to three-dimensional space or quadrilateral grids is straightforward.

###  3.1 A review of the traditional Box-DFM

We denote the computational domain by $\Omega_{h}$ and consider a Delaunay
triangulation $\\{T:T\in\mathcal{T}\\}$ of the domain, where $T$ are non-
overlapping triangular cells such that $\cup_{T\in\mathcal{T}}T=\Omega_{h}$.
The triangular mesh $\\{T:T\in\mathcal{T}\\}$ is called the primal mesh. The
corresponding dual mesh $\\{B:B\in\mathcal{B}\\}$ is obtained through
connecting the midpoints of the sides of $T$ to its centroid, see Figure 1(a)
for illustration, where the boxes $B\in\mathcal{B}$ in the dual mesh are
centered at vertices of the primal mesh.

We define the $P^{1}$-Lagrange finite element space on the mesh
$\\{T:T\in\mathcal{T}\\}$ as

$V_{h}=\\{v\in C(\Omega_{h}):v|_{T}\in P^{1}(T),T\in\mathcal{T}\\},$ (3.1)

where $P^{1}(T)$ is the linear function space on the cell $T$ and
$C(\Omega_{h})$ is the space of continuous functions on $\Omega_{h}$.

We first briefly review the box method for the equi-dimensional model (2.1).
We look for the numerical solution
$p=\sum_{i=1}^{N}p_{i}\psi_{i}(\mathbf{x})\in V_{h}$, such that

$-\int_{\partial B}\mathbf{K}\nabla
p\cdot\mathbf{n}ds=\int_{B}q\,dxdy,\quad\forall B\in\mathcal{B},$ (3.2)

where $N$ is the number of boxes in the dual mesh (also the number of vertices
in the primal mesh), $p_{i}$ is the nodal value of the solution at the vertex
$\mathbf{x}_{i}$, and $\psi_{i}(\mathbf{x})$ is the shape functions serving as
the basis of the $P^{1}$-Lagrange finite element space at $\mathbf{x}_{i}$.
The scheme (3.2) is obtained by integrating (2.1) over every box
$B\in\mathcal{B}$ and applying the Gauss divergence theorem. For brevity, we
omit the discussion of the boundary conditions (2.2) here and henceforth.

We then adopt the traditional Box-DFM as described in [41, 48] to explicitly
take into account the flow contributed by the high-permeable fracture. The
fractures are aligned with the edges of the primal mesh and intersect the dual
mesh at midpoints of the sides of primal cells, as shown in Figure 1(b). Under
the assumption of pressure continuity and through the superposition of
fracture flows, the Box-DFM for the hybrid-dimensional model (2.3), (2.4) is
to find the numerical solution $p\in V_{h}$, such that

$-\int_{\partial B}\mathbf{K}_{m}\nabla
p\cdot\mathbf{n}ds-\sum_{j=1}^{k}ak_{f}\frac{\partial
p}{\partial\nu_{j}}(\mathbf{x}^{F}_{j})=\int_{B}q\,dxdy,\quad\forall
B\in\mathcal{B},$ (3.3)

where $\mathbf{x}^{F}_{j},j=1,2,\ldots,k$ are points where the fractures
intersect $\partial B$, $a$ and $k_{f}$ are aperture and permeability of the
fracture, and $\nu_{j}$ is the unit tangential vector of the fracture at
$\mathbf{x}^{F}_{j}$ pointing from inside of $B$ to its outside.

######  Remark 3.1.

There are multiple ways to derive the scheme (3.3). For example, one can
directly plug the permeability tensor $\mathbf{K}$, established in [57], or
the effective permeability tensor derived in [47], into the model (3.2) to
attain (3.3).  The scheme can also be directly derived from the interface
model (2.3) and (2.4), by establishing separate equations for the sub-boxes
divided by the fractures, summing up the equations belonging to the same box,
substituting the interface condition, and applying the fundamental theorem of
calculus on the fracture edge.

1234567 (a) Unfractured porous media. Solid lines denote the primal mesh,
dashed lines denote a box centered at the vertex $1$ in the dual mesh.

1234567$\mathbf{x}^{F}_{1}$$\mathbf{x}^{F}_{2}$ (b) Porous media with high-
permeable fractures. Blue lines denote fractures, which intersect the box at
$\mathbf{x}_{1}^{F}$ and $\mathbf{x}_{2}^{F}$.

3567$1^{\prime}$$1^{\prime\prime}$$1^{*}$$2^{\prime}$$2^{\prime\prime}$$4^{\prime}$$4^{\prime\prime}$
(c) Porous media with low-permeable barriers. Red lines denote barriers, which
subdivide the box into three sub-boxes. Degrees of freedom of vertices on
barriers are duplicated accordingly to account for pressure discontinuity.

Figure 3.1: The meshes used in box methods for unfractured porous media,
porous media with high-permeable fractures and low-permeable barriers.

###  3.2 The inclusion of low-permeable barriers

Similar to the case with fractures, we assume that the barriers align with the
edges of cells in the primal mesh.

We define the extended $P^{1}$-Lagrange finite element space on the fitting
mesh $\\{T:T\in\mathcal{T}\\}$ as

$\overline{V}_{h}=\\{v\in C(\Omega_{h}\setminus\Gamma):v|_{T}\in
P^{1}(T),T\in\mathcal{T}\\},$ (3.4)

where $\Gamma=\cup_{i}\gamma_{i}$ is the union of all barriers and
$C(\Omega_{h}\setminus\Gamma)$ is the space of continuous functions on
$\Omega_{h}\setminus\Gamma$, namely the functions admit discontinuities only
on $\Gamma$. Compared with the finite element space $V_{h}$ defined in (3.1),
the extended space $\overline{V}_{h}$ in (3.4) is endowed with additional
degrees of freedom at the vertices on barriers to accommodate the
discontinuity of pressure

In the dual mesh, the boxes $B\in\mathcal{B}$ may be subdivided into separate
sub-boxes by the barriers $\Gamma$, as illustrated in Figure 1(c).
Consequently, we define the extended collection of boxes
$\overline{\mathcal{B}}$ as encompassing both the original boxes
$B\in\mathcal{B}$ that remain undivided by $\Gamma$, and the sub-boxes
resulting from such subdivision. As a result, the number of boxes in
$\overline{\mathcal{B}}$ is the same as the degrees of freedom of
$\overline{V}_{h}$. It is worth making a special note about the situation of
the tips of $\Gamma$ to eliminate ambiguity. According to the definition of
$\overline{V}_{h}$, the vertices at the tips of $\Gamma$ have only one degree
of freedom. Correspondingly, the boxes $B\in\mathcal{B}$ that enclose these
vertices are regarded as undivided.

The extended Box-DFM for the hybrid-dimensional model (2.3), (2.5) is to find
the numerical solution $p\in\overline{V}_{h}$, such that

$-\int_{\partial B\setminus\Gamma}\mathbf{K}_{m}\nabla
p\cdot\mathbf{n}ds-\int_{\partial
B\cap\Gamma}k_{b}\frac{p^{+}-p^{-}}{a}ds=\int_{B}q\,dxdy,\quad\forall
B\in\overline{\mathcal{B}},$ (3.5)

where $p^{-}$ and $p^{+}$ are pressure on $\partial B\cap\Gamma$ from the
inside and outside of $B$, respectively. For those undivided boxes
$B\in\mathcal{\overline{B}}\cap\mathcal{B}$, it’s clear that $\partial
B\setminus\Gamma=B$ and $\partial B\cap\Gamma=\emptyset$, thus the equation is
identical to (3.2).

######  Remark 3.2.

People have long noticed [8] that the stiffness matrix of the box method (3.2)
is identical to the stiffness matrix of the linear finite element method; see
a concise proof in [26]. Following the same approach therein, one can prove
that the stiffness matrix of Box-DFM (3.3) is identical to that of FEM-DFM
[32]. Therefore, the stiffness matrices of the box method (3.2) and Box-DFM
(3.3) are symmetric positive-definite, as in the finite element method, which
makes them suitable for many fast solvers. Fortunately, the stiffness matrix
of the extended Box-DFM (3.5) is still symmetric positive-definite. A proof is
attached in Appendix A.

Finally, we combine the schemes (3.3) and (3.5) to formulate our extended Box-
DFM for fracture and barrier networks: Find $p\in\overline{V}_{h}$, such that

$-\int_{\partial B\setminus\Gamma}\mathbf{K}_{m}\nabla
p\cdot\mathbf{n}ds-\sum_{j=1}^{k}ak_{f}\frac{\partial
p}{\partial\nu_{j}}(\mathbf{x}^{F}_{j})-\int_{\partial
B\cap\Gamma}k_{b}\frac{p^{+}-p^{-}}{a}ds=\int_{B}q\,dxdy,\quad\forall
B\in\overline{\mathcal{B}},$ (3.6)

where $\mathbf{x}^{F}_{j},j=1,2,\ldots,k$ are points where the fractures
intersect $\partial B$.

######  Remark 3.3.

If the fracture and barrier intersect at a vertex, the treatment there depends
on whether the fracture penetrates the barrier, or the barrier cuts the
fracture. In the former case, there is only one degree of freedom, one box,
and one equation associated with this vertex. In the latter case, degrees of
freedom are assigned to each sub-box cut by the barriers, and equations are
established for each sub-box.  For an illustration of the differences between
these two approaches, we refer the reader to Example 4 in the numerical
section, where both treatments are implemented and their outcomes
demonstrated.

##  4 Numerical experiments

In this section, we test the effectiveness and range of validity of our
extended Box-DFM  on a set of test cases including published benchmark
problems, see e.g., [19, 5, 59, 24].  If not otherwise stated, the units used
in the examples adhere to SI, i.e., $\mathbf{K}$ in $$\mathrm{m}$^{2}$, $p$ in
$\mathrm{P}\mathrm{a}$, $a$ in $\mathrm{m}$, etc. The results are compared
with those obtained from another extension of Box-DFM developed in [24], which
established separate equations and assigned individual degrees of freedom for
arbitrary types of fractures. To distinguish between the two methods in the
presentation, we refer to our method as ’box-dfm’ and the method developed in
[24] as ’ebox-dfm’, following the naming convention the authors already
established in the literature.

The problems are presented in roughly increasing order of geometrical
complexity. For simplicity, we set the source term $q=0$ in all examples
except the convergence test. Since the method is identical to the traditional
Box-DFM in the absence of barriers, we focus only on benchmarks involving low-
permeable barriers. In all the tests, except for the last one, we assume equal
normal and tangential permeabilities in the barriers.

######  Example 1.

Convergence test

In this experiment, we test the convergence of the method (3.5) to the model
problem (2.3), (2.5) by considering the following setting:
$\mathbf{K}_{m}=\mathbf{I}$, $\frac{k_{b}}{a}=1$,
$\Omega^{-}=[0,\frac{1}{2}]\times[0,1]$,
$\Omega^{+}=[\frac{1}{2},1]\times[0,1]$, $\gamma:\\{x=\frac{1}{2},0\leq y\leq
1\\}$, and

$\begin{cases}p^{-}=\sin(x)\sin(y),~{}q^{-}=2\sin(x)\sin(y),&(x,y)\in\Omega^{-},\\\
p^{+}=p^{-}+\cos(\frac{1}{2})\sin(y),~{}q^{+}=q^{-}+\cos(\frac{1}{2})\sin(y),&(x,y)\in\Omega^{+},\end{cases}$
(4.1)

with the corresponding Dirichlet boundary condition.

We test the algorithm on unstructured grids with different level of
refinement. The grid used for the lowest level of refinement, consisting of
256 triangles, is shown in Figure 4.1. The $L^{2}$-error in porous matrix and
the orders of convergence are presented in Table 1, where $p_{h}$ denotes the
numerical solution and $p$ is the exact solution given in (4.1). From the
table, we can clearly observe the  expected order of convergence of the box
method.

Figure 4.1: Example 1: convergence test. The grid used for the lowest refinement, consisting of $256$ triangles. The red line represents the barrier interface. $i$ | $||p_{h}-p||_{2}$ | Order
---|---|---
0 | 5.39e-04 | -
1 | 1.36e-04 | 1.99
2 | 3.40e-05 | 2.00
3 | 8.50e-06 | 2.00
4 | 2.13e-06 | 2.00
5 | 5.31e-07 | 2.00
Table 1: Example 1: convergence test. $L^{2}$ errors evaluated on
$\Omega^{\pm}$ on different levels of refinement

######  Example 2.

Single barrier

In this experiment, we examine the performance of our method in cases
involving a single barrier. We consider two scenarios: the first involves a
vertical barrier from (0.5, 0.5) to (0.5, 1), see the same test in references
[5, 59]; the second tests a slanted barrier from (0.25, 0.75) to (0.75, 0.25).
In both scenarios, the computational domain is defined as $\Omega=[0,1]^{2}$,
and the permeability tensor in the porous matrix is
$\mathbf{K}_{m}=\mathbf{I}$. The ratio of permeability to aperture of barriers
is given by $\frac{k_{b}}{a}=10^{-5}$. The top and bottom boundaries are
impermeable ($g_{N}=0$), while the left and right boundaries are set with
Dirichlet conditions, $g_{D}=0$ and $g_{D}=1$, respectively.

We conduct the computation using our extended box-dfm and the ebox-dfm
developed in [24] on identical grids. The pressure contours obtained from our
method, alongside the background mesh, are presented in Figure 4.2 (a), (b).
The differences in pressure between the two methods are illustrated in Figure
4.2 (c), (d). The comparison reveals that the discrepancies in the results of
the two methods are tiny. Additionally, we present a comparison of the
pressure profiles obtained by the two methods along specific slices in Figure
4.3. This comparison also exhibits an excellent agreement between two methods.

(a) vertical barrier

(b) slanted barrier

(c) vertical barrier

(d) slanted barrier

Figure 4.2: Example 2: single barrier. Visualization of the solution $p$ (row
$1$) obtained with the proposed method, and the difference $\Delta p$ (row
$2$) to the solution obtained with the ebox-dfm in [24]. In the first scenario
(vertical barrier), the grid contains $253$ vertices and $450$ triangles; in
the second scenario (slanted barrier), it contains $229$ vertices and $404$
triangles.

(a) vertical barrier

(b) slanted barrier

Figure 4.3: Example 2: single barrier. Pressure profiles along the slice (0,
0.75) – (1, 0.75) for the vertical barrier and along the slice (0, 0.5) – (1,
0.5) for the slanted barrier.

######  Example 3.

Regular barrier network

In this experiment, we test the performance of our method on a regular barrier
network. The computational domain is defined as $\Omega=[0,1]^{2}$, and the
permeability tensor of the porous matrix is given by
$\mathbf{K}_{m}=\mathbf{I}$. The domain encompasses six barriers, with
coordinates detailed in [19]. Each barrier has a uniform aperture $a=10^{-4}$
and a permeability $k_{b}=10^{-4}$. The top and bottom boundaries are
impermeable; the left boundary is set as a Neumann boundary with $g_{N}=1$;
and the right boundary is defined by a Dirichlet condition with $g_{D}=1$.

We conduct computations using both our extended box-dfm and the ebox-dfm on
three different grids, with a varying level of refinement. The pressure
contours produced by our method, alongside the background mesh, and the
differences in pressure between the two methods, are illustrated in Figure
4.4. Furthermore, we present pressure slices along the line from (0, 0.1) to
(0.9, 1) in Figure 4.5, which are compared with the reference solution [19]
obtained using the mimetic finite difference (MFD) method with $1,136,456$
elements for the equi-dimensional model.

Unlike the scenarios involving a single barrier, noticeable discrepancies in
the solutions, especially on the coarse and medium refinement grids, between
the two methods are observed. This is possibly due to their difference in
handling the intersection of low-permeable barriers. In our method, tangential
flows on barriers are completely neglected, thus no special condition at
intersections is imposed.  The neglected intersection volume in the ebox-dfm
model may lead to an overestimation of the transfer fluxes across the barrier
intersections in comparison to the equi-dimensional reference, in particular
on coarse meshes. See the subsequent example for more investigations on the
influence of differences in intersection handling.  From the comparison of
slices in Figure 4.5, our results show better agreement with the reference
solution.

(a) coarse grid

(b) medium grid

(c) fine grid

(d) coarse grid

(e) medium grid

(f) fine grid

Figure 4.4: Example 3: regular barrier network. Visualization of the solution
$p$ (row $1$) obtained with the proposed method, and the difference $\Delta p$
(row $2$) to the solution obtained with the ebox-dfm in [24]. The coarse grid
contains $29$ vertices and $46$ triangles; the medium grid contains $205$
vertices and $366$ triangles; the fine grid contains $12,112$ vertices and
$23,822$ triangles.

(a) coarse grid

(b) medium grid

(c) fine grid

Figure 4.5: Example 3: regular barrier network. Pressure profiles along the
slice (0, 0.1) – (0.9, 1). The coarse grid contains $29$ vertices and $46$
triangles; the medium grid contains $205$ vertices and $366$ triangles; the
fine grid contains $12,112$ vertices and $23,822$ triangles. The reference,
provided by the authors in [19], is obtained from the mimetic finite
difference (MFD) method with $1,136,456$ elements for the equi-dimensional
model.

######  Example 4.

Complex fracture-barrier network

In this experiment, we test a complex fracture-barrier network as benchmark
4.3 in [19], where the high-permeable fractures and low-permeable barriers
coexist and intersect.

The computational domain is $\Omega=[0,1]^{2}$, and the permeability tensor of
the porous matrix is $\mathbf{K}_{m}=\mathbf{I}$. The domain contains 8 high-
permeable fractures and 2 low-permeable barriers, with a uniform aperture
$a=10^{-4}$. All fractures have a uniform permeability $k_{f}=10^{4}$, and all
barriers have a uniform permeability $k_{b}=10^{-4}$. For the exact
coordinates of the fractures and barriers, we refer the readers to the
Appendix C in [19]. Two sub-cases are considered. In the sub-case (a), where
the pressure gradient is predominantly vertical, the left and right boundaries
are impermeable, i.e., $g_{N}=0$; and the top and bottom boundaries are
Dirichlet with $g_{D}=4$ and $g_{D}=1$, respectively. In the sub-case (b),
where the pressure gradient is predominantly horizontal, the top and bottom
boundaries are impermeable, i.e., $g_{N}=0$; and the left and right boundaries
are Dirichlet with $g_{D}=4$ and $g_{D}=1$, respectively.

The treatment of intersections between high-permeability fractures and low-
permeability barriers is a non-trivial issue. In [19], a harmonic average of
the permeabilities of the fracture and barrier for the intersection cell is
adopted in methods that explicitly account for the intersections. In the ebox-
dfm, pressure and flux continuity is enforced at intersections. In our
proposed method, we prioritize the barrier at intersections, imposing
discontinuous pressure at the intersection, but we will discuss an alternative
approach in the sequel.

We conduct the computation on a grid containing $1,391$ vertices and $2,680$
triangles using our extended box-dfm and the ebox-dfm for both cases. The
pressure contours produced by our method, alongside the background mesh, and
the differences in pressure between the extended box-dfm and the ebox-dfm, are
presented in Figure 4.6. Additionally, we present pressure slices along the
line from (0, 0.5) to (1, 0.9) in Figure 4.7, comparing them with the
reference solution provided in [19].

In the reference solution, a harmonic average of the fracture and barrier
permeabilities were assigned to the cells occupying the intersection region,
thereby the intersection exhibits a more significant blocking feature. In
reality, it is also possible for the intersection to show a more flow-
preferable conducting feature. Therefore, we attempt an alternative approach
to impose pressure continuity at the intersection and present the
corresponding pressure slices in Figure 4.8. In these figures, we observe a
larger deviation of our results from the reference solutions, but a smaller
deviation from the results of the ebox-dfm is observed. These findings are in
line with the discussions in [19] on the influence of not modelling the
intersection region explicitly.

(a) Flow from top to bottom

(b) Flow from left to right

(c) Flow from top to bottom

(d) Flow from left to right

Figure 4.6: Example 4: complex fracture-barrier network. Visualization of the
solution $p$ (row $1$) obtained with the proposed method, and the difference
$\Delta p$ (row $2$) to the solution obtained with the ebox-dfm in [24]. The
first scenario is the flow from top to bottom, and the second scenario is the
flow from left to right. Both cases contain $1,391$ vertices and $2,680$
triangles. Pressure discontinuity is assumed at the intersections of fracture
and barrier in the box-dfm.

(a) Flow from top to bottom

(b) Flow from left to right

Figure 4.7: Example 4: complex fracture-barrier network. Pressure profiles
along the slice (0, 0.5) – (1, 0.9). The reference, provided by the authors in
[19], is obtained from the mimetic finite difference (MFD) method with
$2,260,352$ elements for the equi-dimensional model. Pressure discontinuity is
assumed at the intersections of fracture and barrier in the box-dfm.

(a) Flow from top to bottom

(b) Flow from left to right

Figure 4.8: Example 4: complex fracture-barrier network. Pressure profiles
along the slice (0, 0.5) – (1, 0.9). Pressure continuity is assumed at the
intersections of fracture and barrier in the box-dfm.

######  Example 5.

Realistic barrier network

In this example, we test the case of a realistic barrier network containing 64
low-permeable barriers whose distribution is based on an image of an outcrop
near Bergen, Norway. This test case is taken from [24], where the authors
transformed high-permeable fractures from the original test case presented in
[19] into a variant featuring low-permeable barriers to validate various
discrete fracture models.

The computational domain is
$\Omega=[0,700$\mathrm{m}$]\times[0,600$\mathrm{m}$]$, and a flow from left to
right is driven by the Dirichlet boundary condition
$g_{D}=1,013,250$\mathrm{m}$$ at $x=0$\mathrm{m}$$ and $g_{D}=0$\mathrm{m}$$
at $x=700$\mathrm{m}$$, along with impermeable boundary conditions at
$y=0$\mathrm{m}$$ and $y=600$\mathrm{m}$$. All barriers in the domain have a
uniform aperture of $a=10^{-2}$\mathrm{m}$$ and permeability
$k_{b}=10^{-18}$\mathrm{m}\mathrm{/}\mathrm{s}$$. The permeability of the
porous matrix is set to $k_{m}=10^{-14}$\mathrm{m}\mathrm{/}\mathrm{s}$$. For
detailed coordinates of the barriers, we refer the readers to [19].

We conduct the computation on three grids with different levels of refinement
and compare the results of our extended box-dfm with those obtained from the
ebox-dfm. The results of tpfa-dfm, mpfa-dfm, and ebox-dfm-mortar methods for
the same test case are available in [24] for a more comprehensive comparison.
From the pressure contours and slices shown in Figure 4.9 and Figure 4.10,
respectively, one can observe the convergence between the results of the box-
dfm and ebox-dfm with mesh refinement. The pressure profiles along (625, 0) –
(625, 600) show that our method’s results on a coarse mesh exhibit better
agreement with the reference solution, which is obtained using the ebox-dfm on
a fine mesh consisting of $307,829$ cells, while the finest mesh of the
refinement series consists of $114,721$ cells..

(a) coarse grid

(b) medium grid

(c) fine grid

(d) coarse grid

(e) medium grid

(f) fine grid

Figure 4.9: Example 5: realistic barrier network. Visualization of the
solution $p$ (row $1$) obtained from the proposed method, and the difference
$\Delta p$ (row $2$) to the solution obtained with the ebox-dfm in [24]. The
coarse grid contains $1,884$ vertices and $3,611$ triangles; the medium grid
contains $4,968$ vertices and $9,646$ triangles; the fine grid contains
$57,904$ vertices and $114,721$ triangles.

(a) coarse grid

(b) medium grid

(c) fine grid

(d) coarse grid

(e) medium grid

(f) fine grid

Figure 4.10: Example 5: realistic barrier network. Pressure profiles along the
slice (0,0) – (700, 600) (row $1$) and along the slice (625, 0) – (625, 600)
(row $2$). The coarse grid contains $1,884$ vertices and $3,611$ triangles;
the medium grid contains $4,968$ vertices and $9,646$ triangles; the fine grid
contains $57,904$ vertices and $114,721$ triangles. The reference is obtained
from the ebox-dfm with $154,814$ vertices and $307,829$ triangles.

######  Example 6.

Three-dimensional test case

In this example, we test the flow in the unit cube $\Omega=[0,1]^{3}$,
containing $9$ regularly oriented barriers. One can find the exact coordinates
of the barriers in case $2$ presented in [10]. The boundary $\partial\Omega$
is divided into three parts. The first part,
$\partial\Omega_{D}=\\{(x,y,z)\in\partial\Omega:x,y,z>0.875\\}$, is imposed
with the Dirichlet condition $g_{D}=1$; the second part,
$\partial\Omega_{N}=\\{(x,y,z)\in\partial\Omega:x,y,z<0.25\\}$, is imposed
with an inflow boundary condition $g_{N}=1$; the remaining boundaries are
impermeable. All barriers have a uniform aperture $a=10^{-4}$ and permeability
$k_{b}=10^{-4}$. The permeability of the porous matrix is different in two
bulk regions. In
$\Omega_{1}=\\{(x,y,z)\in\Omega:x>0.5,y<0.5\\}\cup\\{(x,y,z)\in\Omega:x>0.75,0.5<y<0.75,z>0.5\\}\cup\\{(x,y,z)\in\Omega:0.625<x<0.75,0.5<y<0.625,0.5<z<0.75\\}$,
we take $k_{m1}=0.1$; in $\Omega_{2}=\Omega\setminus\Omega_{1}$, we take
$k_{m2}=1$.

We present the pressure contours and slices in Figure 4.11 and Figure 4.12,
respectively. As before, the convergence of box-dfm and ebox-dfm can be
observed with grid refinement. It can be seen from the slices that the ebox-
dfm is closer to the reference solution on the coarse grid, while our proposed
method is closer to the reference on the fine grid,  where the reference
solution is obtained from the data of the original benchmark study in [10].

(a) coarse grid

(b) medium grid

(c) fine grid

(d) coarse grid

(e) medium grid

(f) fine grid

Figure 4.11: Example 6: three-dimensional test case. Visualization of the
solution $p$ (row $1$) obtained from the proposed method, and the difference
$\Delta p$ (row $2$) to the solution obtained from the ebox-dfm in [24]. The
coarse grid contains $948$ vertices and $3,840$ tetrahedrons; the medium grid
contains $8,316$ vertices and $41,414$ tetrahedrons; the fine grid contains
$104,343$ vertices and $593,347$ tetrahedrons.

(a) coarse grid

(b) medium grid

(c) fine grid

Figure 4.12: Example 6: three-dimensional test case. Pressure profiles along
the line (0, 0, 0) – (1, 1, 1). The coarse grid contains $948$ vertices and
$3,840$ tetrahedrons; the medium grid contains $8,316$ vertices and $41,414$
tetrahedrons; the fine grid contains $104,343$ vertices and $593,347$
tetrahedrons. The reference, provided by the authors in [10], is obtained from
the control volume multi-point flux approximation (MPFA) method with
$1,046,566$ elements for the hybrid-dimensional model [35].

######  Example 7.

Range of validity

So far, we have tested only cases of barriers with low permeability in both
the normal and tangential directions. In this experiment, we investigate the
range of validity of the method for barriers with various tangential
permeabilities [5, 59].

We consider the computational domain $\Omega=[0,1]^{2}$ with matrix
permeability $\mathbf{K}_{m}=\mathbf{I}$. The domain contains four barriers,
$\gamma_{1},\gamma_{2},\gamma_{3}$, and $\gamma_{4}$, whose axes are given by:
(0.3, 0.2) – (1.0, 0.2), (0.0, 0.4) – (0.7, 0.4), (0.3, 0.6) – (1.0, 0.6), and
(0.0, 0.8) – (0.7 0.8), respectively. All barriers have a uniform aperture
$a=10^{-3}$ and normal permeability $k_{n}=10^{-3}$. We set the tangential
permeability of $\gamma_{2}$ and $\gamma_{4}$ equal to their normal direction,
i.e., $k_{\tau}=10^{-3}$. For the barriers $\gamma_{1}$ and $\gamma_{3}$,
their tangential permeabilities take various values: $k_{\tau}=10^{-3},1$ or
$10^{3}$.

Two sub-cases are tested as in [5, 59]. In sub-case (a), the left and right
boundaries are impermeable, i.e., $q_{N}=0$, and the top and bottom boundaries
are set as Dirichlet with $g_{D}=1$ and $g_{D}=0$, respectively. In sub-case
(b), the boundary condition is purely Dirichlet with
$g_{D}(x,y)=(2x-1)(3x-1)$.

We conduct the computation on a grid containing $11,957$ vertices and $23,512$
triangles using our extended box-dfm and ebox-dfm for all scenarios. The
pressure contours obtained from our proposed method for both cases are shown
in Figure 4.13. We present the pressure profiles along the line (0.65, 0) –
(0.65,1) in Figure 4.14, and compare them with the reference solution obtained
from the box method for equi-dimensional model on an extremely refined mesh.

As can be seen from the figures, the model errors are small if the tangential
permeability of the barriers is no greater than that of porous matrix,
provided the aperture of the barriers is small. However, for the scenarios
where barriers have high tangential permeability, the treatment of neglecting
tangential flow in the proposed model is not satisfactory. Therefore, this
extended box-dfm method should not be used for barriers that have high
tangential permeability,  a situation that is expected to be of little
relevance in practice.

(a) case (a)

(b) case (b)

Figure 4.13: Example 7: range of validity. Visualization of the solution $p$
obtained with the proposed method. Both cases contain $11,957$ vertices and
$23,512$ triangles. The tangential permeability $k_{\tau}=10^{-3}$ is adopted
in both cases.

(a) $k_{\tau}=10^{-3}$ in case (a)

(b) $k_{\tau}=1$ in case (a)

(c) $k_{\tau}=10^{3}$ in case (a)

(d) $k_{\tau}=10^{-3}$ in case (b)

(e) $k_{\tau}=1$ in case (b)

(f) $k_{\tau}=10^{3}$ in case (b)

Figure 4.14: Example 7: range of validity Pressure profiles along the line
(0.65, 0) – (0.65,1). The box-dfm and ebox-dfm are computed on a grid
containing $11,957$ vertices and $23,512$ triangles. The reference solution is
obtained from the box method with $269,152$ vertices and $538,164$ triangles
for the equi-dimensional model.

##  5 Summary

In this work, we have proposed a very simple extension of the Box-DFM to
include blocking barriers, by adopting a broken linear Lagrange finite element
space as the trial function space and using divided sub-boxes to apply the
Gauss divergence theorem. The proposed method is identical to the original
Box-DFM in the absence of blocking barriers and only requires minimal
additional degrees of freedom, compared with other variants of Box-DFMs, to
accommodate the pressure discontinuity across barriers. In addition, it
inherits the local mass conservation, flexibility in the computational grid,
and the symmetric positive-definiteness of the stiffness matrix from the
original method.

Several widely practiced benchmark problems are tested. Numerical experiments
exhibit good performance of our method when dealing with blocking barriers. It
should be noted that our method is based on the interface model, which
neglects the pressure jump across highly permeable fractures and the
tangential flow along blocking barriers. Thus, the method is not suitable for
handling cases of fractures with low permeability in the normal direction and
high permeability in the tangential direction. This assertion is verified in
the last numerical example. Moreover, it is not difficult to extend the
treatment for barriers in two-phase flow. In future work, we will study the
performance of such an extended Box-DFM in applications of multi-phase flows.

##  Appendix A The symmetry and positive definiteness of the stiffness matrix
of (3.5)

In this appendix, we show the symmetry and positive definiteness of the
stiffness matrix in (3.5), provided that $\mathbf{K}_{m}$ is symmetric
positive definite. For simplicity, we assume $\mathbf{K}_{m}$ is a constant
matrix and $\frac{k_{b}}{a}>0$ is a constant scalar. However, this condition
can be relaxed to allow for piecewise constant $\mathbf{K}_{m}$ on cells and
piecewise constant $\frac{k_{b}}{a}\geq 0$ on edges of the primal mesh,
provided that no area is enclosed by impermeable barriers.

######  Lemma A.1.

Let $\phi_{i}\in\overline{V}_{h}$ the Lagrange basis and
$B_{i}\in\overline{\mathcal{B}}$ the box associated with the node $i$, for
$i=1,\ldots,\\#\overline{\mathcal{B}}$. Then we have

$\sum_{T\in\mathcal{T}}\int_{T}\mathbf{K}_{m}\nabla
v\cdot\nabla\phi_{i}dxdy=-\int_{\partial
B_{i}\setminus\Gamma}\mathbf{K}_{m}\nabla v\cdot\mathbf{n}ds,\quad\forall
v\in\overline{V}_{h},$ (A.1)

where $\mathbf{n}$ denotes the unit outer normal of $\partial B_{i}$.

###### Proof.

It is proved in [26] that,

$\int_{T}\mathbf{K}_{m}\nabla v\cdot\nabla\phi_{i}dxdy=-\int_{\partial
B_{i}\cap T^{o}}\mathbf{K}_{m}\nabla v\cdot\mathbf{n}ds,$ (A.2)

where $T^{o}$ is the interior of the cell $T\in\mathcal{T}$. Summing over all
$T\in\mathcal{T}$, one will get (A.1). ∎

Based on the above lemma, we can prove the following result.

######  Theorem A.1.

Assume $\mathbf{K}_{m}$ is a symmetric positive-definite constant matrix and
$\frac{k_{b}}{a}$ is a positive constant scalar. The stiffness matrix in the
method (3.5) is symmetric positive-definite.

###### Proof.

Let $\phi_{i}\in\overline{V}_{h}$ the Lagrange basis and
$B_{i}\in\overline{\mathcal{B}}$ the box associated with the node $i$, for
$i=1,\ldots,N$, where $N=\\#\overline{\mathcal{B}}$. One can calculate that,
for the first term in (3.5),

$\begin{split}-\int_{\partial
B_{i}\setminus\Gamma}\mathbf{K}_{m}\nabla\phi_{j}\cdot\mathbf{n}ds=\sum_{T\in\mathcal{T}}\int_{T}\mathbf{K}_{m}\nabla\phi_{j}\cdot\nabla\phi_{i}dxdy=-\int_{\partial
B_{j}\setminus\Gamma}\mathbf{K}_{m}\nabla\phi_{i}\cdot\mathbf{n}ds.\end{split}$

As for the second term in (3.5), it is non-zero only if the nodes $i$ and $j$
are at the same vertex or on the same edge on $\Gamma$. We discuss three
possibilities ($i\neq j$):

  1. 1.

If the nodes $i$ and $j$ are on the same vertex, then

$-\int_{\partial
B_{i}\cap\Gamma}\frac{k_{b}}{a}(\phi_{j}^{+}-\phi_{j}^{-})ds=-\frac{3}{4}\int_{\partial
B_{i}\cap\partial B_{j}}\frac{k_{b}}{a}ds=-\int_{\partial
B_{j}\cap\Gamma}\frac{k_{b}}{a}(\phi_{i}^{+}-\phi_{i}^{-})ds.$

  2. 2.

If the nodes $i$ and $j$ are two vertices of an edge $e$ and on the same side
of $\Gamma$, then

$-\int_{\partial
B_{i}\cap\Gamma}\frac{k_{b}}{a}(\phi_{j}^{+}-\phi_{j}^{-})ds=\frac{1}{8}\int_{e}\frac{k_{b}}{a}ds=-\int_{\partial
B_{j}\cap\Gamma}\frac{k_{b}}{a}(\phi_{i}^{+}-\phi_{i}^{-})ds.$

  3. 3.

If the node $i$ and $j$ are two vertices of an edge $e$ and on the two sides
of $\Gamma$, then

$-\int_{\partial
B_{i}\cap\Gamma}\frac{k_{b}}{a}(\phi_{j}^{+}-\phi_{j}^{-})ds=-\frac{1}{8}\int_{e}\frac{k_{b}}{a}ds=-\int_{\partial
B_{j}\cap\Gamma}\frac{k_{b}}{a}(\phi_{i}^{+}-\phi_{i}^{-})ds.$

In all cases, we have

$-\int_{\partial
B_{i}\cap\Gamma}\frac{k_{b}}{a}(\phi_{j}^{+}-\phi_{j}^{-})ds=-\int_{\partial
B_{j}\cap\Gamma}\frac{k_{b}}{a}(\phi_{i}^{+}-\phi_{i}^{-})ds.$

Therefore, the stiffness matrix of the method (3.5) is symmetric.

Take an arbitrary grid function $\\{v_{i}\\}_{i=1}^{N}$, and denote
$v=\sum_{i=1}^{N}v_{i}\phi_{i}\in\overline{V}_{h}$. One can calculate that,
for the first term in (3.5),

$\begin{split}-\sum_{i,j=1}^{N}v_{i}v_{j}\int_{\partial
B_{i}\setminus\Gamma}\mathbf{K}_{m}\nabla\phi_{j}\cdot\mathbf{n}ds&=-\sum_{i=1}^{N}v_{i}\int_{\partial
B_{i}\setminus\Gamma}\mathbf{K}_{m}\nabla v\cdot\mathbf{n}ds\\\
&=\sum_{i=1}^{N}v_{i}\sum_{T\in\mathcal{T}}\int_{T}\mathbf{K}_{m}\nabla
v\cdot\nabla\phi_{i}dxdy\\\
&=\sum_{T\in\mathcal{T}}\int_{T}\mathbf{K}_{m}\nabla v\cdot\nabla vdxdy\geq
0,\end{split}$

where the equality holds only if $v$ is constant on every subdomain separated
by $\Gamma$ (if no area is enclosed by $\Gamma$, then $v$ must be a constant
on the whole domain.)

To facilitate the discussion on the second term in (3.5), we define
$\mathcal{E}^{\Gamma}=\\{e:e=\partial T\cap\Gamma,\,\forall
T\in\mathcal{T}_{h}\\}$ the collection of edges on $\Gamma$, and let
$[\\![v]\\!]_{a}=v^{+}_{a}-v^{-}_{a},[\\![v]\\!]_{b}=v_{b}^{+}-v_{b}^{-}$ the
jumps of $v$ at $a$ and $b$ on the two sides of $e\in\mathcal{E}^{\Gamma}$,
where $a$ and $b$ are the two vertices of $e$ and $\pm$ are the limits
specified for every $e$. Then, one can calculate that,

$\begin{split}&-\sum_{i,j=1}^{N}v_{i}v_{j}\int_{\partial
B_{i}\cap\Gamma}\frac{k_{b}}{a}(\phi_{j}^{+}-\phi_{j}^{-})ds=-\sum_{i=1}^{N}v_{i}\int_{\partial
B_{i}\cap\Gamma}\frac{k_{b}}{a}(v^{+}-v^{-})ds\\\
&=\frac{1}{2}\sum_{e\in\mathcal{E}^{\Gamma}}\int_{e}\frac{k_{b}}{a}\left([\\![v]\\!]_{a}(\frac{3}{4}[\\![v]\\!]_{a}+\frac{1}{4}[\\![v]\\!]_{b})+[\\![v]\\!]_{b}(\frac{1}{4}[\\![v]\\!]_{a}+\frac{3}{4}[\\![v]\\!]_{b})\right)ds\\\
&=\frac{1}{8}\sum_{e\in\mathcal{E}^{\Gamma}}\int_{e}\frac{k_{b}}{a}\left(2[\\![v]\\!]_{a}^{2}+2[\\![v]\\!]_{b}^{2}+([\\![v]\\!]_{a}+[\\![v]\\!]_{b})^{2}\right)ds\geq
0,\end{split}$

where the equality holds only if the jump of $v$ vanish on $\Gamma$, i.e., $v$
is continuous across $\Gamma$.

Summing up the two inequalities, one can conclude that

$-\sum_{i,j=1}^{N}v_{i}v_{j}\left(\int_{\partial
B_{i}\setminus\Gamma}\mathbf{K}_{m}\nabla\phi_{j}\cdot\mathbf{n}ds+\int_{\partial
B_{i}\cap\Gamma}\frac{k_{b}}{a}(\phi_{j}^{+}-\phi_{j}^{-})ds\right)\geq 0,$

where the equality holds only if $v$ is constant on the whole domain.
Therefore, the stiffness matrix of the method (3.5) is positive-definite, if
the Dirichlet condition on a non-zero measure boundary is specified. ∎

## References

  * [1] Aavatsmark, I., Barkve, T., Bøe, O. and Mannseth, T., 1998. Discretization on unstructured grids for inhomogeneous, anisotropic media. Part I: Derivation of the methods, SIAM Journal on Scientific Computing, 19(5), pp.1700-1716.
  * [2] Aavatsmark, I., Barkve, T., Bøe, O. and Mannseth, T., 1998. Discretization on unstructured grids for inhomogeneous, anisotropic media. Part II: Discussion and numerical results, SIAM Journal on Scientific Computing, 19(5), pp.1717-1736.
  * [3] A.S. Abd and A.S. Abushaikha, On the performance of the node control volume finite element method for modeling multi-phase fluid flow in heterogeneous porous media, Transport in Porous Media, 135, (2020), pp.409-429.
  * [4] R. Ahmed, M.G. Edwards, S. Lamine, B.A. Huisman and M. Pal, Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model, Journal of Computational Physics, 284 (2015), pp.462-489.
  * [5] P. Angot, F. Boyer and F. Hubert, Asymptotic and numerical modelling of flows in fractured porous media, ESAIM: Mathematical Modelling and Numerical Analysis, 43, (2009), pp.239-275.
  * [6] P.F. Antonietti, C. Facciola, A. Russo and M. Verani, Discontinuous Galerkin approximation of flows in fractured porous media on polytopic grids, SIAM Journal on Scientific Computing, 41 (2019), pp.A109-A138.
  * [7] R.G. Baca, R.C. Arnett and D.W. Langford, Modelling fluid flow in fractured‐porous rock masses by finite‐element techniques, International Journal for Numerical Methods in Fluids, 4 (1984), pp.337-348.
  * [8] R.E. Bank and D.J. Rose, Some error estimates for the box method, SIAM Journal on Numerical Analysis, 24, (1987), pp.777-787.
  * [9] G.I. Barenblatt, I.P. Zheltov and I.N. Kochina, Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata], Journal of Applied Mathematics and Mechanics, 24 (1960), pp.1286-1303.
  * [10] I. Berre et al., Verification benchmarks for single-phase flow in three-dimensional fractured porous media, Advances in Water Resources, 147 (2021), p.103759.
  * [11] I.I. Bogdanov, V.V. Mourzenko, J.F. Thovert and P.M. Adler, Effective permeability of fractured porous media in steady state flow, Water Resources Research, 39, (2003).
  * [12] W.M. Boon, J.M. Nordbotten and I. Yotov, Robust discretization of flow in fractured porous media, SIAM Journal on Numerical Analysis, 56 (2018), pp.2203-2233.
  * [13] K. Brenner, M. Groza, C. Guichard, G. Lebeau and R. Masson, Gradient discretization of hybrid dimensional Darcy flows in fractured porous media, Numerische Mathematik, 134 (2016), pp.569-609.
  * [14] K. Brenner, J. Hennicker, R. Masson, P. Samier, Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix–fracture interfaces, IMA Journal of Numerical Analysis, 37 (2017), pp.1551-1585.
  * [15] S. Chen, Discontinuous Galerkin method for hybrid-dimensional fracture models of two-phase flow, 488 (2023), Journal of Computational Physics, p.112244.
  * [16] L. Chen, A new class of high order finite volume methods for second order elliptic equations, SIAM journal on numerical analysis, 47 (2010), pp.4021-4043.
  * [17] T. Chen, C. Clauser, G. Marquart, K. Willbrand and D. Mottaghy, A new upscaling method for fractured porous media, Advances in Water Resources, 80 (2015), pp.60-68.
  * [18] C. D’Angelo and A. Scotti, A mixed finite element method for Darcy flow in fractured porous media with non-matching grids, ESAIM: Mathematical Modelling and Numerical Analysis, 46 (2012), pp.465-489.
  * [19] B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti, I. Stefansson and A. Tatomir, Benchmarks for single-phase flow in fractured porous media, Advances in Water Resources, 111 (2018), pp.239-258.
  * [20] B. Flemisch, A. Fumagalli and A. Scotti, A review of the XFEM-based approximation of flow in fractured porous media, Advances in Discretization Methods: Discontinuities, Virtual Elements, Fictitious Domain Methods, (2016), pp.47-76.
  * [21] N. Frih, V. Martin, J.E. Roberts and A. Saâda, Modeling fractures as interfaces with nonmatching grids, Computational Geosciences, 16 (2012), pp.1043-1060.
  * [22] A. Fumagalli, L. Pasquale, S. Zonca and S. Micheletti, An upscaling procedure for fractured reservoirs with embedded grids, Water Resources Research, 52 (2016), pp.6506-6525.
  * [23] D. Gläser, R. Helmig, B. Flemisch and H. Class, A discrete fracture model for two-phase flow in fractured porous media, Advances in Water Resources, 110, (2017), pp.335-348.
  * [24] D. Gläser, M. Schneider, B. Flemisch and R. Helmig, Comparison of cell-and vertex-centered finite-volume schemes for flow in fractured porous media, Journal of Computational Physics, 448 (2022), p.110715.
  * [25] H. Guo, W. Feng, Z. Xu and Y. Yang, Conservative numerical methods for the reinterpreted discrete fracture model on non-conforming meshes and their applications in contaminant transportation in fractured porous media, Advances in Water Resources, 153 (2021), p.103951.
  * [26] W. Hackbusch, On first and second order box schemes, Computing 41 (1989) 277–296.
  * [27] H. Hægland, A. Assteerawatt, H.K. Dahle, G.T. Eigestad and R. Helmig, Comparison of cell-and vertex-centered discretization methods for flow in a two-dimensional discrete-fracture–matrix system, Advances in Water resources, 32, (2009), pp.1740-1755.
  * [28] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydrosystems, Springer-Verlag, 1997.
  * [29] H. Hoteit and A. Firoozabadi, Multicomponent fluid flow by discontinuous Galerkin and mixed methods in unfractured and fractured media, Water Resources Research, 41 (2005).
  * [30] H. Hoteit and A. Firoozabadi, An efficient numerical model for incompressible two-phase flow in fractured media, Advances in Water Resources, 31 (2008), pp.891-905.
  * [31] T. Kadeethum, H.M. Nick, S. Lee and F. Ballarin, Flow in porous media with low dimensional fractures by employing enriched Galerkin method, Advances in Water Resources, 142 (2020), p.103620.
  * [32] M. Karimi-Fard, A. Firoozabadi, Numerical simulation of water injection in 2D fractured media using discrete-fracture model, In SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, September (2001), SPE-71615.
  * [33] M. Karimi-Fard, L.J. Durlofsky and K. Aziz, An efficient discrete-fracture model applicable for general-purpose reservoir simulators, SPE Journal, 9 (2004), pp.227-236.
  * [34] J.G. Kim and M.D. Deo, Finite element, discrete‐fracture model for multiphase flow in porous media, AIChE Journal, 46 (2000), pp.1120-1130.
  * [35] T. Koch, D. Gläser, K. Weishaupt, S. Ackermann, M. Beck, B. Becker, S. Burbulla, H. Class, E. Coltman, S. Emmert and T. Fetzer, DuMux 3–an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling. Computers & Mathematics with Applications, 81 (2021), pp.423-443.
  * [36] M. Köppel, V. Martin, J. Jaffré and J.E. Roberts, A Lagrange multiplier method for a discrete fracture model for flow in porous media, Computational Geosciences, 23 (2019), pp.239-253.
  * [37] P.P. Leme and A. Spagnuolo, A multiple-porosity model for a single-phase flow through naturally-fractured porous media, Computation and Applied Mathematics, 17 (1998), p.19.
  * [38] L. Li, and S.H. Lee, Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media, SPE Reservoir evaluation & engineering, 11 (2008), pp.750-758.
  * [39] V. Martin, J. Jaffré and J.E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media, SIAM Journal on Scientific Computing, 26 (2005), pp.1667-1691.
  * [40] A. Moinfar, A. Varavei, K. Sepehrnoori and R.T. Johns, Development of an efficient embedded discrete fracture model for 3D compositional reservoir simulation in fractured reservoirs, SPE Journal, 19 (2014), pp.289-303.
  * [41] J.E.P. Monteagudo and A. Firoozabadi, Control‐volume method for numerical simulation of two‐phase immiscible flow in two‐and three‐dimensional discrete‐fractured media, Water Resources Research, 40 (2004).
  * [42] J.E.P Monteagudo and A. Firoozabadi, Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects, SPE Journal, 12 (2007), pp.355-366.
  * [43] J.E.P. Monteagudo and A. Firoozabadi, Comparison of fully implicit and IMPES formulations for simulation of water injection in fractured and unfractured media, International Journal for Numerical Methods in Engineering, 69 (2007), pp.698-728.
  * [44] V.V. Mourzenko, I.I. Bogdanov, J.F. Thovert and P.M. Adler, Three-dimensional numerical simulation of single-phase transient compressible flows and well-tests in fractured formations, Mathematics and Computers in Simulation, 81, (2011), pp.2270-2281.
  * [45] I. Mozolevski, M.A. Murad and L.A. Schuh, High order discontinuous Galerkin method for reduced flow models in fractured porous media, Mathematics and Computers in Simulation, 190 (2021), pp.1317-1341.
  * [46] J. Noorishad and M. Mehran, An upstream finite element method for solution of transient transport equation in fractured porous media, Water Resources Research, 18 (1982), pp.588-596.
  * [47] M. Oda, Permeability tensor for discontinuous rock masses, Geotechnique, 35, (1985), pp.483-495.
  * [48] V. Reichenberger, H. Jakobs, P. Bastian and R. Helmig, A mixed-dimensional finite volume method for two-phase flow in fractured porous media, Advances in Water Resources, 29 (2006), pp.1020-1036.
  * [49] T.H. Sandve, I. Berre and J.M. Nordbotten, An efficient multi-point flux approximation method for discrete fracture–matrix simulations, Journal of Computational Physics, 231 (2012), pp.3784-3800.
  * [50] P. Schädle, P. Zulian, D. Vogler, S.R. Bhopalam, M.G. Nestola, A. Ebigbo, R. Krause and M.O. Saar, 3D non-conforming mesh model for flow in fractured porous media using Lagrange multipliers, Computers & Geosciences, 132 (2019), pp.42-55.
  * [51] N. Schwenck, B. Flemisch, R. Helmig and B.I. Wohlmuth, Dimensionally reduced flow models in fractured porous media: crossings and boundaries, Computational Geosciences, 19 (2015), pp.1219-1230.
  * [52] M. Ţene, S.B. Bosma, M.S. Al Kobaisi and H. Hajibeygi, Projection-based embedded discrete fracture model (pEDFM), Advances in Water Resources, 105 (2017), pp.205-216.
  * [53] A. Vogel, J. Xu and G. Wittum, A generalization of the vertex-centered finite volume scheme to arbitrary high order, Computing and visualization in science, 13 (2010), pp.221-228.
  * [54] J.E. Warren and P.J. Root, The behavior of naturally fractured reservoirs, Society of Petroleum Engineers Journal, 3 (1963), pp.245-255.
  * [55] X. Wu, H. Guo, Z. Xu and Y. Yang, A reinterpreted discrete fracture model for Darcy–Forchheimer flow in fractured porous media, Advances in Water Resources, 179 (2023), p.104504.
  * [56] Y.S. Wu and K. Pruess, A multiple-porosity method for simulation of naturally fractured petroleum reservoirs, SPE Reservoir Engineering, 3 (1988), pp.327-336.
  * [57] Z. Xu and Y. Yang, The hybrid dimensional representation of permeability tensor: A reinterpretation of the discrete fracture model and its extension on nonconforming meshes, Journal of Computational Physics, 415 (2020), p.109523.
  * [58] Z. Xu, Z. Huang and Y. Yang, The hybrid-dimensional Darcy’s law: A non-conforming reinterpreted discrete fracture model (RDFM) for single-phase flow in fractured media, Journal of Computational Physics, 473 (2023), p.111749.
  * [59] J. Zhao and H. Rui, A discrete fracture-matrix approach based on Petrov-Galerkin immersed finite element for fractured porous media flow on nonconforming mesh, Journal of Computational Physics, (2023), p.112718.
  * [60] A. Zidane and A. Firoozabadi, An efficient numerical model for multicomponent compressible flow in fractured porous media, Advances in Water Resources, 74 (2014), pp.127-147.


edge corresponding to $l\in[1,a_{ij}]$. The factors $\prod_{i}m_{i}!$,
$\prod_{ij}a_{ij}!$ and $\prod_{i}2^{a_{ii}}$ come from symmetry, the rest is
the usual Wick formula.

The total number of edges $e$ satisfies constraints:

$e\geq{3\over 2}v,\quad k=e-v\implies e\in\\{k,\dots,3k\\}\,,$

hence in the expression $A_{k}$ the propagator $(p^{ij})_{1\leq i,j\leq n}$
appears at most $3k$ times.

One-loop exactness is equivalent to an infinite sequence of differential
equations:

$A_{1}=0,~{}A_{2}=0,\dots$ (9.47)

Up to symmetry, the number of distinct graphs for $A_{1},A_{2},A_{3}$ is
$5,41,378$ respectively.

For example, 5 graphs appearing in $A_{1}$ are the following:

$\Gamma_{1}$$\Gamma_{2}$$\Gamma_{3}$$\Gamma_{4}$$\Gamma_{5}$$g$$g$$g$$g$$g$$f$$f$$f$$f$$f$$f$

and the expression $A_{1}$ is

${1\over 2}\sum_{i,j}p^{ij}\,\partial_{ij}g-{1\over
2}\sum_{i_{1}j_{1}i_{2}j_{2}}p^{i_{1}j_{1}}p^{i_{2}j_{2}}\,\partial_{i_{1}}g\,\partial_{j_{1}i_{2}j_{2}}f-\\\
-{1\over
8}\,g\sum_{i_{1}j_{2}i_{2}j_{2}}p^{i_{1}j_{1}}p^{i_{2}j_{2}}\,\partial_{i_{1}j_{1}i_{2}j_{2}}f+{1\over
8}\,g\sum_{i_{1}j_{2}i_{2}j_{2}i_{3}j_{3}}p^{i_{1}j_{1}}p^{i_{2}j_{2}}p^{i_{3}j_{3}}\,\partial_{i_{1}j_{1}i_{2}}f\,\partial_{j_{2}i_{3}j_{3}}f+\\\
+{1\over
12}\,g\sum_{i_{1}j_{1}i_{2}j_{2}i_{3}j_{3}}p^{i_{1}j_{1}}p^{i_{2}j_{2}}p^{i_{3}j_{3}}\,\partial_{i_{1}i_{2}i_{3}}f\,\partial_{j_{1}j_{2}j_{3}}f\,.$

## Acknowledgements

We are grateful to Robert Bryant and Joseph M. Landsberg for useful
discussions. We are grateful to Nikolai Perkhunkov for useful advises and help
with managing huge Maple computations. A.O. is grateful to IHES for
invitations and an excellent working atmosphere.

## References

  * [1] P. Etingof, D. Kazhdan, A. Polishchuk, When is the Fourier transform of an elementary function elementary?, Sel. math., New ser. 8 (2002) 27-66.
  * [2] J. L. Brylinski, Transformations canoniques, Dualité projective, Théorie de Lefschetz, Transformations de Fourier et sommes trigonométriques, Asterisque 140-141 (1986), 3-134.
  * [3] Encyclopedia of Special Functions: The Askey-Bateman Project, Tables of integral transforms, Volume II, McGraw-Hill Book Company, Inc. 1954, Section 8.6, formula 23.
  * [4] A. Coffman, A. Schwartz, and C. Stanton, The Algebra and Geometry of Steiner and other Quadratically Parametrizable Surfaces, Computer Aided Geometric Design (3) 13 (April 1996), p. 257-286.
  * [5] I. Dolgachev, Kummer Surfaces: 200 Years of Study, Notices of the American Mathematical Society 67(10):1.
  * [6] E. I. Wilczynski, Projective-differential geometry of curved surfaces, Trans. AMS, 8 (1907) 233-260; 9 (1908) 79-120, 293-315; 10 (1909) 176-200, 279-296.