Direct simulation vs subgrid scale modelling of fluid flow in fractured or fracturing porous media

The accuracy and the limits of validity of the discontinuous pressure model, which describes fluid flow inside a fracture using a subgrid scale approach, is assessed by comparing simulation results with those from direct simulation using Stokes flow. While the subgrid scale approach assumes a unidirectional flow, the Stokes model includes both velocity components. This is at the cost of meshing the interior of the fracture, which is here achieved through a spline-based mesh generation scheme. This scheme explicitly couples the spline representing the discontinuity to the fracture mesh and thereby alleviates the (re)meshing requirements for the interior of the fracture. The subgrid model and the direct simulation of Stokes flow approaches are compared by simulating a typical case containing a pressurised fracture, highlighting the advantages of using a subgrid model for the range in which its assumptions are valid, and showing its capabilities to accurately include the influence of the fracture on the porous material even outside this range.

While it is known that the cubic law and derived models provide accurate results for the overall fluid transport inside the fracture, no comparison has been made yet between the discontinuous pressure model and the actual fluid behaviour inside the fracture. Herein, we briefly summarise the discontinuous pressure model and describe a method for simulating Stokes flow within deforming and propagating fractures with the aim of comparing these two approaches. Three cases will be assessed, namely a stationary fracture, a propagating fracture with a realistic opening height, and a case in which this height has been increased beyond the assumptions normally valid for the discontinuous pressure model.

Porous medium
We consider a two-dimensional domain Ω p , which is composed of a porous material and is intersected by a fracture, see Fig. 1. The porous sub-domain is described by the displacements of the solid, u, and by the interstitial Fig. 1 Schematic overview of the porous domain and of the fracture domain fluid pressure p p . The fracture is represented as a onedimensional discontinuity Γ d , which allows for a jump in the displacement field and in the fluid pressure across the fracture.

Governing equations
The deformations of the porous material are assumed to occur fast compared to the changes in interstitial fluid pressure. This allows to neglect inertia terms, resulting in the quasi-static balance of momentum for the solid-fluid mixture: ∇ · σ s − αp p I = 0 (1) with I the identity matrix and σ s the Cauchy stress tensor inside the solid material. Linear elasticity and plane-strain conditions are assumed. The fluid flow inside the porous material is governed by Darcy's law. Together with the mass balance this gives the following expression for changes in the pressure due to the solid deformations and a slowly moving fluid: with α the Biot coefficient, M the Biot modulus modulus, μ the fluid viscosity and k the intrinsic permeability. The standard boundary conditions apply: σ s · n = τ on Γ τ or u = u on Γ u (3) q · n = q on Γ q or p p = p on Γ p (4) with τ the prescribed traction, q the prescribed fluid flux, and u and p the constrained displacements and interstitial fluid pressure.

Discretisation
The porous domain is discretised using T-splines [3,36], cast into a traditional finite element format using Bézier extraction [25,35]. The discontinuity is represented using interface elements [39,40] and is propagated along a C 0 continuity line through mesh-line insertion [10,11]. The interstitial pressure is discretised using the cubic T-splines N pp and the solid displacements is discretised using the quartic T-splines N s : The temporal discretisation of the mass balance has been carried out using an implicit backward Euler scheme, evaluating all variables of the current time step at time t + t and discretising the velocity terms as˙ = ( t+ t − t )/ t. These discretisations allow the momentum balance from Eq. 1 to be cast into its weak form and discretised as: while the mass balance from Eq. 2 becomes: with m T = [1 1 0] and B s = LN s , where L is the displacement to strain mapping matrix. The arrays f d and q d represent the forces and the fluxes which stem from the discontinuity, and will be detailed in the next section.

Fractures
The interior of a fracture is represented by a twodimensional domain Ω f , shown in Fig. 1b. This domain uses the local coordinate system (x d , y d ), has a total length based on the current fracture length L t+ t , and has a height corresponding to the fracture opening height at the end of the time step h = n T Γ d N s u t+ t , using the displacement jump u t+ t and normal to the discontinuity n Γ d . The interior of the fracture is described through the fluid pressure p d , and the fluid velocity components v, w.

