A dimensionally-reduced fracture flow model for poroelastic media with fluid entry resistance and fluid slip

We develop a model which couples the ﬂow in a discrete fracture to a deformable porous medium. To account for the discrete representation of the fracture, a dimensionally-reduced ﬂuid ﬂow model is proposed. The ﬂuid ﬂow model incorporates both a reduced permeability of the fracture walls due to the skin effect, and a slip of ﬂuid ﬂowing along the permeable fracture walls. Biot’s model for poroelastic media is coupled to a fracture ﬂow model based on a thin-ﬁlm approximation of the compressible Navier-Stokes equations. The fracture ﬂow model incorporates a ﬂuid entry resistance parameter to relate the leak-off through the fracture walls to a pressure jump across the fracture walls, and the Beavers-Joseph-Saffman slip rate coeﬃcient to represent the ﬂuid slip along the fracture walls. The numerical model is based on a thermodynamic framework in which all energy storage and dissipative mechanisms in the problem are identiﬁed, including the mechanisms related to the interface effects. The thermodynamic framework is employed to solve the nonlinear coupled problem up to a speciﬁed energy range through a Picard iteration technique and to study the model and its results. Studies are presented for a range of ﬂuid entry resistance parameters and Beavers-Joseph-Saffman slip rate coeﬃcients, showing the capability of the model to simulate skin and slip effects dimensionally-reduced


Introduction
The study of fluid flow through fractured porous media is relevant to a wide array of applications. A prominent application lies in the field of geomechanics, where both naturally and artificially created fractures exist in subsurface reservoirs. The (computational) modeling of fractured porous media is very challenging on account of the multi-physics and multiscale character of the relevant physical phenomena. In particular, the typically large separation in length scales between the fracture aperture and the scale of the medium in which the fracture is embedded requires a dedicated modeling and simulation approach, especially with regard to the fluid flow inside the fracture.
Over the past decades, a wide range of simulation techniques has been developed to model the flow through fractured porous media. Recent overviews of the state of the art can be found in, e.g., Refs. [1][2][3]. Two categories of models for fractured reservoirs can be distinguished, namely, implicit fracture representation models, in which the fractures are considered as altered properties of the porous continuum, and explicit fracture representation models, in which the fracture and the medium in which it is embedded are treated as separate geometric entities.
Building on Bergkamp et al. [4], this contribution focuses specifically on the interaction between the fluid flow in the fracture and the surrounding porous medium. To model interface phenomena, it is natural to consider an explicit fracture representation. In Bergkamp et al. [4] we present a model and a (staggered) numerical solution strategy that fully couples a free flow and a deformable porous medium over a shared boundary. The free flow is treated as an incompressible Stokes flow and the behavior of the porous medium is described by Biot's equations with Darcy flow. Over the interface, we consider not only conservation of fluid mass and momentum, but also the (partial) slip encountered by a fluid as it flows along a permeable medium, as described by the Beavers-Joseph-Saffman condition [5,6], and the presence of a skin layer of relatively low permeability. To model this skin layer, we relate the fluid flux to a jump in pressure over the interface through Showalter's fluid entry resistance parameter [7,8]. By scaling the fluid entry resistance parameter, we can use the pressure jump to numerically represent a steep pressure gradient over the thin layer of reduced permeability.
In Ref. [4], the flow model inside the fracture is fully-dimensional, i.e., the model is defined on a volumetric domain. Since the mesh resolution in this fully-dimensional model is dictated by the aperture of the fracture, upscaling to situations in which the porous medium domain is many orders of magnitude larger than the fracture aperture is computationally impractical. To enable upscaling, in this contribution we consider a discrete fracture representation, where a fracture is represented by a surface (in three dimensions) that is embedded in a volumetric porous medium. We develop a dimensionally-reduced fluid flow model, i.e., defined on the embedded surface (lower-dimensional manifold), which retains the interface conditions of Bergkamp et al. [4]. In the context of numerical implementation, this work builds on the Enhanced Local Pressure (ELP) method [9,10], which is an Extended Finite Element Method (X-FEM) [11,12] for poromechanical continua, with the ability to represent a jump in fluid pressure across the interface between a dimensionally-reduced fracture and a porous medium.
In this manuscript, we present two key innovations in the computational modeling of fractured porous media: • A dimensionally-reduced flow model is derived for the flow inside the fracture. The commonly used Reynolds equation [13,14] does not take into account fluid slip and leak-off as represented by the interface relations in Ref. [4]. We derive an improved fluid flux relation based on the thin-film approximation of the compressible Navier-Stokes equations [14], that does take these effects into account. Various contributions have been made to address the dimensional-reduction problem, see, e.g., [15][16][17][18][19]. A fracture-flow model that accounts for both fluid slip and fluid entry resistance effects in the interaction with a poroelastic medium has, to the best of our knowledge, not been considered.
• As part of the numerical model, we present a thermodynamic framework for the fractured poroelastic medium model. This thermodynamic framework, which builds on the established work for Biot's poroelasticity model (e.g., [20]), identifies all energy storage and dissipative mechanisms in the problem, including these directly related to the interface conditions of Ref. [4]. The thermodynamic analysis provides detailed insights on the influence of the model parameters, specifically the fluid entry resistance and the slip parameter, on the different dissipative mechanisms. In this work, the thermodynamic framework is also used to develop a Picard iteration technique that solves the nonlinear coupled problem up to a specified, physics-based, energy rate threshold. Although not considered in this work, the identification and formalization of all dissipative mechanisms provides a basis for an extension to fracture propagation using established energy-based criteria (e.g., [21,22]).
These modeling and simulation developments endow our fracture simulation framework with the capability of incorporating fluid slip and surface permeability effects at the scale of reservoir fractures. A detailed numerical study of this novel capability is considered herein. This manuscript is outlined as follows. In Section 2, we commence with the derivation of the porous medium model, the dimensionally-reduced fracture flow model, and the interface conditions employed to couple these models. The thermodynamic framework is then presented in Section 3. Section 4 presents the employed numerical solution technique, consisting of an incremental-iterative solution procedure and a spatial discretization technique based on the ELP method. In Section 5, we demonstrate the proposed method using a range of numerical simulations. First, a verification study of the simulation framework is presented. Subsequently, the influence of the fluid entry resistance parameter and the fluid slip parameter, is studied. Conclusions are presented in Section 6.

