An optimization approach for flow simulations in poro-fractured media with complex geometries

A new discretization approach is presented for the simulation of flow in complex poro-fractured media described by means of the Discrete Fracture and Matrix Model. The method is based on the numerical optimization of a properly defined cost-functional and allows to solve the problem without any constraint on mesh generation, thus overcoming one of the main complexities related to efficient and effective simulations in realistic DFMs.


Introduction
The present work deals with the simulation of the flow in the subsoil, modelled by means of the Discrete Fracture and Matrix (DFM) model. According to this model, underground fractures are represented as planar polygons arbitrarily oriented in a three dimensional porous matrix. The flows considered here are governed by the Darcy law in the three dimensional matrix and by an averaged Darcy law on each fracture plane, with suitable matching conditions at fracture-matrix interfaces and at fracture intersections. The quantity of interest is the hydraulic head, given by the sum of the pressure head and elevation. Single phase stationary flow is considered, with the assumption of continuity of the hydraulic head at both fracture-matrix interfaces and at fracture-fracture intersections and no longitudinal flow is allowed along fracture intersections. This is a simplified model with respect to other DFM approaches, described, for example in [28] or, more recently, in [14], but still representative of realistic configurations, characterized, e.g., by highly permeable fractures. The main focus of the present work is on geometrical complexity aspects, proposing a problem formulation and a numerical approach suitable for complex and randomly generated networks. The described approach can however be extended to different flow models and different coupling conditions. The geometrical complexity of DFM models mainly arises from the multi-scale nature of the resulting domains and from the presence of multiple intersecting interfaces, where the solution displays an irregular behavior. DFM models are proposed as an alternative to homogenization techniques [33], dual and multy-porosity models [17], or embedded discrete fracture matrix (EDFM) models [26,29,23], and are characterized by the explicit representation of the underground fractures, dimensionally reduced to planar interfaces into the porous matrix. As a consequence of the random orientation, fractures usually form an intricate system of intersections, with the presence of fractures with very different sizes spanning several orders of magnitude that generate intersections with huge geometrical complexities as, for example, 2D and 3D geometrical objects with very different dimension and objects with enormous aspect ratios. The research on effective numerical tools for DFM simulations is particularly active, see e.g. [3,1,15,32,12,31,4,18]. One of the key aspects is the meshing of the domain, with a mesh conforming to the interfaces, suitable for standard approaches for the imposition of interface conditions. The generation of a conforming mesh for realistic fracture networks might, in fact, result in an impossible task, for the extremely high number of geometrical constraints. The mesh conformity constraint at the interfaces can be relaxed by using extended finite elements as suggested, e.g., by [24,21]. Different approaches are based on the Mimetic Finite Difference method (MFD) [27], as described, for example, in [2,32], or on Hybrid High Order (HHO) methods as proposed by [16], where a partial non-conformity is allowed between the mesh of the porous medium and of the fractures, or also on Discontinuous Galerkin discretizations, as in [4]. Two or multi-point flux approximation based techniques are described in [34,20] and gradient schemes in [15]. Virtual Element (VEM) based discretizations have also been recently investigated to ease the mesh generation process in complex DFMs, as in [6] where the VEM is coupled to the Boundary Element method, and in [22], in [19] for poro-elasticity problems, or in [5] where an arbitrary order mixed VEM formulation is proposed.
This work presents a development of an optimization-based approach, first proposed for Discrete Fracture Networks [10,11,7,13,8] and recently extended to DFM problems in [12]. This approach avoids any mesh conformity requirement for the imposition of interface conditions, which are instead enforced through the minimization of a properly defined cost functional. The computation of the quantities involved in functional definition does not require any constraint on the mesh. Further, the resolution of the optimization problem via a gradient-based scheme allows to de-couple the problems on each fracture and the problem on the porous matrix, thus paving the way for an efficient parallel implementation of the numerical scheme, similarly to what done in [13,9]. The discretization scheme described in [12] relies on the Boundary Element Method for the discretization of the problem on three dimensional matrix blocks, thus requiring the splitting of the original three dimensional domain into sub-domains not crossing the fractures, and thus implying a partial mesh conformity at the fracture-matrix interfaces. Here, the three dimensional domain is not split into sub-domains and Finite Elements are used for the discretization of the matrix, on tetrahedral elements that can arbitrarily cross Fig. 1 Nomenclature exemplification the fractures. Finite elements on triangular meshes are used for the fractures, with elements not conforming to the tetrahedral mesh and also arbitrarily placed with respect to fracture-fracture intersections. The proposed discretization approach thus greatly improves the usability of the method to general DFM geometries, allowing a trivial meshing process of extremely complex domains, thanks to the complete independence of the mesh from all the interfaces.
The structure of the manuscript is the following: Section 2 describes both the classical and the optimization based formulation of the flow problem in a DFM; the following Section 3 describes the derivation of the discrete problem and the proof of its well posedness; Section 4 shows how an equivalent unconstrained optimization problem is derived, and the gradient based scheme used for problem resolution; Section 5 reports some numerical results and finally some conclusions are proposed in Section 6.

