Embedded discontinuity finite element modeling of fluid flow in fractured porous media

In many geomaterials, particularly rocks and clays, permeability is greatly enhanced by the presence of fractures. Fracture sets create an overall permeability that is anisotropic, enhanced in the directions of the fractures. In modeling the fractures via a finite element method, for example, meshing around these fractures can become quite difficult and result in computationally intensive systems. In this article, we develop a relatively simple method for including the fractures within the elements. Flow through the bulk medium is assumed to be governed by Darcy’s law, and the flow on the fracture by a generalization of the law. This model is embedded in a finite element framework, with the fractures passing through the elements. In this formulation, elements with fractures are given an enhanced permeability in the direction of the fractures. With these enhancements, the material essentially becomes anisotropically more permeable in the direction of fracture sets.


Introduction
There are many instances of fluid flowing through porous solid media. Groundwater flow, oil and natural gas flow in underground reservoirs, geological fault systems, permeation of concrete, and fluids flowing in biological tissues are examples. In all but the last example, the permeability of the system is often greatly affected by localized structures such as fractures, shear bands, dilation bands, and compaction bands. While opening fractures and dilation bands tend to increase fluid flow along the band [16], compaction bands inhibit flow across the band [15]. Some systems, like faults, exhibit more complex behavior, with bands of high permeability outside the fault core due to fracture, and very low permeability across to cataclasis of the particles in the band [1,13]. The result is both increased permeability in the direction of the fault and impaired permeability across it. In this article, we will focus on permeability enhancement in the direction of flow. Modeling the inhibition across the localized regions is saved for another article.
A number of models have been introduced for fluid flow in porous media. At the continuum scale, Darcy's Law is one of the simplest and most widely used. It is an effective model for a wide variety of soils and other materials. In highly permeable materials, where flow may become nonlinear or even turbulent, it may not be accurate. Darcy-Stokes (e.g. [8], among others) or other models may be appropriate. Similarly, modification have been made for clays and other low permeability materials, particularly where the chemical interaction between clay and water molecules becomes a factor. In this article, we will assume Darcy's Law is valid. For complex geometries, a number of numerical techniques can solve the resulting equations in standard ways. The finite element method is the approach we employ in the article.
There are at least two approaches modeling the effect of localized regions on permeability. If there are dense fractures in the medium, modeling each one may prove to be impractical, and the effect could be incorporated into anisotropy of the permeability matrix. Double-porosity models, proposed by Barneblatt [2], which can differentiate between fluid in more permeable fracture sets and less permeable pores, capture fluid flow better than standard single-porosity models in these situations [4,19]. In this article, we are interested in the case of modeling a few, large fractures or other localized zones explicitly. In a finite element context, these fractures may either be meshed explicitly with the bulk material around them, or an intraelement method may be used. While the first approach is certainly viable, there is extra difficulty in meshing around the fractures. The approach is taken in discrete fracture network (DFN) models, for example [17]. Some of these models, for example, [10], neglect bulk permeability, which is appropriate in some applications. While there is significant research in this area, the process is computationally intensive.
Intraelement methods initially were developed for the problems of localization in mechanics. These approaches include enhanced strain methods, for example, [7,11,14], extended finite elements [3], and several other formulations. They have been used to model propagating fractures, shear bands, and other localized structures. It is worth noting, too, that some of these methods have been extended to the case of poromechanics ( [5,12] and references therein), incorporating the effects of fluid flow coupled with solid mechanics.
The focus of this article, however, is to develop an approach just for the evaluation of flow through porous media. As in the methods above, we will assume that the thickness of the fractures is small in comparison to the overall scale of interest. Hence, it will be acceptable to treat the localized region as a line element in two dimensions or a surface element in three. We will develop enhancements to the permeability matrix due to the presence of fractures, and then implement those enhancements into a finite element framework. This article builds on the formulation presented in [6].
The remainder of the article is organized as follows: Sect. 2 will develop the continuum equations for flow with enhanced permeability in localized regions. Section 3 will discuss how those are implemented into the numerical equations. Section 4 will show some example problems, and finally conclusions will be presented in Sect. 5.