Fractured poroelastic medium model
We consider the interaction between a fluid-saturated poroelastic medium and a stationary, i.e., non-propagating, fluidfilled fracture, as depicted in Fig. 1. The poroelastic medium is represented by the closed domain ⊂ R d , for d = 2, 3. The domain has an exterior boundary , with unit normal vector n directed out of . The fluid-filled fracture is embedded in the poroelastic medium and is considered to have a thickness which is relatively small compared to its other spatial dimensions. Due to its high aspect ratio, we consider the fracture to be a (d − 1)-dimensional manifold (i.e., a surface in three dimensions, or a curve in two dimensions). We consider the fracture to be represented by the smooth interface c . The two sides of this interface, + c and − c , represent the fracture walls on either side of the fluid domain, as illustrated in the zoom in Fig. 1. Together, these two sides form an interior boundary of the poroelastic medium (i.e., the total contact surface between the fluid in the fracture and the porous medium), ± c = + c ∪ − c . Throughout this manuscript, we use the superscripts + and − to label the two sides of the interface, and the superscript ± for the union of these sides. Consistent with the external boundary, we define the unit normal vector n on ± c to be directed out of the porous medium . We assign the unique normal vector n c = n − to the interface c , such that the interface normal vector points from − c to + c . The boundary of the interface c (i.e., a curve in three dimensions, or points in two dimensions) is formed by the (d − 2)-dimensional manifold ∂ c , which represents the fracture front.
The porous medium is assumed to be elastic and the continuous pore space is assumed to be fully fluid-saturated. The combined behavior of the porous structure and the pore fluid is governed by Biot's theory with Darcy flow, as reviewed in Section 2.1. The choice to describe the fracture as a smooth (d −1)-dimensional manifold has implications for the description of the fluid flow within the fracture. Specifically, the lack of spatial dimension in the thickness-direction is limiting, in the sense that we have to make assumptions regarding the flow profile across the manifold. This dimensional reduction of the fluid flow model presented in Bergkamp et al. [4] is discussed in Section 2.2. The coupling conditions of Ref. [4] are then adapted to the dimensionally-reduced fracture setting in Section 2.3. An overview of the complete model is presented in Box 1.

The poroelastic continuum
The domain shown in Fig. 1 is occupied by material points x ∈ . We treat the poroelastic medium occupying as a continuum, i.e., we only consider macroscopically averaged quantities. The displacement of the solid structure of the porous medium is denoted by u(x, t) at x ∈ and the corresponding velocity is denoted by u(x, t). The velocity of the fluid occupying the continuous pore space of the porous structure is denoted by v(x, t) and the hydrostatic pore pressure is denoted by p(x, t). The specific discharge, which is a measure for the fluid flow through the porous medium, is represented by q = φ(v −u), where φ is the porosity of the porous medium, and (v −u) is the velocity of the pore fluid relative to the porous structure.
The behavior of the porous structure is assumed to be isotropic, only small deformations and deformation gradients occur, and the stress-strain relations are assumed to be linear and reversible. In accordance with these assumptions, we consider the infinitesimal strain tensor ε = ∇ s u = 1 2 ∇u + (∇u) T . Furthermore, we assume the porosity φ of the medium not to be affected by the deformation field. The porous structure is assumed to be saturated with a single-phase viscous fluid, whose behavior follows Darcy's law [23], where k is the permeability of the porous medium and η is the dynamic viscosity of the fluid. Following the above assumptions, the behavior of the poroelastic medium can be described by Biot's theory [24].
Assuming that there is no mass transfer between the porous structure and the pore fluid, Biot's theory states that mass conservation of the porous medium is described by the storage equation [25], where α is the Biot-Willis coefficient [26] and M is the compressibility modulus, as elaborated in Remark 1. The mass balance is completed with the boundary conditions q · n =q on q , where q is the prescribed filtration velocity on the external boundary q and p is the prescribed pore pressure on the external boundary p , with q ∪ p = and q ∩ p = ∅.
We neglect inertia, convection, and body forces, including gravity. Furthermore, we consider the consolidation process, i.e., the settlement of a porous medium subjected to a load, to be isothermal. The momentum balance then reduces to, ∇ · σ = 0, (4) where σ is the total stress in the porous medium, σ = σ e − αp I , (5) with σ e the effective stress, as proposed by Terzaghi [27,28]. The effective stress is given by Hooke's law, σ e = 2με + λtr(ε)I, (6) where μ and λ are the Lamé parameters, which are related to the Young's modulus, E, and the Poisson ratio, ν, of the drained porous structure, and tr(ε) = ∇ · u. The momentum balance is completed with the boundary conditions where ū is the prescribed displacement of the porous structure on the external boundary u and t is the prescribed traction on the external boundary t , with t ∪ u = and t ∩ u = ∅. By considering directions normal and tangential to a boundary separately, combinations of (7a) and (7b) may also be applied.

Remark 1 (Material parameters of the porous medium)
. The Biot-Willis coefficient [26], α, can be expressed in terms of the bulk modulus of the drained specimen, K , and the bulk modulus of the solid constituent, K s , as The compressibility modulus M is defined as with K f the bulk modulus of the pore fluid.