Problem description
This section is devoted to a brief description of the problem of interest, referring to [12] for a more detailed exposition and well posedness results. Let us consider a polyhedral block of porous material, denoted as D, crossed by a fracture network Ω given by the union of planar polygonal fractures F i , i = 1, . . . , N F in the threedimensional space, i.e. Ω = N F i=1 F i . We further denote by F the set of all fracture indexes. Fractures might intersect, and fracture intersections, also called traces, are indicated as Sm, m = 1, . . . , N S . We assume, for simplicity, that each trace is given by the intersection of exactly two fractures, such that an injective map σ : F ] can be defined between a trace index and a couple of fracture indexes, as σ(m) = {i, j} being Sm =F i ∩F j . Further, S i is the set of indexes of all the traces on fracture F i and S the set of indexes of all the traces in the network. Let us introduce the domainD = D \Ω, thus given by the original block D without the internal fractures. Calling ∂D the boundary ofD, let us denote by Γ ± i the portion of ∂D that matches fracture F i , for i = 1, . . . , N F , the superscript "+" or "−" referring to one of the two sides of the boundary "around" the fracture (see Figure 1); the unit normal vector to Γ ± i is n ± Γi , always pointing outward fromD. A jump operator is introduced for any sufficiently regular vector function v onD, defined as the jump of v along the normal direction to the faces Similarly, for i = 1, . . . , N F we denote byF i the fracture F i without traces, i.e. F i = F i \ m∈Si Sm, and for each trace Sm, m ∈ S i , for any sufficiently regular vector function w i on F i , the jump of the normal component of w i across trace Sm on F i is denoted as: with S ± m the two sides of the portion of the boundary ofF i lying on Sm and n i Sm the normal unit vector to Sm with a fixed orientation on F i . These jump operators are easily extended to functions defined on the whole 3D domain D and on the whole fractures F i , i = 1, . . . , N F , with the ± superscripts still denoting the two sides of the interface Γ i ≡ F i , ∀i = 1, . . . , N F , or Sm, ∀m = 1, . . . , N S . The portion of ∂D not matching any fracture is split in a Dirichlet part Γ D and a Neumann part Γ N , Γ D ∩ Γ N = ∅, where, for simplicity of exposition, we assume homogeneous Dirichlet and Neumann boundary conditions are enforced. Similarly, the boundary of each fracture ∂F i , i = 1, . . . , N F , is split in a Dirichlet and Neumann part, γ iD and γ iN , respectively. If fracture F i lies in the interior of D, then we set γ iD = ∅, and homogeneous Neumann boundary conditions are prescribed on γ iN ≡ ∂F i . If N F = 1, we assume that |γ 1D | > 0, whereas, if there is more than one fracture in the network, we allow γ iD = ∅ for i = 1, . . . , N F . The problem of the equilibrium distribution of the hydraulic head in D can be then stated in strong formulation as: H i = 0 on γ iD , i = 1, . . . , N F (8) where H D is the hydraulic head inD, H i the hydraulic head on F i , i = 1, . . . , N F and f is a volumetric source term. The operator ∇ represents the three-dimensional gradient inD, ∇ i is the two-dimensional gradient on the plane containing fracture , for x ⊂D is a symmetric positive definite matrix representing the transmissivity of the porous matrix and K i (x) ∈ R 2×2 , x ⊂ F i is a symmetric positive definite matrix representing the tangential transmissivity of the fracture F i on its tangential plane. Finally, n Γ N is the outward unit normal vector to Γ N , and for a given index i = 1, . . . , N F , nγ iN the outward unit normal vector γ iN on the plane of fracture F i . Here, for simplicity, we have considered only the source term on the fractures deriving from the exchange with the porous matrix and homogeneous boundary conditions, but the extension to a more general case is immediate. Conditions (3) and (4) express the continuity of the solution at fracture-matrix interfaces and at fracture intersections, respectively, whereas Equation (5) enforces the balance of fluxes at the traces.
Let us now introduce the following functional spaces: first, on each fracture Also, on each trace Sm, m = 1, . . . , N S we set the spaces U m = H − 1 2 (Sm) and H m = H 1 2 (Sm). We introduce the following variables: thus representing a sort of internal Robin boundary condition on the traces; and, thus again being a linear combination of the jump of the co-normal derivative of H D across interface Γ i and the trace of H D on Γ i , and V i the dual of V i . We remark that, as H D ∈ H 1 Ω (D) the hydraulic head is continuous across interfaces We also define the bilinear forms: Then, problem (3)-(9) can be written in weak formulation as: find being (v, w) ω the scalar product in L 2 (ω). The coupling conditions in weak form are given by: for all i = 1, . . . , N F , and for all m = 1, . . . , Parameters α > 0 and β > 0 ensure stability of the problems written independently on each fracture and on the three dimensional domain. This is required to obtain a discrete formulation suitable for parallel computing. If N F = 1, or equivalently S = ∅, well posedness of the problem on the unique fracture F 1 is guaranteed by the assumption on |γ 1D | > 0.
Problem (12)-(16) is well posed. To show this, let us introduce the function space H 1 Ω + (D) defined as: and thus incorporating the matching conditions at the interfaces. Let us then write the following problem: Problem (17) is well posed, as it can be easily seen that H 1 Ω + (D) is an Hilbert space with the scalar product, [12]: Problem (12)-(16) is equivalent to problem (17); indeed, recalling that, (14)- (15) are satisfied by construction. Moreover summing (13) for i = 1, . . . , N F and (12), using (16) and the definition of U m i and Q i , for i = 1, . . . , N F , m ∈ S i , we get (17). We propose a reformulation of problem (12)- (16) well suited for discretization on non conforming meshes and parallel computing, based on a PDE constrained optimization approach. To this end, we introduce a Table 1 Labels used for the dimension of discrete variables cost functional expressing the error in the fulfilment of the interface conditions as continuity and flux conservation: The solution to problem (12)-(16) is obtained as the minimum of functional J(H D , H F , U S ) constrained by the PDE equations on the 3D domain and on the fractures: constrained by (12) and by (13) 3 Discrete formulation The PDE constrained optimization formulation is specifically designed to allow for an easy discretization of the problem using non conforming meshes and to obtain a discrete problem suitable for effective resolution using parallel computing resources. The imposition of the interface constraints expressed by equations (14)-(16) with a standard approach requires some sort of mesh conformity at the interfaces: either a perfect matching of the nodes on the meshes to enforce conditions by means of degrees of freedom equality constraints, or the weaker condition of alignment of mesh edges with the interfaces, to use mortaring techniques. In contrast, the imposition of interface conditions through the functional only requires the computation of integrals on the traces, as shown below, and thus meshes can be arbitrarily placed with respect to the interfaces, see Figure 2 for an example of non-conforming meshes in the rock matrix and on the fractures. Further, the minimization process allows to decouple the problems on the fractures and on the three-dimensional domain, for parallel computing.
The discretization strategy proposed in this work is based on the use of standard finite elements on tetrahedra for the three dimensional domain and finite elements on triangles for the fractures. Let us then denote by T D δ D the tetrahedral mesh on D, characterized by a mesh parameter δ D , by T i δ F i a triangular mesh on F i , i = 1, . . . , N F , with mesh parameter δ Fi , and by T i δ Γ i a possibly different triangular mesh on F i , with mesh parameter δ Γi . We further introduce a discretization of the one-dimensional traces, different on each fracture, denoted by T i δ Sm,i , with mesh parameter δ Sm,i , i = 1, . . . , N F , m ∈ S i . We denote by h D the finite dimensional ap- i,k , with N m ui dofs and m i,k a basis function. Tables 1-2 summarize the labels used for the dimensions of the discrete variables, the name used to denote the basis functions and the notation used, in the following, for the matrices collecting integrals of these basis functions. We build arrays of dofs by collecting column-wise the dofs of each discrete function and with abuse of notation we denote the dof array with the same symbol of the corresponding function, thus having arrays We define arrays u m , m = 1, . . . , N S , as u m = [(u m i ) T (u m j ) T ] T for i, j = σ(m) with i < j, and we further collect column wise arrays h i , q i , u m i and u m forming: where m 1 , . . . , m Si are the indexes in S i ordered increasingly.
The discrete version of functional J is the following: obtained replacing the discretized variables and using L 2 norms. The discrete functional can be written in matrix form, computing the integrals of the basis functions and collecting the values into matrices. Considering the first norm in J , we have: and we can define three matrices as follows, for each The computation of matrix G i DF is not straightforward, as the two involved variables are defined on different meshes. In particular, the intersection of the three dimensional tetrahedral mesh with the fracture plane needs to be computed. This operation defines a polygonal tessellation of F i which is then sub-triangulated, thus generating a triangular interface mesh. This sub-triangulation process can be performed without any mesh quality requirement, as the resulting mesh is used only for quadrature purposes. The computation of the elements in G i DF is finally performed first computing the intersection of the elements of the interface mesh with the elements in T i δ F i , and subsequently the required integral on the intersection region. Element neighbourhood information is used to efficiently perform the task. The computation of the interface mesh is a quite complex and expensive task. Also in this case element neighbourhood information is used for efficiency, and further can be performed independently fracture by fracture and thus in parallel, which is of paramount importance for the applicability of the method to complex geometries.
We can proceed similarly with the remaining terms of the functional J ; to this end, for m = 1, . . . , N S , and all the possible couples of indexes i, j such that Fig. 2 Polygonal tessellation on a sample fracture given by the intersection of a tetrahedral mesh with the fracture plane (blue in the right panel) overlapped with the fracture triangular mesh (red in the right panel) such that We can collect these matrices, defined locally at the various interfaces into global matrices to derive a compact form of the functional. Let us define matrix For all i = 1, . . . , N F , let us assemble matrices for increasing values of m ∈ S i and j, k = σ(m), j < k, i.e.: The matrix formulation of J then reads as: We can re-write also the discrete constraint equations in matrix form. We follow a standard procedure and we define matrix where the integral on F i is performed on the interface mesh, generated by the intersection of the tetrahedral mesh with each fracture.
ii for all m ∈ S i . These matrices are used for the definition of matrix B ∈ R N h F ×Nu , defined grouping column-wise matrices B i , for i = 1, . . . , N F . Matrix D ∈ R N h F × Nq is built as follows: with integrals computed on the intersection of the mesh T i where integrals are computed intersecting the mesh T i δ Γ i for variable q i with the triangulated interface mesh given by the intersection between the tetrahedral mesh with the fracture F i .
, the discrete formulation of the constrained minimization process is: being b D ∈ R N h D the array resulting from the forcing term.
Let us now introduce the following matrices: and let us collect column-wise variables q, u into variable w, then optimality conditions for problem (22)-(24) are given by the following linear system: Well posedness of problem (22)-(24) derives from non singularity of the saddle point matrix M.
Lemma 1 Let matrices A, B be defined as in (25). Let A be full rank, let L = [A −B], and let Z be a matrix obtained collecting row-wise column vectors z k , k = 1, . . . , Nw, Nw = Nu + Nq forming a basis of ker(L), then matrix Z T GZ is positive definite.
Proof We start observing that matrix A is full rank as both matrices A D and A F are full rank under the assumption that α, β > 0. Then dim(ker(L)) = Nw. To construct a basis of ker(L), let us take e k , the k-th vector of the canonical basis of R Nw , and let us set z k = (A −1 Be k , e k ). According to the index k, e k might correspond to a non-null function q i for some i = 1, . . . , N F or a non-null function u m i for some i = 1, . . . , N F , m ∈ S i . In both cases will show that z T k Gz k > 0.