Governing equations
The fracture is assumed to not contain any solid material and to have smooth, but porous walls. It is assumed that the fluid flow inside the fracture adapts fast to changes compared to the fluid pressure inside the porous material, allowing inertial terms to be neglected. This has previously been shown to be a valid assumption through simulations in which the inertial terms were included [18]. These assumptions allow the fluid behaviour to be described through the Stokes equations: At the boundaries on the top and bottom walls of the fracture we have a no-slip condition: and, using the tangential vector t d , the normal vector n d and the interface permeability of the walls k i , the mass conservation condition across the fracture walls reads: which enforces the velocity at the fracture walls to match the changes in fracture height and fluid leak-off. This fluid leak-off is also imposed on the porous material by imposing the fracture outflow through the top wall: and similar for the bottom wall. By altering this interface permeability, cases ranging from fluid-blocking fractures to fractures containing a near to continuous pressure between them and the surrounding porous material can be simulated. At the inlet of the fracture, a fully developed parabolic velocity profile is imposed: (15) with Q in the total fracture inflow. Finally, the traction used to couple the fluid pressure inside the fracture to the stresses inside the porous material is given by: with the solid part, τ s , of the traction assumed to follow an exponential traction-separation law for the normal component of these tractions, and being zero for the tangential component. It is noted that these boundary conditions are commonly used in fracture flow models. These boundary conditions were used both for the discontinuous pressure model and for the direct simulation of the flow inside the fracture to enable a direct comparison of the results. However, the formulation for the direct simulation of the fracture is such that other boundary conditions can easily be substituted, for instance the Beavers-Joseph-Saffman condition [4,33], which is commonly used for Stokes flow over porous objects.

Discontinuous pressure model
The discontinuous pressure model [12,13,30] assumes the pressure gradient inside the fracture to only depend on x d , while the pressure variations in y d are negligible. It is furthermore assumed that the tangential velocity v is much larger than the normal velocity w, and that changes in v are dominated by y d , such that ∂v/∂y d >> ∂v/∂x d . These assumption originates from the large difference between the fracture length and opening height, and reduce Eq. 9 to: from which the velocity profile inside the fracture is obtained as: Combining this velocity profile with the mass conservation inside the fracture, Eq. 11, integrating over the fracture height, and substituting the normal velocity at the walls with the boundary condition from Eq. 13 results in: which is discretised using Eqs. 5, 6, and p d = N d p d as: with the interface permeability term using a lumped integration scheme to prevent non-physical oscillations in the fracture outflow velocity [19,40] and the other terms using a standard Gauss integration scheme. If this lumped integration scheme was not employed, these spatial oscillations would dominate the solution near the fracture tips for courser meshes [19]. Equations 7, 8, and 20 fully describe the behaviour of the fractured porous material, with the coupling terms between the fracture and porous material given by: using the rotation matrix R to convert between the fracturelocal and global coordinate systems [15], and the integral over Γ ± d indicates that the integration is carried out for the top as well as the bottom walls of the fracture. These equations are combined in a single monolithic scheme, solving simultaneously for the displacements u t+ t , interstitial fluid pressure p t+ t p , and discontinuity pressure p t+ t d .

Direct simulation
For the direct simulation of the fluid flow inside the fracture the interior of the fracture needs to be discretised, and this discretisation needs to be adapted to changes in opening height and fracture length. To achieve this we start of with one-dimensional NURBS which define the discontinuity in the porous domain, and create an additional set of one-dimensional NURBS to represent the fracture height, as shown in Fig. 2a. These two NURBS are then used to create a rectangular mesh in the parametric domain (Fig. 2b), which in turn is mapped onto the physical domain through x d = ξ and y d = η · h/2 for the horizontal fractures considered in this paper (Fig. 2c). Since this discretisation explicitly uses the NURBS defining the discontinuity, the fracture mesh is easily updated once the fracture propagates. Furthermore, the use of a parametric domain as in-between allows the integration to be done over η, thereby incorporating changes in the fracture opening height without the need to regenerate the mesh.
To aid the implementation of the boundary conditions between the porous and fracture domains, the velocity components are split in an interior and boundary part. This, combined with the mesh generation from two separate onedimensional NURBS, discretises the interior of the fracture as: Fig. 2 Overview of the steps performed for the discretisation of the fracture with N 3 and N 4i the sets of cubic and quartic splines used for the height discretisation of the fracture, and N 4b the splines used for the height discretisation that are nonzero at the top and bottom of the fracture. N 1 s indicates that only the first spline is used for the inlet velocity discretisation, and the other splines N 2+ s for the interior discretisation. The set of splines N 2+ s directly implements the no flow boundary condition at the fracture tip by using only the discontinuous splines, whereas N pp and N s use all splines along the discontinuity, thereby allowing for non-zero vertical velocities and pressures at the fracture tip.
Using the discretisations from Eqs. 23-25 the momentum balances (9 and 10) and mass balance (11) are cast into a weak form and discretised as: with the fracture height h = n T d N s u t+ t , and the spatial derivatives in the physical space, required for the mapping x d = ξ , y d = ηh/2, given as: The boundary conditions from Eqs. 12-13 are enforced through their weak forms, and are discretised as: using the normal vector n = [n 1 n 2 ] and the tangential vector t = [t 1 t 2 ], both taking into account large deformations. Unlike the discontinuous model, using a lumped integration scheme for the interface permeability terms is no longer possible due to the non-square matrices. As a result, possible oscillations in the fracture outflow can occur. Finally, the coupling terms used in Eqs. 7 and 8 are given by: The porous domain-Eqs. 7 and 8-and the fracture domain-Eqs. 26-28 and Eqs. 31-32-are solved in an iterative manner until both have achieved a converged solution at t + t, thereby using a fully implicit time discretisation scheme. First, an iteration of the porous domain with a Newton-Raphson scheme is carried out. To approximate the influence of the fracture flow, an undrained assumption is used in which the pressure inside the fracture is altered due to changes in the fracture height [17,28]. Using Eq. 19 the pressure changes inside the fracture can be approximated as: This allows to approximate the derivatives of the coupling terms in Eqs. 33-34 with regard to the interstitial fluid pressure and displacements as: using the matrix N ds defined as u = N ds u, and a stabilising factor 0 < γ < 1. This factor is applied to alter the tangential terms related to the fracture, and prevent oscillations from the iterative scheme. A factor γ = 0.5 worked well for all cases described in this paper, and prevented unstable oscillations that occurred for higher values of γ . After a Newton-Raphson iteration has been carried out using the above tangential matrices, the fluid velocity and pressure inside the fracture are resolved using the newly obtain interstitial pressure p t+ t p and the displacements u t+ t from the most recent Newton-Raphson iteration. Once these velocities and pressure have been determined, the error of the porous domain is checked based on the newly obtained fracture pressure (since the fluid flow within the fracture is linear and the opening height has not changed, the fracture domain is exactly resolved at this point). This error is first used in a linear line-search scheme to improve the convergence. If the error exceeds the convergence criterion afterwards, another iteration of the porous domain is carried out, followed by again resolving the fracture flow. This scheme is summarised through Algorithm 1.
Once convergence has been reached in the porous domain, the stresses ahead of the discontinuity are checked, and if they exceed the fracture criterion σ yy > f t the discontinuity is propagated for a single element. Due to the dependence of the fracture mesh on this discontinuity, this automatically extends the fracture mesh to include the new discontinuity length. The fracture propagation is assumed to occur during the time step, and therefore more iterations are performed to obtain a converged solution at t + t using the new discontinuity length. Finally, once a converged solution is achieved and no fracture propagation occurs, the complete system is considered to have achieved a correct solution and the next time step is initiated.