Dimensionally-reduced fracture flow model
We embed a stationary fluid-filled fracture, represented by the smooth interface c , in the fluid-saturated poroelastic medium , as shown in Fig. 1. The smooth interface c is occupied by material points x ∈ c . Across the interface, both the displacement of the solid structure u and the pore pressure p are considered to be discontinuous. On the interface, we define the jump in the displacement field as u = u + − u − , and the average displacement as {u} = 1 2 (u − + u + ). To simplify notation, in the remainder we will express the displacement jump as = u , with normal component n = · n c , and non-normal components s = − n n c . For the pore pressure we similarly define the jump as p = p + − p − , and the average as {p} = 1 2 (p − + p + ). The behavior of the fluid in the fracture is characterized by the fluid pressure p c (x, t) for x ∈ c . The fluid pressure in the fracture is defined separately from the pore pressure, to allow for a discontinuous pressure over the fracture walls, similar to Refs. [9,10]. This enables modeling of, e.g., the skin effect, as elaborated in Bergkamp et al. [4].
In this section, we derive the governing equations for the flow in the dimensionally-reduced fracture. We commence by treating the domain on which the fluid flow model is defined as fully-dimensional (d-dimensional), with the restriction that its thickness is considered to be relatively small compared to its other (d − 1) dimensions. The fully-dimensional representation of the fracture flow domain is depicted on the right in Fig. 1, with the coordinate ξ associated with n c , i.e., aligned with the thickness direction of the fracture. Within this fully-dimensional representation of the fracture, we consider the thin-film approximation of the compressible Navier-Stokes equations [14]. The governing equations are subsequently integrated over the thickness of the fracture to yield the governing equations for the dimensionally-reduced fracture flow. The case considered here is non-standard, in the sense that we substitute the Beavers-Joseph-Saffman condition [5,6] in our derivation, allowing us to prescribe the fluid-slip encountered by the fracture fluid as it flows along the porous fracture walls.
We first consider the thin-film approximation of the compressible Navier-Stokes equations. Assuming the thickness of the fracture to be substantially smaller than the other dimensions, the thin-film approximation of the mass and momentum balance read (see, e.g., Ref. [14] for a detailed derivation) where ρ is the density of the fluid, ∂ n (·) = ∇(·) · n c is the gradient component normal to the fracture, ∇ s (·) = ∇(·) − ∂ n (·) n c is the surface gradient (see, e.g., Ref. [29]), v c is the velocity of the fluid in the fracture (with normal component v c,n = v c · n c and non-normal component v c,s = v c − v c,n n c ), and ς [m/s] is the injection rate per unit surface area. The pressure in the fully-dimensional fracture is considered to be constant in the direction normal to the fracture midplane and is therefore equal to the pressure in the dimensionally-reduced fracture. Hence, we denote the pressure by p c in both situations, observing that ∂ n p c = 0.
The thin-film balance laws (8) are complemented by an equation of state to relate the density to the pressure. We assume the density to be constant, i.e., ρ = ρ 0 , and the rate of density to be proportional to the rate of pressure as ρ = ρ 0 K −1 fṗ c (see Remark 2). Note that, due to rapid pressure changes, the mass density rate can be non-negligible, despite the variations in the mass density itself being very small. The mass balance (8a) can then be simplified to Since the pressure is assumed to be constant through the thickness of the fracture, integration over the thickness yields where the flux through the fracture is defined as and where the total leak-off into the reservoir, = q · n c , and the fracture opening rate, ˙ n = u · n c , follow from the integration of the last term on the left-hand side of (9): The integrated form of the fracture mass balance is discussed in Remark 3.
To attain the non-normal component of the fluid flow velocity through the fracture, v c,s , we combine the momentum balance (8b) with the Beavers-Joseph-Saffman condition. The Beavers-Joseph-Saffman condition is valid at the two sides of the interface (the fracture walls) [5,6,30], where it relates the traction acting on the fluid in the direction tangent to the fracture to the tangential fluid slip as Following the thin-film assumptions, the traction acting on the fluid can be expressed as such that Integration of (8b) using these boundary conditions yields a quadratic flow profile within the fracture as The flux through the fracture follows by integration over the thickness of the fracture as The traction acting on the fluid along the fracture walls follows by back-substitution of (16) into (13) and (14) as The flow profile (16) is illustrated in Fig. 2 (solid lines). Reference results (markers) are obtained using a steady state finite element approximation of the Navier-Stokes equations. For this reference result, the flow of water through a channel of length 10 cm and height 1 mm, with a pressure drop of 10 kPa/m, is considered. Fig. 2a illustrates the influence of the slip parameter, β, and conveys that in the limiting case of β → ∞ and no wall motion, a non-slipping Poiseuille flow is obtained (see Remark 4). Note that Fig. 2a presents the flow profile in dimensionless form. The dimensionless form is derived in Remark 5, where the dimensionless slip parameter β = β n / √ k is introduced. It is observed that fluid slip can be a result of a decrease in the slip parameter, β, a decrease in the fracture aperture, n , or an increase in the permeability of the porous medium. Fig. 2b considers the no slip case, but with a shear deformation of the fracture walls. This figure conveys that the flow relation (16) correctly accounts for the fracture wall motion. The governing equations for the dimensionally-reduced fracture flow are summarized in Box 1. Note that, to close the system of equations for the pressure in the fracture, we assume the flux through the fracture front (tips in the twodimensional case) to vanish, i.e., Q c,s = 0 on ∂ c .

Remark 2 (The equation of state).
The Tait equation is a generally accepted compressibility model for liquids [31]. In its integrated form, the Tait equation is given by (19) where K f is the bulk modulus and K f its variation with pressure, both measured in the reference state (ρ 0 , p 0 ).
Linearization of (19) around the reference state yields For the operating conditions considered in this work, the pressure difference can be assumed to be much smaller than the bulk modulus of the fluid, i.e., p − p 0 K f . Under these conditions, it is valid to only consider the leading terms in (20): To exemplify the operation conditions considered in this work, we consider water at 100 • F compressed at 100 bar. The bulk modulus, K f , of water is then approximately 2.3 GPa (23 kbar), and the (dimensionless) bulk pressure gradient, K f , is approximately 3.4 (see Ref. [32]). Since the pressure is orders of magnitude smaller than the bulk modulus, the abovementioned assumption holds. The relative linearization error for the density in (21) at these conditions is 0.4%. It is also noted that in general, the bulk modulus and its derivative are temperature-dependent. Temperature-dependence of the density can, however, reasonably be neglected in the operating regime considered in this work (see, e.g., Ref. [33]).

Remark 3 (Integrated fracture flow mass balance).
Integrated over the fracture, the mass balance (10) readṡ where use has been made of the Q c,s = 0 fracture front condition. The components in (22) are defined aṡ where V comp , V leak−off and V aperture are the volume rates associated with the compressibility of the fluid, the leak-off, and the fracture aperture, respectively. Moreover, I represents the volumetric inflow rate. The net fracture volume rate, which we define as V net = I −V comp −V leak−off −V aperture , equates to zero.
where v pois c,s is the thickness-averaged Poiseuille flow velocity in equation (24). The fracture flow profile (16) can then be expressed aŝ This dimensionless flow profile is illustrated in Fig. 2.

The interface coupling conditions
The fluid flow in the fracture, c , is coupled to the porous medium, . We employ the coupling conditions of Ref. [4], which are here adapted to the dimensionally-reduced interface setting described above.
Conservation of mass Over the interface, conservation of mass is enforced. In line with the assumptions made in Section 2.2, the density of the fluid is assumed not to vary over the interface between the fracture fluid and the porous medium, such that the conservation of mass can be expressed as v + (28) Using the identities in Remark 6, these conditions can alternatively be expressed using the average and jump operators as where v c,n = v c · n c . Using the above-defined definitions for the fracture aperture and total leak-off, we obtain Since it holds that a positive leak-off represents a total flow rate from the fracture into the porous medium.
Conservation of momentum Following Ref. [4], the momentum balance across the interface is represented by the equilibrium of traction where the traction vectors acting on the fluid are given by equation (18). In terms of the average and jump operators, the momentum balance can be expressed as where, using equation (18), it follows that The fluid-entry resistance condition In our model, we incorporate the option to model the presence of a skin layer of reduced permeability at the interface (see, e.g., Ref. [34]). Following Ref. [7], in the continuum setting we model the influence of this layer by relating the flux across the fracture walls to the pressure difference between the porous medium and the flow in the fracture as where the parameter γ [kg/m 2 s] is referred to as the fluid entry resistance. Using the identities in Remark 6, these conditions can be written as where use has been made of the assumption that the pressure inside the fracture is constant across its thickness, i.e., p c = 0. Note that the total leak-off (31) can be expressed as Slip condition In contrast to Ref. [4], the fluid slip condition is here not imposed as a separate interface condition. Instead, the Beavers-Joseph-Saffman condition (13) is used for the derivation of the fracture flow profile discussed above. See, e.g., Refs. [5,30,6], for discussions on the validity of the Beavers-Joseph-Saffman condition. The validity of the finite element implementation of this model is also considered in the simulations in Ref. [4]. It is noted that, in the absence of deformation rates of the porous medium, the Beavers-Joseph-Saffman condition resembles the Navier slip condition [35].
All interface conditions are summarized in Box 1.