Non-null function q i
Let us start considering the first case and in particular let us assume that q i = ϕ i,j for a certain index j = 1, . . . , Nq i . Let us consider two different scenarios: the case N F = 1, i.e. a porous medium with a single fracture and the more general case N F > 1.
-N F = 1 As q 1 = 0, for the non singularity of A D , it is h D = 0, and in particular it can be either h D = const or h D = const.
• h D = const = 0 As we assumed homogeneous Dirichlet conditions, h D = const = 0 is possible only if Γ D = ∅. In virtue of definition (11), for the consistency and conformity of the method, we have βh D|F1 = q 1 , which is possible only if q 1 is constant on F 1 . By (13), being |γ 1D | > 0, it is: and thus h 1 = 0, and in particular h D|F1 − h 1 L 2 (F1) > 0.
• h D = const Let us set s 1 := q 1 − βh D|F1 , and let us consider equations (12) and (13), that become: where the source term s 1 appears with opposite sign. We thus have If h D = const = 0, proceeding similarly to the case N F = 1, we have βh D|Fi = q i = 0 and h D|Fi − h i L 2 (Fi) > 0. If instead h D = const, we proceed in the following way: since q i = 0 we have h i = 0, whereas it is h j = 0, for all j = 1, . . . , N F , j = i. Choosing, in particular, one index j such that fracture F i and F j intersect in a trace Sm, we have Let us now consider u m i = m i,j for some i = 1, . . . , N F , m ∈ S i , and for an index j = 1, . . . , N m ui , depending on the value of k. Also in this case it can be easily shown that we have h i = 0, whereas we have h D = 0 and hp = 0 for all p = 1, . . . , N F , p = i, thus having again a non null functional value and thus z T k Gz k > 0. Being G positive semi-definite by definition, it is x T Gx ≥ 0 and x T Gx = 0 if and only if x ∈ ker(G), [25] and being z T k Gz k > 0, z k ∈ ker (G) for z = 1, . . . , Nw. The space Z = span{z 1 , . . . , z Nw } is thus a subspace of Im(G), and each vector y ∈ Z can be written as y = Zv, for a vector v ∈ R Nw , v = 0. Then v T Z T GZv > 0. The proof follows from Lemma 1 applying a classical argument of quadratic programming (see Theorem 16.2 in [30]).