Governing equations for porous flow in localized media
We consider a saturated porous medium with fluid flowing through it. The bulk permeability follows Darcy's Law, which we will write as v bulk ¼ Àk bulk rh; where v bulk is the superficial velocity, k bulk the permeability in the bulk material, and h the total head. As is standard in this formulation, we assume that the flow velocity is small and ignore the velocity head.
Next, we add a fracture or set of fractures or other localized regions to the medium. As mentioned previously, these regions have an implied width that is much less than the width of the overall structure. Hence, we will model these regions as being one dimension lower than the structure in which they are embedded and having an implied enhancement to permeability. We will use the term fractures from here on out, with the understanding that this implies to any localized band with enhanced permeability along its direction.
The enhancement to the permeability will also follow a Darcy-like formulation, witĥ v ¼ ÀkðHÞî ð1Þ where quantities with the^indicate enhancements of the permeability due to the fractures. Hence, we can writeî ¼ ðrh Á lÞl where l is the unit vector in the direction of the fracture. Here, H is the thickness of the fracture surface or deformation band. Clearly, the enhancement of the permeability will vary with crack dilation, though we will not investigate the nature of this relationship in this article. The enhancement to the permeability is only active on the fracture, hence the total flow velocity on the fracture can be written where d S is a Dirac delta function acting on the surface.

Differential equations
The governing differential equation, in this form, is not much different than the continuum case.
where G is a fluid source term, such as a fluid injection site. When we assume the validity of Darcy's Law v ¼ Àkrh ð4Þ the equations can be solved. The difference lies only in the modifications to the terms k.

Finite element implementation
The differential equation is solved approximately using the finite element method. We choose total head h as a primary variable. Other researchers have, equivalently, tracked pore pressure. To accommodate the enhancements, we add nodes to the element where the fracture surface intersects the element edges. By connecting these nodes, the element becomes divided. While many techniques have been developed to handle this division, the simplest for the case is to treat each sub-element separately and perform the integration on each side, as shown in Fig. 1. There are a few ways these divisions can be accomplished. For simplicity, we do not add nodes to the band inside the element, allowing a piecewise linear gradient on the fracture surfaces. While not necessary, this approach simplifies the formulation. Since the gradient along the fracture surface will be treated as linear, using high-order interpolation in the bulk elements will not improve the rate of convergence. Of the possible ways to split the element, we choose to split an initially triangular element. The advantage of starting with triangular elements is that there will always be one side with one of the initial nodes and the other side with two. We choose to treat this second side as a quadrilateral sub-element, as in Fig. 1a , though it could also split into two triangular sub-elements (Fig. 1b). The resulting shape functions (Fig. 2) follow the shape functions for each sub-element type. If the corresponding node is not part of both sub-elements, the shape function is zero on the sub-element it does not correspond to. If it is on both sub-elements, it uses the linear shape function for the triangular sub-element and bilinear for the quadrilateral sub-element. Quadrilateral elements can also be divided, but with different results depending how the element is split. If two nodes are on each side of the fracture, as shown in Fig. 1c, two quadrilateral sub-elements may be used. If one corner is cut off, as in Fig. 1d, the triangular subelement may be used on the corner side, and the side with three nodes may be divided into three separate triangles. These are not the only options, but the case is more complex than for triangles.
Elements at the ends of the fracture may simply be split into two triangular sub-elements, as shown in Fig. 3. These elements do not contain the fracture, but allow for compatible head across the element.