Remark 6 (The average and jump operators).
Various identities hold for the jump and average operators. Consider a function, a, which is discontinuous over the interface. The value of a on either side of the interface c , i.e., − c and + c , can be expressed as Porous medium mass balance: Darcy's law: Porous medium momentum balance: Terzaghi effective stress / Hooke's law: Fracture flow (Reynolds flow): Overview of the model for the poroelastic medium with a dimensionally-reduced fracture.
respectively. For the product of two discontinuous fields, a and b, the jump operator evaluates to Similarly, the average of this product can be expressed as It is noted that since it by definition holds that n − = n c and n + = −n c , it follows that n = −2n c and {n} = 0. Consequently, it follows that an = −2{a}n c and {an} = − 1 2 a n c .

Thermodynamic framework for the fractured porous medium
We here present a thermodynamic framework for the model presented above. This framework will be used to derive an energy-based convergence criterion for the iterative solver considered in Section 5, and can serve as a basis for an extension to fracture propagation based on the energy release rate.
We restrict our analysis to isothermal consolidation processes, in which case it holds that [20] U +Ḟ =Ẇ +Ẇ irr , (42) which states that the sum of the internal energy rate, U , and the dissipation rate (due to friction), Ḟ , is equal to the sum of the mechanical power Ẇ and the power associated with the irreversible processes, Ẇ irr . The energy rate balance (42) can alternatively be expressed as that the rate of total energy of the system is equal to zero: To demonstrate that the model presented above satisfies this energy rate balance in the case of a stationary fracture, we consider the rate of internal energy of the system to be decomposed aṡ where, for an isotropic linear elastic fluid-saturated porous medium, the internal energy density, U porous , can be expressed in terms of the strain, ε, and the fluid content, ζ , as [20,36] In this expression the undrained Lamé parameter, λ u , is related to the drained Lamé parameter, the Biot-Willis coefficient and the compressibility modulus as λ u = λ + α 2 M. From the internal energy it directly follows that where (46a) directly results in Terzaghi's effective stress relation in combination with Hooke's law, and where (46b) expresses the fluid content in terms of the strain and pore pressure as Under the assumption that the pressure is much smaller than the bulk modulus of the fluid, i.e., p c K f (see Remark 2), the internal energy density rate for the flow in the fracture is related to the pressure in the fluid bẏ The dissipation functional associated with the frictional fluid can be decomposed aṡ Using Darcy's law (Box1.2), the rate of dissipation associated with the flow through the porous medium is given bẏ which is non-negative by definition. Adhering to the assumptions made for the Reynolds flow in Section 2.2, the viscous dissipation of the fluid inside the fracture, which is also non-negative, can be written aṡ where, upon substitution of equation (16), we obtaiṅ Note that the Couette term is associated with fluid flow driven by wall motion. Since the flow profile in the fracture is affected by the wall slip, the Couette term depends on the parameter β. The Poiseuille dissipation is driven by the pressure gradient along the fracture, and is not directly dependent on the wall slip parameter. The dissipation rates due to the skin effect and the Beavers-Joseph-Saffman wall slip in equation (49) are elaborated aṡ where use has been made of the interface conditions (13) and (35). Equation (53a) conveys that the rate of dissipation due to the skin effect depends quadratically on the pressure drop across the skin, with its scaling being dependent on the surface entry resistance parameter, γ . The dissipation due to wall slip in equation (53b) depends quadratically on the slipping speed between the porous medium and the flow inside the fracture, with its scaling depending on the slip parameter, β. It is important to note that since all contributions to equation (49) are non-negative, the rate of dissipation for the considered model is non-negative, i.e., Ḟ ≥ 0.
To complete our energy balance, we furthermore define the rate of mechanical work aṡ Finally, the irreversible work is defined aṡ where the minus sign in front of the first integral is a result of the outward-pointing normal vector. Substitution of equations (44), (49), (54) and (55) into the energy rate relation (43) then yieldṡ which, using the formulation summarized in Box 1, is equal to zero in correspondence with equation (43). See the proof presented below for details.
Proof that the energy rate balance is satisfied. To proof that the total energy rate balance (43) is satisfied, we substitute the following expressions in equation (56): • Multiplication of the momentum balance (Box1.3a) with a displacement rate u and integrating over the domain where, after integration by parts, use has been made of the boundary conditions (Box1.3b), (Box1.3c) and (Box1.3d).
• Multiplication of the porous medium mass balance (Box1.1a) with the pressure p and integrating over the domain • Multiplication of the mass balance in the fracture (Box1.5a) with the fracture flow pressure p c and integrating by parts where use has been made of the boundary condition (Box1.5c) and curvature effects have been neglected.
• Multiplication of the flux through the fracture (Box1.5b) with the fracture pressure gradient ∇ s p c and integrating over the fracture domain c yields After substitution of the above and (13) into (56) we obtain: Finally, using that it follows from (16) that in combination with equation (34), it follows that (see Remark 8) which upon substitution in (61) proofs that ˙ = 0.

Remark 7 ((ε, p)-formulation of the internal energy).
Using the definition of the fluid content ζ in equation (47), the internal energy density (45) can alternatively be expressed in terms of the strain and the pore pressure as: Remark 8 (Interface integral identity). Consider two arbitrary functions, a and b, which are discontinuous at the interface c . The integral of the product of these two functions over the two sides of the interface, i.e., ± where use has been made of the identities in Remark 6.

Numerical model
In this section we discuss the computational solution procedure for the fractured poroelastic medium model summarized in Box 1. We commence with the introduction of the incremental-iterative solution procedure in Section 4.1. In Section 4.2, we discuss the mixed finite element discretization, where the partition of unity method [37,38] is employed to introduce (fracture) discontinuities and (tip) singularities. Finally, in Section 4.3, based on the energy balance discussed in Section 3, convergence criteria for the incremental-iterative solution procedure are proposed.
# Solve the linear coupled problem (Section 4.2)

Incremental-iterative solution procedure
To solve the coupled problem in Box 1 in terms of the porous medium displacement, u : where T is the total simulation time), we consider the fully-coupled incremental-iterative solution outlined in Algorithm 1. In comparison to Bergkamp et al. [4], the velocity field in the fracture is not considered as a primary field variable (i.e., it can be evaluated based on the above-mentioned field variables through equation (16)). This simplifies the satisfaction of the interface coupling conditions (Section 2.3), making it straightforward to solve for the fully coupled system. A staggered solution procedure, as employed in, e.g., Ref. [4], is therefore not considered convenient in the dimensionallyreduced fracture setting of this manuscript.
To discretize the problem, we consider time steps of size t, such that t = ı t, with time index ı = 0, . . . , n t and n t the number of time steps (excluding the initial condition corresponding to ı = 0). Within each time step, we employ Picard and time derivatives are approximated bẏ To linearize the fully coupled (non-linear) system of equations in Box 1, within each Picard iteration the fracture aperture is fixed to the value at the previous iterate as After n t time steps we obtain the discrete solution