Unconstrained optimization problem
We can proceed formally, replacing the constraint equations into the functional, to obtain an unconstrained minimization problem. We have h = A −1 (Bw + b), from which we obtain: The unconstrained minimization problem then reads min w w T Gw + 2g T w (27) or equivalently Gw + g = 0. Matrix G is symmetric positive definite, given the equivalence of (27) with (22)- (24). The unconstrained minimization problem can thus be solved with a gradient based iterative method, such as the conjugate gradient method. The steps of the method are as follows: The computation of quantity y k = Gd k , at each step k can be performed as follows: k . If β = 0, which is possible as long as there is a non empty portion of the Dirichlet boundary for the three dimensional domain, i.e. |Γ D | > 0, the computation ofh k , λ k at each step can be performed independently and in parallel on each fracture and on the three dimensional domain, thus easily allowing to use parallel computing resources for efficient resolution of the scheme, thanks to the block diagonal structure of A and A F . If β > 0 then problems on the fractures can be decoupled from the problem in the bulk domain as follows: at step k > 0, beinḡ and similarly for λ k .

Numerical results
In this section we provide some numerical results in order to show the applicability of the present approach to flow simulations in porous media crossed by arbitrarily complex networks of fractures. All the simulations are performed using linear Lagrangian finite elements on T D δ D for h D , linear Lagrangian finite elements on T i δ F i for h i , piecewise constant basis functions on T i δ Γ i for q i and piece-wise constant basis functions on T i δ Sm,i for u m i on each trace Sm, on each fracture F i , i = 1, . . . , N F , m ∈ S i .