Results
To compare the discontinuous pressure model with the Stokes flow model, a typical case is simulated, shown in Fig. 3. The domain is 4 × 10 m, with a horizontal fracture originating at the left boundary. In Sections 4.1 and 4.3 the fracture has a length L f rac = 2 m and is not allowed to propagate. In Section 4.2 the fracture has an initial length L f rac = 0.5 m and is allowed to propagate. A fluid inflow Q in = 10 −6 m 2 /s is imposed at the fracture inlet. The porous material is characterised by a Young's modulus E = 20 GPa, a Poisson's ratio ν = 0.2, porosity n f = 0.2, intrinsic permeability k = 10 −16 m 2 , Biot coefficient α = 1.0, and a bulk modulus K s = 10 GPa. The fluid has a viscosity μ = 10 −3 Pa s and a bulk modulus K f = 1 GPa. The discontinuity is characterised by the interface permeability k i = 10 −10 m/Pa s and an exponential traction-separation relation with tensile strength f t = 1 MPa and a fracture energy G c = 1 kN/m.
The temporal discretisation has been carried out using an implicit Euler scheme, using timestep size t = 1 s. The spatial discretisation has 10 × 5 elements away from the discontinuity, with additional refinement layers inserted near the centre. Each of these refinement layers uses half the element length of the previous layer and consists of 4 vertical elements, as shown in Fig. 4. The fracture mesh for the Stokes flow used 10 elements to discretise the fracture height.

Non-propagating fracture
The simulations using the discontinuous pressure model exhibit a quadratic convergence rate, as shown in Fig. 5b. Due to the approximations made for the tangential matrix, and the iterative scheme alternating between the porous and fracture domains, only a linear convergence rate was attained for the simulations for Stokes flow. This difference in convergence rate, combined with the subgrid model only requiring a single domain to be resolved, resulted in much faster simulations using the discontinuous pressure model.
The pressure inside the porous material and the discontinuity are shown in Fig. 4. The results from the Stokes flow simulations show a near to constant pressure over the fracture height, and this pressure corresponds well to the value computed in the discontinuous pressure model. Hence, the model with the Stokes flow and discontinuous pressure model obtain the same interstitial fluid pressure surrounding the fracture, and the displacements of the porous material result in a similar fracture opening height of approximately h = 0.1 mm at t = 200 s.  This pressure in the discontinuity is shown in Fig. 6 for simulations using additional mesh refinement layers. While both simulations attain the same pressure for finer meshes, the discontinuous pressure model yields slightly more accurate results for coarser meshes. This is also seen for the total fluid transport inside the fracture, Fig. 7, where the Stokes flow model shows larger oscillations compared to the discontinuous pressure model for coarser meshes near the fracture tip.
The outflow from the fracture is shown in Fig. 8. No oscillations are observed in Fig. 8b due to the lumped integration scheme applied to the interface permeability terms in the discontinuous pressure model. Since the lumped integration scheme could not be applied to the Stokes flow model, strong oscillations in the fracture outflow can be observed for this model except for higher levels of refinement. However, these oscillations average out and do not influence the pressure inside the discontinuity or the total fluid flux inside the fracture.
Finally, the horizontal velocity inside the fracture obtained through the Stokes flow simulations and the velocity obtained by post-processing the discontinuous pressure model results given in Fig. 9 show no differences between the two models. It can therefore be concluded that the Stokes flow and discontinuous pressure model yield the same results, while the Stokes flow model is slower to converge and requires a finer mesh.