Finite element discretization
The finite element method is used to discretize the coupled problem in Box 1. The weak formulation on which this discretization is based is derived in Appendix A. We consider Taylor-Hood elements for both the porous medium problem defined over the domain and for the fracture flow problem defined on the interface c . We employ C 0 -continuous piecewise (bi)quadratic polynomials for the displacement/velocity field, and C 0 -continuous piece-wise linear polynomials for the pressure fields. To represent the discontinuity in the porous medium fields across the fracture, the displacement field and pore pressure field are enriched using the partition of unity method [37,38]. Our implementation of the partition of unity method is based on the phantom-node concept of Ref. [40]. We use a level set description (e.g., Ref. [12]) to determine the geometry of the fracture and the integration procedure of Ref. [41] to evaluate integrals on the cut elements. Following Ref. [42], we also enrich the solution spaces for the porous medium problem with tip enrichments to represent the singular behavior of the fields at the fracture tips. We express the discrete fields (denoted by a superscript h, which refers to the mesh-size parameter) as: follows from the linear system of equations The block matrices in this linear system are defined as and the right-hand-side vector is given by The components of the matrices and vectors in (74) and (75) are defined in Box 2. The derivations are presented in Appendix A.

Energy rate convergence criterion
The energy rate balance (43) introduced in Section 2 can be written aṡ with the vector of energy rate balance components defined aṡ In Section 3 it is demonstrated that in the infinite dimensional setting (no spatial discretization errors), the energy rate balance (77) is equal to zero by definition. In the finite dimensional setting, the fracture flux equality (60) is not satisfied identically, because the surface gradient of the pressure is, in general, not in the considered test space for the fracture flow flux. The energy rate associated with this discretization error is given bẏ In the discretized setting we augment the array (78) with this error term, i.e., such that the sum of this vector evaluates to zero up to numerical precision (e.g., numerical integration errors) when the nonlinear system of equations has been solved up to numerical precision.
To satisfy the global energy rate balance (77) it suffices to augment ˙ with Ė h . It is important to note, however, that although the global energy rate balance is satisfied in the discrete setting, this does not hold for the local balance laws. The employed Galerkin-based discretization using inf-sup stable elements [43] does allow for local fluctuations in the balance laws. The employed elements are known to be stable with well-understood asymptotic mesh-convergence behavior in the context of the Biot problem [44,45]. We opted to also employ an inf-sup stable discretization of the two-field fracture flow model (similar to Ref. [46,47]), based on the observation that this mixed formulation reduces the local fluctuations in the flux field compared to a single-field formulation (which would allow for inter-element discontinuities in the flux field if a C 0 -discretization of the pressure field is employed). It is noted that the rigorous inf-sup stability proofs for the considered mixed elements do not trivially extend to the situations with discontinuities and tip-singularities, such as considered in this work. We do observe all our simulations to yield stable results, provided that sufficiently fine time step and mesh sizes are used. During the Picard iterations, i.e., when the nonlinear system of equations is not yet solved up to numerical precision, the total energy rate balance is violated. This makes the energy rate balance a suitable criterion to test convergence of the Picard iterations. We propose two energy rate balance relations to test convergence: The first criterion checks whether the total energy rate is smaller than a specified tolerance, tol (with unit W ). The second criterion checks whether, for all components, the variation between two Picard steps is smaller than the specified tolerance. This second criterion is necessary, as it may occur that the first criterion is satisfied, because errors in two (or more) of the components in the energy rate array (80) cancel out.

Numerical simulations
This section outlines the simulations performed to analyze the behavior of the numerical model, as derived in Section 4. An additional benchmark simulation is presented in Appendix B. The simulations considered here pertain to the reopening of an initially closed stationary, i.e., non-propagating, fracture in a permeable medium, in response to the injection of fluid. To study the interaction between the fracture flow and the permeable medium, we study the fracture volume rate balance and the energy rate balance of the coupled system for varying fluid entry resistance and Beavers-Joseph-Saffman slip rate parameters. To be able to provide a comprehensive analysis of all relevant volume rate and energy rate terms involved, we employ the two-dimensional geometry depicted in Fig. 3 for all presented test cases. In this setting, the porous domain is For notational brevity, we here and in the remainder omit the superscript h indicating discretized fields. Throughout the test cases, we utilize the model parameters listed in Table 1, from which we can derive various other parameters, as denoted in Remark 9. The considered parameters are selected to be of an order that is representative for applications of the model in the subsurface setting, specifically, for fractures located in relatively permeable sandstones (see, e.g., Ref. [48]). At the center of the fracture, fluid is injected at a volumetric rate I , which we model as a point source, i.e., ς = I δ(x, 0), with δ denoting the Dirac function centered at the origin of the coordinate system x. Note that both the fluid entry resistance parameter γ and the Beavers-Joseph-Saffman slip rate parameter β are not listed, as these are varied throughout the presented test cases. Where applicable, results are reported per unit depth of the two-dimensional domain.
We commence our study by analyzing the behavior of the model under the influence of both an imposed fluid entry resistance and a wall slip in Section 5.1. With the exception of the mesh verification study discussed below, all presented results are based on a 60 × 61 elements mesh. Using the Taylor-Hood elements discussed in Section 4.2, this leads to a system with 11,105 unconstrained degrees of freedom. A time step size of t = 1 s is used, and the tolerance for the Picard iterations is set to 1 W, which is approximately three orders of magnitude lower than the considered injection power. After considering the verification aspects in Section 5.1, we will vary the fluid entry resistance and Beavers-Joseph-Saffman slip rate parameter in Section 5.3 and Section 5.4, respectively. Table 1