Problems with known solution
We first propose two simple problems with known analytical solutions, labelled Problem 1 and Problem 2, having the same domain and type of boundary conditions. A cubic domain with unitary edge length is considered; the bottom face is on the plane z = − 1 2 with respect to a reference system Oxyz, and the cube is crossed by a single fracture F 1 placed on the plane z = 0, see Figure 3, left. The problems are set as follows: with f = −1 for Problem 1 and f = 0 for Problem 2, β = 1, K D = K 1 = 1 for both problems. Dirichlet boundary conditions are set on cube faces on planes z = − 1 2 , z = 1 2 , x = 0, x = 1, Neumann boundary conditions on cube faces on planes y = 0 and y = 1. Boundary conditions on fracture edges are prescribed accordingly to the boundary conditions on cube edges. Dirichlet and Neumann boundary conditions are derived from the analytical solution, which is h = 1 4 (x 2 +y 2 )+ 1 2 |z| for Problem 1 and h = 1 2 (x 2 −y 2 )+z for Problem 2. The two problems here considered also share the same meshes. In Figures 3, right and 4, we display the colormap of the solution of Problem 1. The mesh for the three dimensional domain is non conforming with the fracture plane and independent from the mesh on the fracture, as shown in Figure 3, right. In Figure 5 we report the behaviour of the error with respect to the mesh size both in L 2 and in H 1 norm for Problem 1. The three dimensional mesh parameter ranges between 0.02 and 4 × 10 −5 , the mesh on the fracture between 0.3 and 5×10 −3 . Due to the non conformity of the mesh and to the irregular behaviour of the solution across the interface, sub-optimal convergence trends are obtained. The obtained slopes for the error are compatible with the bounded regularity of the solution h ∈ H 2 (D). Optimal convergence curves are however recovered if the solution across the interface is smooth. In fact, if we consider Problem 2, having a smooth solution, optimal convergence trends are recovered, as reported in Figure 6.

