Modelling and discretization of flow in porous media with thin, full-tensor permeability inclusions

When modelling fluid flow in fractured reservoirs, it is common to represent the fracturesas lower-dimensional inclusions embedded in the host medium. Existing discretizationsof flow in porous media with thin inclusions assume that the principal directions of theinclusion permeability tensor are aligned with the inclusion orientation. While this mod-elling assumption works well with tensile fractures, it may fail in the context of faults,where the damage zone surrounding the main slip surface may introduce anisotropy thatis not aligned with the main fault orientation. In this paper, we introduce a generalizeddimensional reduced model which preserves full-tensor permeability effects also in theout-of-plane direction of the inclusion. The governing equations of flow for the lower-dimensional objects are obtained through vertical averaging. We present a framework fordiscretization of the resulting mixed-dimensional problem, aimed at easy adaptation ofexisting simulation tools. We give numerical examples that show the failure of existingformulations when applied to anisotropic faulted porous media, and go on to show theconvergence of our method in both 2D and 3D


Introduction
Modeling and simulation of flow in porous media with faults, fractures, and other thin inclusions representing discontinuities is central to a wide range of subsurface engineering applications, including geothermal energy exploitation (Bödvarsson and Tsang, 1982), shale gas extraction (Cao et al., 2016), carbon sequestration (Johnson et al., 2009), and energy storage (Nagelhout and Roest, 1997).
The inclusions are characterized by a high aspect ratio, and permeability significantly different from that of the host medium; hence, they severely affect flow patterns. This poses a challenge for traditional simulation models, which are based on upscaling of finescale details into an equivalent permeability (Oda, 1985;Farmer, 2002;Liu et al., 2016;Saevik et al., 2013). We instead focus on an alternative approach, which explicitly represents the inclusions in the mathematical and simulation models and thereby to a large degree avoids challenges relating to parameter upscaling. To avoid elongated cells at the inclusion in the computational grid, it is common to represent the inclusions as codimension one objects embedded in the host medium Nordbotten et al., 2019). The intersection of inclusions further gives rise to line and point intersections of co-dimension two and three. Each of these objects (matrix, inclusions, and intersection points and lines) are represented as independent subdomains separated by interfaces. We refer to this representation of the geometry as mixed-dimensional.
Governing equations for fluid flow in lower-dimensional representation of the inclusion can be derived by integration in the direction orthogonal to the inclusion. This leads to a decomposition of the governing equations into an in-plane component that represents flow within the inclusion, and an out-of-plane component that couples flow between the inclusion and the host medium. While the in-plane flow has been modeled with both linear and non-linear, as well as both isotropic and non-isotropic flow models (Martin et al., 2005;Reichenberger et al., 2006;Brenner et al., 2017Brenner et al., , 2018, existing models for the coupling term are limited by an assumption on orthogonal flow between inclusion and host. Reduced order models for flow were also developed for aquifers, leading to the same set of equations, see for instance Bear (1979), Yortsos (1995), and Nordbotten and Celia (2011). These existing models will be denoted as "local" in the following, meaning that each partial differential equation (PDE) contains only quantities associated with the subdomain where the PDE is defined.
Local models generally work well when the inclusion is a joint (tensile fracture). However, inclusions with a more complex geological history may have significantly more complex flow properties in the out-of-plane direction. For instance, the damage zone in the vicinity of faults may exhibit shear fractures, slip surfaces, and/or deformation bands, as summarized in Fossen et al. (2007). These features introduce secondary permeability anisotropy in the damage zone as they tend to have preferred orientations, as shown by both field studies (Fossen et al., 2005;Johansen and Fossen, 2008) and core analysis (Hesthammer et al., 2000). This leads to preferential flow directions that are neither parallel nor orthogonal to the main plane. This type of flow cannot be represented by existing models that employ dimension reduction. To the Authors' best knowledge, the only attempt to modeling faults and their surrounding damage zones in a mixed-dimensional framework can be found in Fumagalli and Scotti (2019). However, they still apply local formulations to model the damage zones as lower-dimensional objects which are connected on one side to the fault and on the other side to the rock matrix, hence conceptually seeing the whole fault zone as a multilayer object. An alternative approach would be to implement the fault core as a transmissibility multiplier and the damage zone by modifying the grid permeability in the cells adjacent to the model faults, as illustrated in Wilson et al. (2020). In the following, we will consistently refer to the thin inclusions as faults, notwithstanding that all methods presented herein can be applied to models of fractures and other thin inclusions, however, we expect that the methods proposed are of more importance for faults.
The contribution of this paper is two-fold: First, we present a generalized dimensional reduced model that can preserve full-tensor permeability effects also in the out-of-plane direction of the fault. The resulting reduced equations have a form similar to that of traditional models, however the more general coupling structure leads to additional terms both in the in-plane and out-of-plane equations. These terms, as well as our whole novel formulation, will be denoted as "semi-local" in the following, emphasizing the fact that the new PDEs will contain quantities that, while physically in the same location, from a modeling perspective reside outside the subdomain where the PDE is defined, specifically the internal boundary between the subdomain and its higher dimensional neighbor.
Multiple discretization schemes have been proposed for the local dimensionally-reduced models, including methods based on finite volumes (Helmig et al., 1997;Karimi-Fard et al., 2003;Sandve et al., 2012), mixed finite elements (Martin et al., 2005;Boon et al., 2018), virtual elements  and mimetic methods (Formaggia et al., 2018). A comparison study for all these discretizations of flow in fractured media can be found in Flemisch et al. (2018) and Berre et al. (2020) for 2D and 3D flow, respectively. However, the additional terms arising in our formulation bring the semi-local model outside the scope of previously proposed discretization methods. The second contribution of the paper is therefore the derivation of discretization schemes for semi-local models. We achieve this in two stages: First, based on the unified framework for discretization of mixed-dimensional problems with local interface laws presented in Nordbotten et al. (2019), we present conditions under which any standard discretization scheme for fixed-dimensional problems can be extended to mixed-dimensional problems with semi-local interface laws. Second, we present a concrete discretization approach based on finite volume methods.
The paper is organized as follows. In Sec. 2, the mathematical model is presented, first for a domain with a single fault, and then for a general faults configuration. Thereafter, in Sec. 3, the unified discretization is formulated. After presenting simulation results in Sec. 4, concluding remarks are given in Sec. 5.

Flow modelling in faulted porous media
In this section, the mathematical model for flow in faulted porous media is presented, first for a porous domain containing a single fault (Sections 2.1 and 2.2), and then for a general network of faults (Section 2.3). For the general case, we also provide the weak formulation of the interface problem (Sections 2.4-2.5), which will be useful from the perspective of implementation.

Domain with a single fault
We start by considering two three-dimensional porous media Ψ 1 and Ψ 2 , each of them with its Neumann and Dirichlet boundaries and , respectively. The two threedimensional domains are separated by a fault Ψ 3 , which is a thin, almost two-dimensional object of thickness (in the following will be denoted as the aperture), but which is Figure 1. Representation of the fault as a thin three-dimensional domain Ψ 3 (left) and as a two-dimensionl manifold Ω 3 (right). The boundary of Ψ adjacent to Ψ 3 is denoted by Ψ 3 Ψ , for = 1, 2, while is the normal vector which is always pointing outwards from Ψ , for = 1, 2, 3.
currently represented as three-dimensional. We note that Ψ 3 need not be planar, i.e.
need not be constant. We denote by Ψ 3 Ψ , for = 1, 2, the boundary of Ψ adjacent to Ψ 3 . Furthermore, let be the normal vector which is always pointing outwards from Ψ . It thus follows that 3 = − on Ψ 3 Ψ . A representation of the fault as a thin three-dimensional domain Ψ 3 is illustrated in the left of Fig. 1. Darcy flow in the three-dimensional porous medium is then governed by the following set of equations ( = 1, 2, 3): Here, is pressure, q is the Darcy flux, is a source, and is a second-order tensor representing the absolute permeability divided by fluid viscosity. Equation (1) represents mass conservation, while equation (2) is Darcy's law. Equation (3) represents flux continuity conditions on Ψ 3 Ψ , where 3, represents flow from Ψ 3 to Ψ , thus by flux continuity it follows that 3, = − ,3 . Finally, equations (4)-(5) are boundary conditions on Ψ and Ψ , repectively.
Before deriving the governing equations for the lower-dimensional manifold, we discuss the decomposition of the permeability tensor within the fault. Existing local laws for faults as embedded thin inclusions assume that the principal directions of the local permeability tensor are aligned with the fault orientation, as illustrated in Fig. 2.a. Hence, more general orientations of the principal permeability directions, shown in Fig. 2.b-2.c, cannot be represented by existing models. To be concrete, we let the permeability on Ψ 3 have the following decomposition in terms of a coordinate system aligned with the fault orientation: Here, 3, is a 2 × 2 second-order tensor representing the within-fault permeability and 3,⊥ is a scalar representing the normal permeability. The off-diagonal term 3, is a twovector representing the symmetric off-diagonal components of 3 ; for local interface laws, these off-diagonal components are assumed to be negligible, i.e. 3, = 0 (Nordbotten and Celia, 2011;Berre et al., 2020). The inclusion of this anisotropic term leads to significant complications in the modeling and discretization, and is the main topic of this work. With this structure of the fault permeability, the Darcy flux for the fault can be decomposed , where the 2-vector tangential component 3, and the scalar normal component 3,⊥ have the following form: Here, ∇ and ∇ ⊥ = represent the in-plane and out-of-plane components of the gradient for the fault, respectively. Figure 3. Illustration of the mixed-dimensional geometry. Ω 3 is connected to the higher dimensional neighbors Ω through the interfaces Γ ,3 , for = 1, 2. Note that Ω 3 , Γ ,3 and Ω 3 Ω are all coinciding in physical space.

Model reduction
To proceed, we apply integration over the perpendicular direction to achieve a dimension reduction of the fault. This replaces Ψ 3 with a lower-dimensional domain Ω 3 (see right of Fig. 1). Note that we use Ψ to represent the equi-dimensional geometry, that is all Ψ are 3D, and Ω to denote the mixed-dimensional geometry. We also introduce two interfaces Γ ,3 on each side = 1, 2 of Ω 3 , as illustrated in Fig. 3. The interfaces physically represent the half zone comprised between the fault and either side of the surrounding matrix. In a mixed-dimensional setting, they have no perpendicular extent, and serve as connectors between two objects of different dimensions. Note that, due to the dimension reduction of the model, Ω 3 , Γ 1,3 , Γ 2,3 , Ω 3 Ω 1 and Ω 3 Ω 2 are all coinciding in physical space. Furthermore, we define the integrated Darcy flux (2) 3 and the average pressure Here, we use subscripts to index the domains, and superscripts (when necessary for clarity) to emphasize the effective topological dimension of the domain, e.g. (1) and (7) along the perpendicular direction, the out-of-plane component of the gradient is converted into a jump operator as follows: The governing equations for the fault are then obtained from equations (1), (7), (4) and (5) by integrating in the perpendicular direction. Moreover, since the fault is assumed to be thin, we assume that the permeability is constant across the perpendicular direction.
Together with the definitions above, this results in where we have also introduced the integrated source term and boundary flux We emphasize that the differential operator ∇ 3 in eqs. (11)-(12) operates on the manifold Ω 3 . Compared to traditional upscaled models, see for instance Nordbotten et al. (2019), additional terms ,3 appear in equation (12), analogous to the flux terms ,3 in equation (11). This two-vector term, which is not present in previous work, represents the withinfault flux induced by pressure differences between the fault and the surrounding matrix, and is defined for either side of the fault as where the permutation variable ,3 is positive if the coordinate systems of Ω 3 and Ω 3 Ω coincide, and negative otherwise.
To complete the model, we derive a constitutive law for ,3 . This is obtained by integrating equation (8) in the perpendicular direction, that is The left hand side of equation (17) is approximated using the trapeizodal rule, that is where continuity of the flux across the boundary between the fault and the surrounding matrix is applied. The first term at the right hand side of equation (17) is approximated Finally, the second term at the right hand side of (17) is resolved using the jump operator defined in equation (10) as follows: By incorporating eqs. (18), (19) and (20) into equation (17), we identify the flux ,3 having the following form: Here, the first term on the right-hand side represents the local component of the constitutive law, while the second part is the semi-local contribution that induces a flux across Γ ,3 due to the pressure gradient within the lower-dimensional manifold Ω 3 .
Inspecting equations (16) and (21), we see that both the normal permeability 3,⊥ and the off-diagonal permeability 3, are in the reduced model naturally interpreted as properties of the interface Γ ,3 . In the continuation, we will thus generalize the model as derived above, and index these quantities with the interface, i.e. 3, → 3, , and 3,⊥ → 3, ,⊥ are assigned independently to either side of the fault.
We remark that due to the model reduction, the within-fault permeability 3, and the normal permeability 3, ,⊥ scale with the aperture and its inverse, respectively, while the off-diagonal permeability 3, , remains as in the equi-dimensional model. In order to present equations (22)-(27) without reference to this small parameter, these scalings have been incorporated directly into the material constants. Thus, the mixed-dimensional permeability 3 is related to the equi-dimensional 3 as follows We point out that, when one reduces multiple dimensions at once, these scalings get exponents corresponding to the number of dimensions below the ambient dimension. We also emphasize that the normal and off-diagonal permeabilities are in principle not a property of the fault itself, but instead a property which belongs to the internal interface Γ ,3 between the fault and either side of the higher-dimensional neighbors. This represents an important extension of the existing local laws for fractured porous media, making the model also applicable to faulted porous media, since it allows for capturing the anisotropic character of the fault damage zone. Moreover, since different values of 3, , and 3, ,⊥  = . Hence, we can write for all

Domain with a general network of faults
whereˇis the set of neighbors of Ω of lower dimension, e.g.ˇ1 = {Ω 2 , Ω 3 , Ω 4 }. It is easy to show that as long as the mixed-dimensional permeabilities are diagonally dominant in the sense of , ,⊥ det , > , , · , , ,

Mixed-dimensional formulation of the fault-matrix flows
While equations (29)-(34) constitute a full semi-local model, they are stated in a form which is not immediately amenable for discretization. This and the following subsection explore the model in more detail, with a goal of rewriting the equations in a form that can be handled by standard discretization schemes with only minimal adaptations. A discretization approach based on this reformulation is then given in Section 3.
In order to simplify the exposition, we will introduce a mixed-dimensional notation following Nordbotten et al. (2019). In particular, we will denote the collection of pressure functions as = 1 , ..., | | , and similarly the collection of all fluxes (both in domains and across boundaries) as = 1 , ..., | | , 1,1 , ..., | |,| | . It is sometimes convenient to refer explicitly to only the domain or boundary fluxes, and we will therefore sometimes abuse notation and simply write = ( , ). We refer to these as mixed-dimensional functions, and consistently denote them with calligraphic font. We adopt the natural convention that when evaluating a mixed-dimensional function at a point, say ∈ Ω , then we simply evaluate the function on that domain, so that ( ) = ( ). In a similar sense, we denote the disjoint union of domains as = ( Ω ) , Γ , .
With this notion of mixed-dimensional functions, the extension of the divergence and gradient operators to the mixed-dimensional setting is natural. First, we extend the concept of continuous functions by requiring that for to be continuous, then it must hold that, for all Γ , , · = , . Then, for any point ∈ Ω we define while for any point on an interface ∈ Γ , we define Now we can write equations (29) -(34) simply as: where we have also introduced the collection of sources = 1 , ..., | | and the collection of boundary fluxes = 1 , ..., | | . Here, the material coefficients are now all part of the mixed-dimensional permeability , which is defined such as that for any mixeddimensional gradient = = ( , ), it holds that for any point ∈ Ω : while for any point on an interface ∈ Γ , , it holds that ( ) ( ) = , ,⊥ , + , , , · .
It is then also sometimes convenient to write equation (39) in matrix form, that is for = ( , ) and = = ( , ), one has: Equation (44) highlights the contribution from the semi-local terms in the mixed-dimensional version of Darcy's law.

Weak formulation as an interface system
The semi-local terms in equations (29)-(34) lead to coupling terms between domains that are local in physical space, but non-local in the mixed-dimensional representation of the geometry. A critical example are the fault and its sides, which, from the perspective of implementation, we would prefer to only interact via the interfaces Γ , , and not directly, as is the case for the last term in equation (30).
Thus we are motivated to consider a reformulation of the governing equations before considering numerical discretizations. We proceed by first performing an LU decomposition of equation (44) as follows: where and are defined, respectively, as: and Ω is the Schur-complement defined as Note that, since ΓΓ consists only of scalar values ( , ,⊥ ), this reformulation only depends on the trivial inversion of scalars.
In the following it will be helpful to discuss the components of the mixed-dimensional gradient and divergence, and we therefore additionally define the "full jump" such that for any point ∈ Ω it holds that while the "half jump" ★ is simply the restriction of to Γ , . We then write (with the natural extension of ∇ and ∇·): We now proceed by (formally) eliminating internal domain variables, in order to obtain a problem only posed on interfaces. We note that equations (38) and (39) can now be written as the first order system: where use of equation (45) has been made. By writing out equation (50) in local notation for each Ω and by stating equation (51) explicitly as two equations, we obtain the following set of equations: This reveals that equations (52) and (53) form a locally well-posed system (of standard Darcy type) on each Ω , and we can therefore consider = ( ) for any given .
We formalize this concept by introducing the (continuous) solution operators for the standard elliptic value problem on Ω , S Ω , defined as: where is the solution to where Ω is the global boundary and = 1 |Ω | if Ω ∩ Ω ≠ , and zero otherwise. Using this solution operator, we see that the solution to equations (52) and (53) can be stated as functions of (and a set of number of numbers 0 corresponding to the domains where Ω ∩ Ω ≠ ) as: Inserting = ( , 0 ) etc. into equation (54), we have now reformulated the faultmatrix problem into a pure interface problem. From the perspective of implementation, we desire to consider the interface problem in the weak sense, and we therefore multiply by test functions and integrate to obtain the problem: Find ∈ 2 (Γ) such that, for and ( , 0 ) = 0 if Ω ∩ Ω ≠ . We point out that the inner products in equation (62) are bounded from a formal perspective, since for ∈ 2 (Γ), then ∈ 1 (Ω ), and both −1 ΓΓ ΓΩ ∇ and tr will lie in (at least) 2 (Γ , ).
Finally, we emphasize that equations (61)-(62) are attractive from the perspective of implementation, since the inner products appearing are easy to evaluate, and the solution operators S Ω can be approximated by any standard method, as we will detail in the next section.

Discretizations of flow for faulted porous media
The equations derived in Section 2.5, and in particular the interface problem of equation (61), form the starting point for the discretization approach laid out in this section. We present the general discretization framework in Section 3.1, and discuss implementational aspects in Section 3.2.

Unified discretization
Equation (61) provides a solution operator for the arbitrary standard method used to solve the elliptic boundary value problem (52)-(53) on Ω . To be concrete, we consider each domain Ω and its Neumann boundary Ω = Ω ∪ ∈ˇΩ Ω as endowed with a numerical discretization. Then, the solution operator S can be stated as S : where (Ω ), (Ω ), and ( Ω ) are the discrete representations of 2 (Ω ), 2 (Ω ) , and 2 ( Ω ), respectively, and is the topological dimension of Ω . In particular, S takes as input sinks, vector sources, and Neumann data and returns as output pressures, pressure gradients, and pressure traces. Most discretization schemes for elliptic equations can provide such a solution operator; we discuss the concrete implementation in the next subsection.
To discretize the flux coupling term , , we introduce a mortar-like grid T , on the interface Γ , on which the boundary flux , will be defined. The flux variables are represented as piecewise constant on the mortar grid T , , thus , ∈ 0 (T , ) ⊂ L 2 (Ω ).
In order to allow communications between subdomains, and thus explicitly relate the degrees of freedom of the numerical methods S and the mortar grids T , , we introduce projection operators, namely Π (Ω ) and Π 2 (Ω ) . The former is the compound operator projecting from the coupling variables on the mortar grids to the subdomain degrees of freedom, that is while the latter conversely moves from the numerical variables to the coupling variables, that is Now, following the variational formulation derived in Sec. 2.5, we exploit equation (62) in order to provide discretization-independent framework for faulted porous media. This takes the form: for given numerical discretizations S , find , ∈ 0 (T , ), for all ∈ and ∈ˆsuch that ★ , subject to discrete constraints (for all ∈ ): where the dummy variables , and have the interpretations of sources, forces, and fluxes due to interactions with other domains, respectively. In contrast, the variables , , and are the pressures, pressure gradients, and pressure traces after projection onto the grids T , .
The interpretation of this scheme is as follows. Eq. (67) resolves the internal differential equations in each subdomain, eq.(68) is the projection of variables from the flux grids to the numerical boundary (and source) data, while equation (66) simply states that the flux , between the fault and its surrodundings should satisfy Darcy's law. In the following section, we present the strategy for implementation of this approach and give details for a specific numerical scheme.

MPFA discretization
It is of interest to consider the requirements put on the subdomain solution operators S in some more detail. From the variational formulations stated above, we see that for a discretization on a generic subdomain Ω to interact with the interface Γ , we need to provide operators which: 1. Handle Neumann boundary data on the form Π (Ω ) , for all interfaces Γ where Ω is the higher-dimensional neighbor.
3. Provide a discrete operator tr so that Π 2 (Ω ) can project the pressure trace from Ω to Γ where Ω is the higher-dimensional neighbor.
4. Provide a pressure so that Π 2 (Ω ) can project the pressure to all Γ where Ω is the lower-dimensional neighbor.
6. Provide a pressure gradient so that Π 2 (Ω ) can project the pressure gradient to all Γ where Ω is the lower-dimensional neighbor.
The four first requirements are readily available for any discretization scheme for elliptic equations. Specifically, we have based our solution operators on a cell-centered finite volume method termed the multi-point flux approximation (MPFA) (Aavatsmark, 2002;Nordbotten and Keilegavlen, 2020). Treatment of vector source terms (item 5) is not as natural in primal discretization schemes such as finite elements, but is easy to include in most flux-based discretization methods such as e.g. mixed finite elements. We have employed the approach introduced in Starnoni et al. (2019), which treats the vector source term as part of the discrete divergence operator, and thereby provides an expression of the fluxes in terms of jumps in cell-centers vector sources. Finally, the pressure gradients are discretized as piece wise constant on each cell from an interpolation of the face cells fluxes (item 6). We implemented our model in PorePy, an open-source software for simulation of multiphysics processes in fractured porous media (Keilegavlen et al., 2021).
To better understand the structure of the discrete coupling, it is instructive to write out the coupled system for two subdomains Ω ℎ and Ω separated by an interface Γ (see Fig. 6). Let ℎ and , be the vectors of cell-center pressures in Ω ℎ and Ω respectively, and let be the vector of discrete mortar fluxes in Γ . The discrete coupled system in absence of external sources can then be represented on the generic form The first two rows of the system (69) represent the discretised differential equations in Figure 6. Illustration of a coupling between subdomains. Ω ℎ and Ω are the higher and lower subdomains respectively, Γ is the interface between the two subdomains, Ω ℎ is the portion of the boundary of Ω ℎ as seen from Γ , Π (Ω ) is the projection operator from coupling variables on the mortar grid to each of the subdomains degrees of freedom ( = ℎ, ), and Π 2 (Ω ) is the projection operator from numerical variables to coupling variables.
each subdomain, while the third row is the discretized Darcy's law in the direction perpendicular to the interface. Here, ℎ and are the fixed-dimensional discretizations on the subdomains, ℎ is the discretization of Neumann boundary conditions on Ω ℎ , is the discretization of source terms in Ω , is the discretization of the vector source term on Ω , is the discretized ΩΓ −1 ΓΓ product on Γ , and Π (Ω ℎ ) and Π (Ω ) are the projection operators from coupling variables on the mortar grid to each of the subdomains degrees of freedom. Furthermore, ℎ provides a discrete representation of the pressure trace operator on Ω ℎ , gives the pressure unknowns on Ω , gives the reconstruction of the pressure gradient on Ω , and Π 2 (Ω ) is the projection operator from numerical variables to coupling variables. Finally, is the discretized inverse normal permeability on Γ .
We conclude by making two remarks: firstly, there is no direct coupling between Ω ℎ and Ω and secondly, global boundary conditions are left out of the system.

Numerical examples
We validate the semi-local model and our implementation by a suite of numerical examples. First, we consider a case with a single fault, and show how the semi-local model can capture the effects of anisotropic off-diangonal permeabilities, while the local model fails to do so. Second, we probe the robustness of our discretization on more complex

Comparison to the equi-dimensional model
In this first example, we compare our reduced model to an equi-dimensional model. The aim is to highlight the enhanced modelling capabilities of our formulation with respect to the standard local formulation. With reference to this latter point, we present results of two test cases: the first one where the fault has the same off-diagonal permeability on both sides (see Fig. 2.b), and a second one where different permeability structures are assigned to each side of the fault (see Fig. 2.c).

Case 1: homogeneous permeability
We consider a 2D square domain of side = 1 cut by a horizontal fault of aperture = 1 located in the middle of the domain. In the mixed-dimensional setting we therefore have two 2D domains Ω 1 and Ω 2 and one 1D fault Ω 3 , as illustrated in Fig. 7.
The hydraulic conductivity is isotropic homogeneous for the 2D matrix, that is = I, with = 1, 2, while for the fault we consider the following equi-dimensional full tensor: For simplicity, we take 1 = 2 = . Boundary conditions consist of an applied difference in hydraulic head along the vertical direction and no-flow conditions elsewhere.
In particular, the inlet pressure ℎ is specified on the portion of the bottom boundary where 0.25 < < 0.75 , while the outlet pressure ℎ is specified on the portion of the where Δ is the size of the fault element in the reduced model, and , is calculated from the equi-dimensional model as the mean value of the two fault cells at each location : Convergence results are shown in Fig. 8a. As Fig. 8a clearly shows, our formulation presents about first-order convergence rate, while the local formulation does not converge. This is due to the strong anisotropy of the fault, which is not captured by the standard local formulation. As a result of the anisotropy of the fault, the flow will take a preferential direction towards one of the two inlets, therefore breakig the symmetry of the local formulation. This is better observed in Fig. 8b showing the pressure distribution along the fault for the three models. As Fig. 8b clearly shows, the semi-local and the equidimensional models coincide, while the local formulation exhibits an erroneous symmetric profile.

Case 2: dual permeability
As a further illustration of the enhanced modeling capabilities of the semi-local model, we modify the setup used in the previous section to have different permeability structures on the two sides of the fault. This is relevant for modeling of geological faults, where the two sides of the fault may undergo different damage processes. To that end, we divide the fault into an upper and lower part (see Fig. 9) and assign different permeability structures to the two sides, that is for = 1, 2: In particular, values of , , and ,⊥ are the same as those given in Table 1, while ,1, and ,2, take values of 50 and 80 / , respectively. The aperture of the fault is set to = 2 and we use the same boundary conditions as in Case 1.
Convergence results for the local and semi-local models are shown in Figs. 10a-10b, with the reference solution again computed from an equi-dimensional model with a grid with 40k cells. As in the previous case, the local model fails to converge, while the semilocal model exhibits first order convergence up to the last refinement step. Here, the mesh size is of the same order of the fault aperture, thus further error reduction cannot be expected due to the modeling error in the dimension reduction.

Self-convergence
In this section, we test the robustness of the method on more challenging fault configurations in 2D and 3D.

2D case
We consider the same test case as Case 1 in Boon et al. (2018). The domain is a unit square including a network of five faults (Fig. 11a). Of these five faults, one cuts the square domain into two 2D subdomains, denoted as Ω 1 and Ω 2 , respectively. The faults are numbered for = 3, .., 7 and are of two kinds: Ω 3 and Ω 4 are conductive, that is 3 = 4 = ,1 , while the other three are blocking, that is 5 = 6 = 7 = ,2 . The hydraulic conductivity is isotropic homogeneous for the 2D matrix, with 1 = 2 = , while for the faults we consider an equi-dimensional full tensor with , = 0.1 , , for = 3, .., 7. Boundary conditions consist of an applied difference in hydraulic head along the vertical direction and no-flow conditions elsewhere. Data for the simulations are reported in Table 2 The convergence results, shown in Fig. 11b, indicate a rate of at least first order.
The test thus confirms the performance of our method also in cases that involve faults that are intersecting and have low permeability. Both these features are highly relevant in a geologic setting where fault may have complex geometry and reduced permeability compared to the host rock.

3D case
As a final verification, we consider a 3D case with multiple intersecting faults. The setup is based on Case 2 in the benchmark study described in Berre et al. (2020). The domain is a unit cube including a network of 9 faults, whose intersections divide the cubic domain into several subdomains, as illustrated in Fig. 12a. These 3D subdomains are grouped into two regions, where we assigned different permeabilities ,1 and ,2 , both  Figure 11. 2D self-convergence test: (a) mixed-dimensional geometry and (b) convergence of the average error in pressure within the faults.    [148, 282, 384, 814, 1536, 2298, 3456]) and report the average 2 error in pressure along the faults.
Convergence results are shown in Fig. 12b, indicating first order convergence on average. This confirms the consistency of our implementation also for 3D problems with complex fault geometries.

Conclusions
We presented an improved framework to modelling and discretizing flow in generally anisotropic porous media with thin inclusions, within the context of mixed-dimensional partial differential equations. Our model considers a full permeability tensor for the inclusions, resulting in additional terms arising in our formulation as compared to existing local discretizations. We expect our model to be important for modeling of flow in faulted porous media, however the methods proposed herein can be in any case applied to models of fractures, in fact our full-permeability model naturally reduces to the existing models of fracture-matrix flow when the off-diagonal components of the inclusion permeability tensor are set to zero.
We provided numerical examples showing convergence of the method for both 2D and 3D faulted porous media. In particular, we provided numerical evidence that, as opposed to existing local discretizations, our model is capable of simulating the anisotropic behaviour of the faults near damage zone.
We remark that, in the spirit of flux-mortars coupling schemes, our formulation is independent of the discretization methods used to discretize the flow equations in the porous matrix and the faults. However, we only showed results obtained using a multipoint flux finite volume approach. Nevertheless, the formulation also applies to other discretization methods, e.g. mixed finite elements.