Model behavior with fluid entry resistance and wall slip
We consider the coupled domain depicted in Fig. 3. The behavior of the coupled model is controlled by the model parameters listed in Table 1, along with an imposed fluid entry resistance, γ = 1 × 10 10 kg/m 2 s, and a finite fluid slip along the fracture walls, β = 0.01 [-]. At t = 0 s, fluid is injected into the center of the fracture at a constant rate, after which a balance will be found between opening of the fracture and fluid leaking into the formation. Fig. 4 shows the response of the fracture to the injection of fluid, at three time instances. Results are depicted shortly after injection has started, at t = 1 s, at an intermediate time step, t = 10 s, and when the solution is approaching steady state, at t = 100 s. In Fig. 4a, we see the pressure response to injection. The blue line represents the fluid pressure in the fracture, p c , the light gray line represents the average of the pore pressures at the two sides of the fracture {p}, and the dark gray line represents the jump in pore pressure between the two sides of the fracture p (as defined in Section 2).
Due to symmetry, p + = p − , {p} can be interpreted as the pressure at either wall of the fracture, and p = 0 Pa. At t = 0 s, the fracture is completely closed. As fluid is injected into the center of the fracture, the pressure in the fracture starts to increase, as depicted in Fig. 4a. Initially, the inflow of fluid causes the pressure in the center of the fracture to peak, which is a direct result of the injected fluid not yet being able to flow through the initially closed fracture. Due to the compressibility of the porous medium and the aperture-dependent conductivity of the fracture flow, pressure diffusion through the fracture and the pores is a time-dependent process. The delay in fluid flowing into the porous medium is exacerbated by the fluid-entry resistance. This effect remains visible as the system reaches steady state, as evidenced by the pressure jump between the fracture and the porous medium which is still present at t = 100 s. As the pressure in the fracture increases, the fracture opens, as shown in Fig. 4b.  Fig. 4a. This causes a suction of pore fluid into the fracture, corresponding to a dip in pore pressure at the fracture walls [49]. Note that since all pressures are taken relative to the reference pressure of 0 Pa, the negative pressures are not negative in the absolute sense. Eventually, the increase in fracture fluid pressure throughout the fracture causes the fracture to open along its entire length. Since fracture propagation is not considered in our model, the fracture remains closed at the tips. At t = 100 s, the fracture is completely opened by the fluid, where the pressure decay between the center and the tips of the fracture is caused by leak-off to the   Fig. 4. Here, we also observe that initially, the injection of fluid into the center of the fracture yields an increase in pore pressure near the injection point, a dip in pore pressure ahead of the pressure diffusion waves propagating through the fracture, and an instantaneous response in porous medium deformation. As the fracture fluid pressure increases and propagates throughout the fracture, the pore pressure increases and the porous medium deforms further. As also observed in Fig. 4a, there will remain a jump between the fracture fluid pressure (as shown in Fig. 4a) and the pore pressure surrounding the fracture, due to the influence of the imposed fluid entry resistance. Fig. 7 shows the fracture mass balance and the energy rate balance throughout the simulation. The fracture mass balance in Fig. 7a shows all components defined in Remark 3. The mass balance contribution associated with the compressibility of the fluid is not visible in the figure, as it is negligible at the considered pressures (around 1 MPa). It is observed that initially the leak-off contribution is small. This is a consequence of the fluid not yet being transported through the fracture. As the fracture opens, the surface through which fluid can leak off increases, resulting in a gradual increase in the contribution of leak-off to the mass balance. Toward the final time of the simulation, when the fracture has opened completely, the majority of the injected fluid leaks off into the reservoir. Following equation (Box1.1d), the observed final total leak-off of  approximately 0.9 L/s is in agreement with the observed pressure difference between the fracture and the porous medium in Fig. 4a, which is on average equal to 1.1 bar. Fig. 7b displays the energy storage and dissipation contributions discussed in Section 3. It is observed that the net energy rate, ˙ , equates to zero as expected. Since the injection rate is constant throughout the simulation, the initial dip in the irreversible work associated with the injection process is a direct consequence of the initial pressure drop because of the pressure diffusion through the fracture and the reservoir. In line with the observations on the leak-off behavior from the mass balance, the dissipation associated with the Darcy flow inside the porous medium also gradually increases, and becomes the dominant dissipation mechanism toward the end of the simulation. Since a steady leak-off process develops, the dissipation mechanisms associated with the Poiseuille flow and the skin effect also approach a non-zero steady state. The observed reduction in the slip dissipation is a consequence of the increasing aperture of the fracture. Finally, note that since the shear motion of the fracture walls is negligible, the contribution of the Couette dissipation is not observable in the figure.

Verification
In this section, we elaborate on our choice of time step size, verify our choice of mesh size, and analyze the convergence of our Picard iteration technique.
To capture the relatively large gradients in pore pressure that occur near the fracture, a minimal time step size is required [50]. For a backward-Euler time integration scheme, this minimal time step size is given by where the consolidation coefficient c is defined as To demonstrate the influence of the employed mesh size, h = 1 m, on the accuracy of our numerical approximations, we compare our results to results obtained with a coarser mesh, with mesh size 2h, and to results obtained with a finer mesh, with mesh size h/2. Comparing the contour plots of the porous medium solutions p and |u|, we notice no observable differences. Therefore, these plots are not shown here. Comparing the fracture solutions p c , {p} and n , as depicted in Fig. 8, the only variations we observe are in the smoothness of the solutions. We observe these variations, for example, when focusing on the maximum and minima in p c . However, these variations are not considered to be significant to our study. The lack of dependence of the solutions on the mesh size is partly due to the use of tip enrichments, which enhance the accuracy of the approximate solution near the fracture tips. We continue our mesh study by considering the balance laws. The fracture volume rate balance results show no observable differences between mesh sizes and are therefore not shown here. The variations in the energy rate balance are minor, and in line with the inaccuracies expected to be introduced in the discretized setting. To demonstrate these variations, we focus specifically on the non-physical energy rate associated with the fracture momentum balance, Ė . Under mesh convergence, Ė should equal zero for all time instances. Looking at Fig. 9, we can conclude that this is indeed the case. Although our Taylor-Hood elements struggle to capture the initial peak in fracture fluid pressure, Ė decreases over time, and goes to zero under mesh convergence. Important to note here is that the range of values for Ė depicted in Fig. 9 is an order smaller than the range of most of the energy rates plotted in Fig. 7b. Based on the small variations between the results for the various mesh sizes, we can conclude that our mesh size of h = 1 m is indeed adequate. Since computation time is not a limiting factor in our simulations, we chose not to coarsen our mesh.
The convergence behavior of the Picard iterations is illustrated in Fig. 10 for the first 10 time steps. This figure conveys that the total energy rate (81), which should go to zero, is indeed gradually reduced by the Picard iterations. The observed linear convergence rate is in agreement with the theory on fixed point iterations (see, e.g., Ref. [39]). The observed difference in finally achieved energy rate is a consequence of the fact that the second criterion in (81), which monitors the variations between the various terms in (80), is critical in the case considered here. For all simulations presented herein, the tolerance in (81) is set to 1 W, meaning that for the results in Fig. 10, all energy rate contributions in (80) vary less than 1 W in the last Picard iteration. The observed convergence behavior, with the Picard solver converging within a maximum of 6 iterations, is representative for all conducted simulations.