DFN problem
In the second example, a more complex and realistic DFN is considered, embedded in a cubic domain, with barycentre in the origin of a reference system Oxyz and edge length equal to two, as shown in Figure 7, on the left. The embedded DFN consists of 20 randomly placed fractures, forming 62 traces, with a number of traces per fracture ranging between 4 and 10. Traces intersect, forming angles as narrow as 11.5 degrees, whereas the minimum angle between the normals of couples of intersecting fractures is 17.3 degrees. A unitary pressure Dirichlet boundary condition is imposed on fracture edges lying on the planes z = 1, and a zero pressure Dirichlet boundary condition is imposed on the cube face on the plane z = 0, all other fracture edges and cube faces being insulated. An inflow is thus obtained through some fracture edges, and outflow occurs through the bottom face of the cube. Figure 7, on the right, shows the computed solution on the 3D domain and on the fractures, through a section of the three dimensional domain, along with the used mesh, characterized by a mesh parameter equal to δ D = 0.01 for the tetrahedral mesh and to δ Fi = 0.05, i = 1, . . . , 20 for the triangular mesh on all the fractures. On this mesh, the total number of unknowns, N h + Nq + Nu is 5541, and the minimization problem is solved using 90 iterations to reach a relative residual of 10 −8 , resulting in a functional value of 0.069, on the considered mesh. Considering a refined mesh, with mesh parameter 1.25 × 10 −4 for the tetrahedral mesh and to 1.25 × 10 −2 for the triangular mesh, the total number of unknowns rises to 28686 and the number of iterations to reach the same relative residual is 125 and the functional value is reduced to 0.057. We remark that the minimum of the discrete functional is greater than zero, as a consequence of the non conformity of the mesh. These results show the viability of the proposed approach in dealing with complex domains.

Conclusions
A new discretization strategy for the simulation of the flow in arbitrarily complex DFM geometries has been presented and validated. The method is based on standard finite element discretizations for both the three dimensional domain and the fractures, and the meshing can be performed independently on each geometrical entity, thus actually overcoming any mesh related issue for DFM simulations. The resulting discrete problem is well posed and can be efficiently solved via a gradient scheme. The proposed numerical tests validate the method and show its applicability to realistic DFM configurations. Although the proposed method can be easily implemented for parallel solution, optimal parallel solver and suitable well balanced partitioning strategies, yielding to efficient parallel solvers, should be investigated but are out of the scope of the present work.