Element internal flow vector and stiffness matrix
Since we have chosen a piecewise linear interpolation of the head along the fracture, the gradient along the fractureî can be written where the total head at nodes A and B at the ends of the fracture in the element are h A and h B , respectively, and x A and x B are the coordinates of those nodes. Hence, the first fraction represents the magnitude of the gradient and the second fraction the direction. This quantity can be modified aŝ where Q e ¼ 0 : and h e is the element head vector.
The element internal flow vector may then be written as The first term arises from the bulk permeability in the medium and is equivalent to the internal flow vector for a medium with only bulk permeability obeying Darcy's law. This can be integrated in the standard manner, using Gaussian quadrature on each sub-element.
The extra term comes from the flow along the fracture. Since the head is linear on this surface, a single integration point in the center is sufficient to calculate this integral. An issue that arises here is that B e is discontinuous across the where o=ol is the derivative in the direction of l. For all the nodes not along the discontinuity, the shape function is zero. For the nodes with endpoints on the discontinuity, the shape functions are linear in the direction of the discontinuity path. Therefore, the value of the derivative is a constant ±1/L s , depending on whether it is increasing or decreasing in the assigned direction of l.
Inserting this into Equation (11), we reach Hence, we recover a symmetric formulation for the problem.
As we have already factored out the element head vector, the element stiffness may be easily determined as The stiffness matrix is symmetric. Again, the first term arises from the bulk permeability, and the second from the enhancement to the permeability along the crack. For the given models, the system of equations is linear, and hence, the stiffness matrix is constant. More complex constitutive models will generally result in a nonlinear system of equations.

Three-dimensional formulation
While this article focuses on two-dimensional implementation, the formulation may be extended to three dimensions. In this case, even for the simplest of elements, there is no guarantee that a linear head gradient will exist on the fracture surface, as this would require that the surface is parameterized by only three nodes. In both brick and tetrahedral elements, it is possible that a plane may cut a quadrilateral surface in the element. Still, the hydraulic gradient may be written where l 1 and l 2 are mutually orthogonal directions that are also orthogonal to the normal to the fracture surface. In the finite element implementation, this may be written Fig. 2 Three of the shape functions for the chosen element. a The triangular sub-element uses a linear interpolation on that sub-element and is zero on the quadratic sub-element for the node that is not adjacent to that sub-element. Similarly, b the quadratic sub-element is bilinear, on that sub-element, and is zero on the triangular sub-element if the node is not adjacent to that sub-element. c For nodes on both sub-elements, the linear shape function is used on the triangular sub-element, and bilinear of the quadratic sub-element, forming a continuous, though not smooth, shape function The internal force vector may then be written where Q i = l i T B e . As in two dimensions, B e may be discontinuous across the surface, but is continuous in the directions along the surface. Hence there is no ambiguity in the definition of Q i . The surface integral may be computed with an appropriate number of integration points.

Example 1
The first example is to verify the model. The problem is a rectangular region, measuring 6 m vertically and five horizontally (Fig. 4a). The bulk permeability is 5 9 10 -5 m/s. There is a head difference of 21 m from top to bottom, and no flow is allowed through the sides. There is a fracture that runs vertically from the middle of the bottom of the sample to the centroid of the sample, withk = 5 9 10 -4 m 2 /s. This is verified against the same geometry, but with an area of finite width with increased permeability (Fig. 4b). Initially, this area has a width of 0.1 m and a permeability of 5 9 10 -3 m/s. Multiplied by the cross-sectional area, this will have the same permeability enhancement as the discontinuity, less the small amount of bulk area taken away. The discontinuity is tested with both structured (Fig. 4a) and unstructured (Fig. 4c) meshes.
The resulting head distributions are shown in Fig. 5. The postprocessor does not actually plot the head values on the fracture, reading each element as a three-node triangle. To verify convergence, we examine two quantities, the head one-half meter to the left of the top of the fracture, and the total flow through the body. These are shown in Tables 1 and 2, respectively. As expected near the discontinuity, some refinement is necessary to get an accurate value of the head. To verify refinement, in the second refinement for the mesh with the finite area of increased permeability, the thickness of that area was halved and the permeability doubled. The three types of meshes converge toward a single result.
The total flow also converges, and more rapidly. As an integrated quantity, it is expected that convergence is more rapid, and coarser meshes give adequate results.
The volumetric flow rate on the crack is shown in Fig. 6 for the coarsest mesh and the finest mesh of each type. Since the head is piecewise linear, the flow rate on the crack is constant on each element. Good agreement can be seen comparing the structured, unstructured, and finite width crack approaches, especially on refinement. The fracture picks up more fluid as it moves down, exhibiting square root or logarithmic type increase versus position.
The error of the flow rate on the crack was studied, compared with a very fine mesh. This is shown in Table 3.
In addition, the convergence rate for both the structured and unstructured meshes was studied by using a very fine mesh as an ''exact'' solution. The L 2 norm of the total head appears to halve as the mesh size is halved, as shown in Table 4.
The overall superficial velocity, being related to gradient of the head, converges somewhat more slowly, but also appears to converge. The error is shown in Table 5.