Influence of fluid entry resistance parameter γ
In this section, we consider the influence of the fluid entry resistance parameter γ on the fractured reservoir problem. We do this by comparing the fluid entry resistance considered in Section 5.1, γ = 1 × 10 10 kg/m 2 s, to a lower fluid entry resistance, γ = 1 × 10 8 kg/m 2 s, and a higher fluid entry resistance, γ = 1 × 10 12 kg/m 2 s. A decrease in fluid entry resistance is expected to increase leak-off from the fracture into the reservoir, resulting in a decrease in pressure jump between the fracture and the reservoir, and vice versa. Fig. 11 shows the fracture solution at t = 100 s for the considered range of fluid entry resistances. Looking at Fig. 11a, we observe that for a relatively low γ , the fracture fluid pressure and the pore pressure at the fracture walls overlap, i.e., there is no pressure jump between the fracture and the porous medium. As γ increases, the pressure jump is observed to increase as well, which is to be expected, since fluid is restricted from flowing into the reservoir. Looking at Fig. 11b, we observe that an increase in fluid entry resistance results in a larger fracture aperture n . This can be explained by the increase in fracture fluid pressure, while the stiffness of the reservoir is kept constant. We note that for the high fluid-entry resistance case, i.e., γ = 1 × 10 12 kg/m 2 s, the observed maximum fracture opening of approximately 2.8 mm approaches that of an impermeable fracture in a solid medium (∼ 3.0 mm, computed by a linear elastic finite element simulation using the domain and mechanical boundary conditions of Fig. 3 and a pressure of 31 bar). The observed difference can be attributed to the fact that the considered case is not completely impermeable, as for example seen from the non-zero pore pressure in Fig. 11a. It is also noted that there is a substantial influence of the finite domain size on the computed fracture opening, disallowing a comparison with the analytical infinite domain solution [51].
Based on Fig. 11, we can conclude that γ controls the balance between leak-off to the porous medium (which is related to the pressure jump over the walls of the fracture via the fluid entry resistance condition (Box1.1d)) and opening of the fracture. As γ increases, the fracture flow velocity tangential to the fracture plane, v c (ξ ), decreases, but the fracture aperture n increases. As a result, the fluid flux through the fracture is not fundamentally affected by the fluid entry resistance parameter. Fracture flow profiles are therefore not shown here.
The dependence of the fracture flow mass balance and the energy rate balance on the fluid entry resistance parameter is illustrated in Fig. 12. The above-mentioned influence of the fluid entry resistance on the leak-off behavior is clearly observed in Fig. 12e, which shows the case of a relatively high fluid entry resistance. Compared to the lower fluid entry resistances shown in Figs. 12a and 12c, the high resistance case results in a significant reduction in fluid leak-off. The increased fracture opening as discussed above is a direct consequence of this behavior.
The increase in fracture opening, and the corresponding increase in fracture fluid pressure, are also observable from the energy rate balances in Figs. 12b, 12d and 12f. For the high fluid entry resistance case in Fig. 12f a notably higher injection power is observed. This is a direct result of the higher pressure, as the injection volume rate is kept constant between the considered cases. In line with this observation, with a skin layer, the internal energy rate in the porous medium increases significantly compared to the case where the fluid enters the reservoir without additional resistance at the fracture walls. It is also observed from the energy rate balance that by increasing the fluid entry resistance parameter, the dissipation rate associated with the skin effect increases significantly, and starts to dominate the dissipation associated with the Darcy flow in the porous medium. A final observation from the energy rate balances in Fig. 12 is that by increasing the fluid entry resistance, two distinguishable stages in the injection process become apparent (see Fig. 12f). In the first stage, approximately up to t = 17 s, the fluid transports through the fracture. This redistribution of the fluid through the fracture happens in the absence of significant pressure fluctuations. At approximately t = 17 s, the pressure front inside the fracture hits the fracture tips, limiting the possibility for further redistribution of the fluid inside the fracture. As a result, continued injection after this point translates to an increase in fracture opening and corresponding increases in fracture pressure and injection

Influence of Beavers-Joseph-Saffman slip rate coefficient β
In this section, we consider the influence of the Beavers-Joseph-Saffman slip rate coefficient β on the fractured reservoir problem. To study the influence of β, we compare results obtained for β = 10 −2 [-], with results obtained with a lower and a higher Beavers-Joseph-Saffman slip rate coefficient, β = 10 −4 [-] and β = 1 [-], respectively. As β → ∞, we expect the fluid slip along the fracture walls to decrease, until a no-slip situation is reached. Fig. 13 shows the fracture solution at t = 100 s for the considered range of Beavers-Joseph-Saffman slip rate coefficients. From Fig. 13a we can conclude that for a relatively low slip rate coefficient, β = 10 −4 [-], the fluid pressure is constant throughout the fracture. However, if we increase β, we observe an effect on the pressure dissipation throughout the fracture, caused by the increased friction imposed on the fluid by the fracture walls. The change in pressure dissipation throughout the fracture also has a slight effect on the fracture aperture profile, as observed in Fig. 13b.
The effect of the Beavers-Joseph-Saffman slip rate coefficient becomes most apparent when looking at the flow profiles in the fracture, as depicted in Fig. 14. As expected, the flow profile can be scaled between full-slip and no-slip at the fracture walls, by varying β. As β increases, the slip along the walls decreases. Since the flux through the fracture stays almost constant (apart from a slight influence from the change in pressure dissipation throughout the fracture, which causes variations in leak-off to the porous medium), a decrease in wall slip results in a higher peak velocity in the center of the fracture. It is noted that the observed flux is in good agreement with the observed leak-off and fracture aperture rate.
For example, for the no slip case, β = 1 [-], the observed pressure gradient in Fig. 13a at x c = 10.5 m corresponds to a flow rate of approximately 0.26 L/s (using Remark 4), which is in good agreement with the flow profile in Fig. 14b. Fig. 15 displays the fracture mass balance and energy rate balance for the various slip parameters. It is observed that the slip behavior has a limited influence on the mass balance. For the high slip case, as shown in Fig. 15a, there is less resistance for the fluid to transport through the fracture, and hence the complete fracture can be opened by fluid faster than for the cases with a higher wall friction in Figs. 15c and 15e. This results in a more rapid increase of the leak-off contribution at the start of the injection, as observable in Fig. 15a. This different initial behavior is also observable in the energy rate balance, for example by the lower injection power at the start of the simulation. This is a result of the fact that a reduction in fracture wall friction allows fluid to be transported through the fracture more easily, which consequently reduces the pressure build-up. Another noticeable difference in the energy rate balance is that the Poiseuille dissipation, which is associated with the friction at the fracture walls, reduces as the slip parameter reduces. For the considered intermediate value of the slip parameter in Fig. 15d, energy is dissipated through the slipping process. Further reduction of the slip parameter, as  Fig. 15b, results in the dissipation of the fracture flow to disappear altogether, as fluid can then be transported through the fracture without resistance.