Propagating fracture
Fracture propagation is shown in Fig. 10. Similar to the non-propagating fracture, the discontinuous pressure model yields slightly more accurate results for coarser meshes, but both models converge towards the same result upon mesh refinement. However, for both models the accuracy is not governed by the fracture flow model, but by the element-wise fracture propagation associated with the use of interface elements.
The fluid velocity inside the fracture after 20 min is shown in Fig. 11. For the initial 0.5 m of the fracture no cohesive zone model was used, whereas after this initial length the traction-separation law is present, since this

Opening height
Finally we consider a case in which the fracture opening height is artificially increased by imposing an additional offset [18], so that the opening height h = h 0 (2 − x d ) +n T d N s u . By using this initial opening height fractures which are outside of the usual height range can be simulated, and the limits in which the discontinuous pressure model becomes invalid can be investigated. It is noted, however, that the cases presented in this section are not realistic and merely serve to illustrate the limits of the discontinuous pressure model.
The velocity inside the fracture is shown in Fig. 12 using h 0 = 2 m. While the vertical velocity component is clearly present in the Stokes flow simulation, this component is not included in the discontinuous pressure model as it assumes a unidirectional flow in the interior of the fracture. This is also seen in the horizontal velocity profiles in Fig. 13. While the results for h 0 = 2 cm and h 0 = 20 cm match, the combination of no-slip and interface permeability boundary conditions allows for a horizontal flow component due to the steep fracture walls. However, even though the velocity profile and flow direction inside the fracture is different, the total fluid transport inside the fracture still agrees between both models.
Another effect of the increased fracture opening is the fluid flowing much easier inside the fracture compared to the porous material. The pressure drops shown in Fig. 6 Fig. 13 Velocity profiles inside the fracture at t = 10 s for the discontinuous pressure model (red dashed line, obtained through post-processing) and Stokes flow (blue solid line) at x = 0.5 m are the result of fracture with a 0.1 mm opening. The increased opening heights provided in this section allow for a lower pressure gradient inside the fracture, resulting in an almost constant pressure throughout the complete fracture. As a result, the fluid flow towards the porous material is no longer limited by the transport within the fracture, but solely governedby the interface permeability and the transport inside the porous material. Therefore, even though the fluid velocity looks significantly different between the two models, the effect of the fracture on the surrounding porous material is the same.

Conclusions
Two models have been discussed for the simulation of fluid flow inside pressurised and propagating fractures. The discontinuous pressure model assumes a negligible influence of the fluid velocity normal to the fracture walls, allowing the interior of a two-dimensional fracture to be described as a one-dimensional line. In contrast, the described Stokes flow model simulates both velocity components inside the fracture, but requires the interior of the fracture to be discretised.
A comparison between the two models shows that for a typical fracture case the same results are obtained. Due to the low fracture opening height relative to its length the tangential fluid velocity is much higher compared to the velocity normal to the fracture, justifying the assumptions made for the discontinuous pressure model. Furthermore, a comparison between the convergence rate using the discontinuous pressure model with the Stokes flow model shows the markedly better convergence of the discontinuous pressure model. Finally, the discontinuous pressure model allows for a lumped integration scheme, suppressing fracture outflow oscillations. This is not possible for the Stokes flow model. Finally, simulations have been shown using an artificially large opening height, i.e. outside of the range in which cubic-law based model is justified. The velocity inside the fracture now clearly differs between both models, with the Stokes flow showing a two-dimensional flow pattern, whereas the discontinuous pressure model only exhibiting fluid flow in the tangential direction. However, due to the high opening height both models result in a near-constant pressure within the fracture, yielding the same fracture outflow. Therefore it can be concluded that even though the discontinuous pressure model provides incorrect results for the interior of the fracture, it is able to accurately describe the overall influence of the fracture on the surrounding porous material even for fractures which are such that the underlying assumptions have been violated.