Example 2
The second example contains a diagonal crack, as shown in Fig. 7a. The crack runs 3.5 m up the right-hand side to 1.5 m from the bottom. The material properties for both the bulk material and the fracture are the same as in the first example.
The total head (Fig. 7b) reveals that much of the potential energy is lost across the fracture, as much of the fluid flows down the fracture. The total flow rate through the sample is 11.62 cm 2 /s.

Example 3
We extend the formulation to include curved fractures. In the third example, the boundaries and boundary conditions are the same, but the crack is a parabolic one extending from the left side to the bottom, as shown in red in Fig. 8. Within each element, the crack is considered to have a linear path as developed. This approximation becomes more accurate as the mesh is refined.
The endpoints are found in each element by matching the equation of the curve to the equations of the element edges. In this way, the curve can be tracked as it enters and exits each element, similar to [7]. The procedure does not handle the case of a curve exiting and reentering the element, but in practice this case has not been an issue, especially for fine meshes where the curves are close to linear in each element.
The resulting head distribution is shown in Fig. 8. The head at a node near the crack, located at (0.8926, 4.375), and the total flow for the problem are shown in Table 6. Convergence is somewhat slower than for straight cracks, but adequate.

Example 4
The final example is an example of fluid flow around a pair of sheet piles, as shown in Fig. 9. The sheet piles are 20 m apart. Inside the piles, the soil has been excavated to 10 m  below the water table, and the resulting water flow is pumped out. The soil has been modeled to 30 m outside the sheet piles, and it was assumed that the head at that distance was equal to the head at the top of the water table outside of the sheet pile. It was found that increasing the extent of the soil modeled beyond this distance did not significantly affect the solution.
A fracture runs diagonally through the ground as shown. The permeability enhancement is 10 -3 m/s, while the bulk permeability is 5 9 10 -5 m/s. The bedrock at the base is considered impermeable.
The resulting head distribution is shown in Fig. 9. The total volumetric flow rate in the sheet piles is 4.986 9 10 -4 m 3 /s. With no fracture, the amount is 4.490 9 10 -4 m 3 /s. Hence, even without connecting to the exit surface, we can see that localized flow can have a significant influence on overall flow in a system.

Conclusion
In summary, we have developed an enhanced finite element for capturing fluid flow through porous media with fracture surfaces. The method allows fractures to be inserted into the model independent of the mesh, making the meshing process simpler. The procedure, by subdividing elements with fractures, is conceptually simpler and easier to implement than an extended or assumed enhanced strain type approach.
The model has been validated by comparing it to a localized zone of finite width. The method compares well to both structured and unstructured meshes. It is worth noting that despite the fact that some sub-elements can have poor aspects ratios, the performance of the method is unaffected.
A mesh modification algorithm has been created which, following the global curve of the mesh, locates the endpoints of the curve in each element. Within each element,   Fig. 6 Enhanced flow on the fracture as a function of position (y-coordinate). The left side shows results for the coarsest mesh of each type, the right for the most refined mesh    the discontinuity has been considered to be straight. This is an approximation, but as the mesh is refined, the fracture length and path follow more closely the true curve. This article focuses on Darcy's Law and Darcy-like models for fractured flow. Using this constitutive model results in a linear, symmetric set of equations that are easily solved. The method can, however, be extended to more complex constitutive models, both in the bulk and along the fractured surface. The resulting equations will generally be nonlinear.
Though beyond the scope of this article, the method can be extended to create a jump in the potential field across   Fig. 9 Geometry and total head calculation for a sheet pile example with a fracture the discontinuity. This extension would mimic the behavior seen on some fractures, where permeability is inhibited across the crack due to very fine particle size near the localization zone. The method may also be coupled to mechanical processes, allowing the cross-sectional area of the localized region to change with deformation, bringing about changes in flow rate. Attention will need to be paid to the stability of coupled formulations. Stability in coupled poromechanics has been dealt in several ways, as discussed in [8,9,18], among others. Such extensions are the focus of ongoing research.