Conclusions
We have developed a model which couples a discrete fracture to a porous medium. A dimensionally-reduced fluid flow model is proposed to couple the flow in the (zero-thickness) fluid domain to the surrounding porous medium through appropriate interface conditions applied on the walls of the fracture. Specifically, we allow for modeling the influence of a skin layer of reduced permeability that may be present on the fracture walls by scaling the fluid entry resistance parameter, and we consider a slip rate parameter (Beavers-Joseph-Saffman slip rate coefficient) to model the slip of the fracture fluid flowing along the fracture walls.
By reducing the dimensionality of the fracture we allow for upscaling to a realistic subsurface setting in which there is a considerable discrepancy in length scales between the fracture aperture and the porous medium domain size. To numerically model the coupling between the discrete fracture and the formation, we apply an extended finite element method (X-FEM), employing the Enhanced Local Pressure (ELP) method to model the jump in pressure over the fracture walls induced by the fluid entry resistance parameter. As part of our numerical model, we identify all energy storage and dissipative mechanisms in the fractured poroelastic medium problem in a thermodynamic framework. We use this thermodynamic framework to establish convergence of the Picard iteration technique which is employed to solve the nonlinear coupled problem. Furthermore, the thermodynamic framework provides detailed insight into the influence of the model parameters on the different dissipative mechanisms.
Both the influence of the fluid entry resistance parameter and the Beavers-Joseph-Saffman slip rate coefficient have been analyzed for a range of values. The reduction of the dimensionality of the fracture is shown not to be limiting to the modeling of skin and fluid slip effects. The numerical model is shown to be robust with respect to the corresponding model parameters, covering both cases of negligible fluid entry resistance and zero slip at the fracture walls, and cases where the influence of both skin and slip effects are considerable.
Both the fracture volume rate balance and the energy rate balance in the system are studied. The fracture volume rate balance is shown to be especially sensitive to the scaling of the fluid entry resistance parameter, where a relatively high parameter value corresponds to significant clogging of the fracture walls, resulting in severely limited leak-off from the fracture into the reservoir. Looking at the energy rate balance, a reduction in leak-off is shown to result in more elastic energy being stored in the poroelastic structure. The Beavers-Joseph-Saffman parameter is shown not to have a significant effect on the fracture volume rate balance. However, looking at the balance of dissipative mechanisms within the system, adjustments to the slip at the fracture walls are shown to cause a shift in balance between the various dissipative mechanisms within the fracture fluid flow. The numerical observations in this study are obtained for a single planar fracture, modeled in two dimensions. In a qualitative sense, we expect that the observations are also representative for situations that are not fundamentally different in terms of physical behavior (e.g., non-planar fractures with moderate curvature). Extension of the model to multiple (interacting) fractures, possibly in three dimensions, would require the incorporation of additional model components (see, e.g., Refs. [52,53]). The extension of the observed results to such cases is therefore not trivial. With respect to the extension of the developed model to the three-dimensional setting, apart from the standard computational effort increase, the main challenge lies in the complexity of the implementation.
The presented parameter study demonstrates that the developed model is capable of mimicking a wide range of parameter values. In practice, the physical model parameters (specifically the Beavers-Joseph-Saffman slip parameter and fluid entry resistance parameter) would either be inferred from experiments, or by upscaling from smaller scale models (e.g., a pore scale model). Especially in situations where the fracture cannot be inspected, this is a challenging task. The availability of a model enables the calibration of model parameters based on observables that are sensitive to these parameters.
As both the fluid entry resistance parameter and the Beavers-Joseph-Saffman slip rate coefficient relate to the (physical) properties of the fracture walls, they are not expected to vary independently. The relation between the two interface effects is a topic of future study, along with the study of non-Newtonian fracture fluids and the study of clogging of the fracture walls via the transport of particles through the flow. Finally, the current model provides a solid basis for the extension to fracture propagation, for which common methods generally rely heavily on the reduced-dimensionality of the fracture. When studying the propagation of fractures, the developed thermodynamic framework is envisioned to aid the analysis of the trade-off between fracture propagation and dissipation of energy to the pore fluid through leak-off to the porous medium.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
The porous medium momentum balance The weak form of the porous medium momentum balance is obtained by multiplying the momentum balance (Box1.3a) with a test function w ∈ U 0 and integrating by parts, which gives Substitution of the boundary condition (Box1.3c) and interface conditions (Box1.3d) then yields where use has been made of the identity in Remark 8. Substitution of the constitutive relations, (Box1.4a), (Box1.4b) and (34) gives where the operator B is defined in (A.3a), and where A (u, w) = 2μ∇ s u : ∇ w + λ(∇ · u)(∇ · w) dV , The fracture flow mass balance The weak form of the momentum balance for the flow inside the fracture is obtained by multiplication of equation (Box1.5a) with a test function r c ∈ H 1 s ( c ) to obtain where Gauss' theorem for manifolds has been used with the assumption that curvature effects can be neglected, and the boundary condition (Box1.5c) has been substituted. Using the operators (A.3e) and (A.7c) gives

Appendix B. Comparison with fully-dimensional model
We compare results obtained for the presented discrete-fracture model to results obtained for the fully-dimensional fracture model discussed in Bergkamp et al. [4]. This comparison is presented for benchmarking purposes and to highlight the differences in applicability between the two models. To compare the models, we employ the test case introduced by Ambartsumyan et al. [55], which concerns the injection of fluid into a fractured porous medium. A schematic representation of the domain is shown in Fig. B.16.
In Bergkamp et al. [4], the test case of Ref. [55] is presented for the fully-dimensional fracture setting. The coupled domain employed for this setting is depicted in Fig. B.16a. The domain has dimensions of approximately 1 m × 1 m.
For details regarding the geometry of the domain, the reader is referred to Bergkamp et al. [4]. In the fully-dimensional model, the fracture flow described by the Stokes equations (in this appendix denoted by the subscript s) is coupled to the surrounding poroelastic formation described by Biot's theory, with Darcy flow. In this coupling, the same interface conditions are considered as in Section 2.3, including the Beavers-Joseph-Saffman slip condition and Showalter's fluid entry resistance parameter.
In the discrete-fracture setting, the fracture is reduced to the midline of the fully-dimensional fracture, with an initial aperture, 0 n , equal to the initial fracture opening in the fully-dimensional setting. The domain employed for this setting is depicted in Fig. B.16b. To model the initial aperture, for the simulations in this appendix, the fracture opening in normal direction, n = u · n c , is modified to n = 0 n + u · n c , for all relations in Box 1.
In the comparison, we will employ the material parameters considered in Bergkamp et al. [4], which are taken from Ambartsumyan et al. [55]. The material parameters are listed in Table B    case, a triangulated mesh is constructed on a reference geometry with an elliptic cavity, which is then warped onto the geometry visualized in B.16. For the dimensionally-reduced model the same concept is applied, but then based on a rectilinear mesh in which the discrete X-FEM crack runs through the elements. Following the insights from the verification studies in Ref. [4] and Section 5.2, the time steps and mesh sizes for both simulations have been selected small enough as to make the simulation results insensitive to these numerical parameters.
The results in Figs. B.17 and B.18 convey that there is a close quantitative correspondence between the fully-dimensional and dimensionally-reduced model for both time instances. Although not directly observable from these figures, also the flow velocities, both in the reservoir and inside the fracture, are found to be in very good agreement. This indicates that although the fracture cavity is relatively large (potentially violating the assumptions made in Section 2.2), the dimensionally-reduced model is capable of closely matching the fully-dimensional model. The most notable difference is the observed offset of the pressure field in the direction normal to the fracture surface. This is a direct consequence of the physical initial cavity in the case of the fully-dimensional model. When measured from the fracture surface, the pressure decay into the reservoir closely matches the dimensionally-reduced result, but the "height" of the initial cavity then results in an offset. A similar effect is observed for the displacement field. Away from the fracture surface, the displacements closely match, since the effect of the minor difference in the surface to which the pressure loading is applied is minimal. It is observed, however, that a slightly higher opening is obtained in the dimensionally-reduced setting, on account of the fact that the "height" of the reservoir is increased compared to the case where the initial cavity is excluded from the domain.
Since the observed differences scale with the "height" of the cavity, theoretically, the dimensionally-reduced results could be attained by choosing a very small initial cavity opening. From the perspective of the finite element discretization for the flow considered in Ref. [4], reducing the cavity opening is not practical, as the opening would bound the size of the elements inside the fracture. This would then lead to a dramatic increase in the number of elements needed inside the cavity. This problem is not present in the dimensionally-reduced model, as in this model the thickness direction is not discretized using finite elements, but is instead accounted for in the derivation of the dimensionally-reduced flow model in Section 2.2.