Benchmarking Model Approaches for Thin Structures in a Porous Matrix

Open and closed faults, that are highly permeable or impermeable relative to the porous matrix in which they are embedded, show a strong influence on the fluid flow regime. The modelling of such thin structures may face severe problems, as high demands on storage requirement and execution time. Meshes require extremely fine resolution within the thin object and in its vicinity. The problems are much less severe, if the thin objects are represented in a lower dimension than the model domain. The resulting combined 2D/1D or 3D/2D models constitute a multi-physics approach. For basic benchmark problems we examine the accuracy of such mixed dimension approaches. We compare with the full dimensional model and with analytical solutions. Specific focus lies on the construction of streamlines, for which the streamfunction approach and the particle tracing technique are applied. Finally advantages and disadvantages of the various numerical approaches are summarized.


INTRODUCTION
There are two oppositional situations, in which thin structures appear in porous media: either they are highly permeable, or they are almost impermeable for fluid flow. In the former case we usually speak of open fractures, fissures or cracks, that allow preferential flow within the thin structure in contrast to the surrounding matrix. In the latter case we speak of dikes or veins, that form an obstacle for flow [1]. This study concerns both of the two cases.
Fracture networks within a porous matrix have increasingly come into the focus of research [2,3,4,5,6]. They have been studied extensively, concerning flow fields and hydraulic properties, experimentally, theoretically and numerically. Here we are dealing with numerical approaches. Various different approaches, using Finite Elements, Finite Volumes or analytical elements can be used nowadays. Of particular interest are approaches, in which the thin object is represented in a lower dimension. Aim of the study is to outline advantages and disadvantages of such approaches, with particular respect to the calculation of streamlines.
As thin objects are very small at least in one spatial dimension, the straightforward approach may lead to numerical problems. The contrast between the microscale of the fracture on one side and the macroscale of the matrix requires extensive mesh refinements. The discretization of the thin structure requires small scale finite elements or volumes, i.e. an extreme mesh refinement within the structure and in its in vicinity. In case of many structures or in a 3D domain, the numerical set-up will easily suffer problems because of high storage demands and long execution times.
Thus, it is a good idea to use a lower dimensional representation of the fractures, i.e. as 2D surfaces in the 3D case, or as 1D polygons in the 2D case; for an illustration see Figures 1. This is referred to as the reduced dimension approach [7], or as embedded discrete fracture model (EDFM) [8]. The use of objects of reduced dimension has been proposed and applied in several publications [9,10]. The resulting mixed-dimensional modeling of domains [11] is a multi-physics approach as it combines fluid regimes of different dimensions. Here we deal with the 2D case and examine a numerical approach by which thin structures, faults, are treated as 1D sub-domains. In the text we refer to it also as 2D/1D model. Computations deliver the hydraulic head distributions, potential and streamfunction. From these variable distributions' streamlines can be calculated, which get special attention in this study. Various modeling approaches concerning streamlines are implemented and compared. These may be briefly classified in two groups: streamfunction contouring and particle tracing.
Holzbecher [12] outlined advantages using the streamfunction recently. The streamfunction approach provides a visual information on the amount of flow, provided the levels are taken equidistantly. Another advantage is that the streamfunction as part of the complex potential is utilized directly, by application of contouring. Algorithms for contouring do not produce additional numerical errors.
In contrast particle tracing (or tracking) methods operate in time-steps, using the velocity field obtained from the hydraulic head or potential. Depending on the timestep size an error is introduced, which is propagated in the following tracing steps. For small time steps tracing algorithms deliver good approximations of the streamlines. This also comes at a cost due to higher demands of execution time. Streamlines constructed by tracing methods depend on the user's choice of starting points, and are thus not good representations of flow tubes, as they are delivered by streamfunction contouring.
In the sequel computations are described that are based on various modelling approaches. Numerical experiments are performed for benchmark cases of a porous matrix with a single fracture. We consider different angles between the direction of the fracture and the regional flow. For a given gradient of hydraulic head we examine the dependence of the flow on the ratio of hydraulic conductivities. Particular focus lies on the head gradient along the permeable fracture and the streamfunction jump across the fracture.

Differential Equations
For the porous matrix the validity of Darcy's Law is assumed, which states that the Darcy velocity vector is proportional to the gradient of the hydraulic head ℎ as shown in equation (1): with hydraulic conductivity as proportionality constant. With a scalar value for isotropic media are considered. Utilizing the continuity equation for an incompressible fluid in steady state, the differential equation hydraulic head h results: For flow in the thin structure the validity equation (2) is also assumed. For low permeability faults the conductivity has a significantly lower value. Vice versa permeable fractures may be porous and show an increased conductivity. However, if they are not porous a proportionality between velocity and head gradient can be assumed. For open fractures the formula can be derived, where s denotes the spatial variable along the fracture. The proportionality constant is given by: with fluid density , acceleration due to gravity g, kinematic viscosity η and fracture aperture [14]. Formula (3) results from the Hagen-Poiseuille profile for laminar flow between two plates, which is an analytical solution of the Navier-Stokes equations [15]. For the flux per unit width of the fracture follows as shown in equation (4): which is often referred to as cubic law. Using experimental and field data Klimczak et al. [16] examined the validity of the cubic law. For the matrix-fault system we take a spatial dependence of , in the matrix and in the fracture: The specific discharge potential in a porous medium is defined by = ⋅ ℎ [17]. φ fulfills the Laplace equation ∇ 2 = 0. In 2D is related to the streamfunction Ψ, which is given in equation (6) by the Cauchy-Riemann equations and obeys the differential equation (see also: Crevilién-Garcia [17]): The specific discharge potential and streamfunction can be gathered in a two-valued complex function Φ = + Ψ, the complex potential. The complex potential is a function from the complex plane with complex variable = + , onto the complex plane.
If the fault is represented as an 1D object, can be related to the specific discharge potential in the matrix as noted in equation (8): The flux in a fracture of aperture d is thus given in equation (9): which corresponds to the differential equation: For sake of generality space and time are normalized in the sequel. Unit length is the halflength of the fracture ℓ, i.e. we use the spatial transformations → /ℓ → /ℓ . Unit velocity is the ambient velocity ∞ , i.e. → / ∞ . The corresponding time transformation is → ⋅ ∞ /ℓ. There are six parameters in the given formulation: the matrix conductivity , the fracture conductivity , the velocity at infinity ∞ , the fracture length 2ℓ and aperture , as well as the angle between the fracture orientation and the direction of the regional flow.

Analytical Solutions
For an impermeable obstacle of length 2ℓ in a flow field with velocity ∞ that is tilted by the angle to the main flow direction, Churchill & Brown [19] presented an analytical solution in form of a complex potential Φ(z): which can also be written as: The first term of the expression represents baseflow, the second the change due to the dike. The solution has a jump in the real part of the complex potential. The velocity field is given by gradient: The solution can be adjusted for the case, where a line object with infinite permeability is located within the 1D flow field. In the solution of Churchill & Brown the real and the imaginary part of the complex potential exchange their roles as potential and streamfunction, respectively. Thus, following Sato [9] one can switch to the conjugate complex of the complex potential and obtains the formula The overbar denotes the conjugate complex. Formula (14) can also be written as The explicit formulation is: A derivation of the formulation (16) is given by Weijermars & van Harmelen [20]. This analytical solution has a jump in the imaginary part of the complex potential, i.e. for the streamfunction, and the real part is continuous.
For a situation, in which the fracture is permeable with conductivity , according to Sato [21] the jump condition at the fracture is:

Benchmark Cases
In order to examine various approaches to modelling of thin structures in a porous environment, two benchmark geometries were constructed. The cases are intentionally simple, consisting of a single straight thin structure within a regional flow field. Figure 2 depicts the flow field of two different constellations: an impermeable thin dike and a highly permeable fault.
a surface plot of the real potential φ (yellow for high, blue for low) 3.
contour lines for the real potential φ (black) 4.
a contour plot of the streamfunction Ψ (white) 5.
The constellation in Figure 2 was chosen as a benchmark case for the numerical simulation. There is an angle = /4 between the orientation of the thin structure and the regional direction of the flow. As second benchmark we chose the situations in which the fracture is aligned with the flow field (α=0). In all benchmarks the model region is a square of 10 units side length. The fracture length is 2. The fracture aperture is = 0.01. A hydraulic gradient of 0.07 was chosen for the tilted fracture case and of 0.2 for the case in which the fracture direction coincides with the regional flow direction. Parametric sweeps were performed for / , the ratio between the hydraulic conductivity of the matrix and the fracture.

Head Modelling Approaches
Three approaches for hydraulic head modelling are examined. The implementation of the full dimensional approach in 2D is based on equation (2) with spatial dependent conductivity according to equation (5), in the sequel abbreviated as 2Dh. The mixed dimension approach is implemented in two ways. It is set up manually by combining a 2D model for the matrix according to equation (2) with a 1D model for the fracture according to equation (10), abbreviated as mDh (mixed dimension for head). Another approach utilizes the build-in fracture mode option of the COMSOL Multiphysics software, which is operating with the pressure p as dependent variable (mDp). For the constant density fluids subject in this study the head-pressure relation is given by ℎ = / g. In this implementation the 2D Darcy mode for the porous matrix is combined with the 1D Darcy-fracture mode for the fracture.

Streamline Modelling Approaches
Various methods are used to simulate streamlines. Two fundamentally different approaches can be distinguished: streamfunction contouring and flow path tracing. Advantages and disadvantages of both approaches have been outlined recently [12]. The streamfunction delivers visualizations with equal stream tubes, which provide information on the flow budget that cannot be obtained otherwise. However, the streamfunction exists only in special situations, as for example for potential flow. In contrast flow path tracing can be performed for every given velocity field, wherever it is derived from. For the streamfunction we compare the fully 2D approach for fracture and matrix with mixed dimensional approaches. In the mixed dimensional setup, the 1D fracture is part of the boundary between 2D regions, in which the streamfunction for the matrix is solved. For the benchmark constellations the construction is sketched in Figure 3. In order to implement the jump condition for the streamfunction at the fracture B1 the model region has to be split in two subdomains, as shown in the figure on the left. The line B2-B1-B3 is separating the two sub-domains, in both of which the streamfunction is computed as solution of equation (7). At the boundaries B2 and B3 the match between solutions in the upper and lower subdomains is required by Dirichlet conditions. At boundary B1 a jump is required according to equation (17), as well using Dirichlet conditions. The calculations of the streamfunction in the two subdomains is thus coupled due to the boundary conditions at the interface B2-B1-B3. There is a link to the potential model that determines the jump at B1. There is no back-coupling from the streamfunction to the potential.
with the space coordinates and and the corresponding given velocity components and . Based on the finite difference discretization of equations (18) and (19).
A flow path can be followed from a starting position ( 0 , 0 ) reaching the position (x, y) within timestep . The method delivers approximations of the final position and thus may suffer from error propagation. The accuracy of the approximation depends on the timestep size . For an introduction into the topic see Holzbecher [22]. Generalizations for various kinds of meshes in 2D and 3D have been developed [23,24]. In this study flow path tracing is performed with high accuracy, using automatic timestep refinement if needed, as the focus of the study is the comparison of methods that are implemented and performed with best technical options.

Numerical Model Set-up
The different mathematical approaches are implemented using COMSOL Multiphysics software [13]. The code allows the coupling of various physics modes in different dimensions and domains and is thus appropriate for the aim of this study. Using COMSOL Multiphysics fracture models have already been presented in various references [25,26,27,28,29].
COMSOL Multiphysics is versatile software for the modeling of coupled systems of partial differential equations by the method of Finite Elements. We use standard quadratic Lagrange elements throughout. The meshes are refined at the fracture. A typical mesh is visualized in Figure 3 on the right. Table 1 gives an overview on the meshes used for the different model approaches. The standard mesh was refined in some model runs in order to check the accuracy of the models, i.e. if there is an influence of mesh size on the results. The different mixed dimension approaches were run on the same meshes, in the table denoted as mixed 2D/1D. In all models the 1D fracture was discretized by 120 equidistant elements. The models that are completely set-up in 2D require meshes that are at least one order of magnitude larger than those used in the mixed dimension approaches. The discrepancy will be even higher for networks in which a multitude of fractures is to be considered. The mathematical solvers were run in fully coupled mode. The exception is the streamfunction calculation within the 2D/1D approach, which was run separately as there is a one-way coupling with the potential model only. The linear equations were solved by the MUMPS solver with the default parameter configuration implemented in COMSOL Multiphysics.

Potential & Hydraulic Head
The hydraulic head distribution in the permeable fracture is depicted in Figures 4. It shows the decrease of the gradient with the increase of the / ratio, as depicted in the legend. Figure 4 depicts results obtained by three methods: the 2D approach (2Dh), the manual construction of the 1D fracture (mDh), and the COMSOL built-in fracture mode (mDp). Solid lines represent results from the 2D model. Markers indicate results of the mixed dimension 2D/1D approaches The hydraulic head gradient is decreasing with the ratio / . This behavior is expected, as due to the increased conductivity in the fracture a lower head gradient is required to induce flow of same magnitude. For increasing conductivity ratio, the convergence towards the analytical solution can be expected. According to equation (13) the analytical solution has a singularity at the fracture ends.
Lines and markers represent the different numerical approaches. Obviously, all markers lie exactly on the solid lines, which in this case shows that all three approaches deliver perfectly coinciding results -for both benchmarks and for all conductivity ratios.

Complex Potential Jumps
As outlined above the permeable fracture model has a jump of the streamfunction at the fracture, the impermeable dike has a jump in the real potential. As the two situations are related to each other, it is sufficient here to discuss only one of them in detail. We choose the permeable fracture case and thus the streamfunction jump is scrutinized here. Figure 5 shows the streamfunction jump as computed for both benchmark constellations and some choices for / . Depending on the ratio of hydraulic conductivities the solid lines represent solutions, obtained by the streamfunction model with a jump condition. Markers indicate the solutions obtained by the 2D streamfunction model.
For low ratios the solutions, obtained by both methods differ significantly, but they converge for high ratios. If the permeability contrast is 1000, only small differences can be observed. The numerical solutions converge towards the analytical solution for the infinitely permeable fracture, which is depicted in the figure for comparison. The analytical expression for the streamfunction jump is calculated using the conjugate complex according to formula (14).

Streamlines
Streamlines for both benchmark cases were computed by various methods: 1.
by flowpath tracing in the 2D model, either using the head or the streamfunction results (2Dh) 4.
by flowpath tracing from the combined 2D/1D model using the head distribution (mDh) 5.
by flowpath tracing from the combined 2D/1D model using the pressure distribution, calculated by the Darcy and Darcy-fracture modes (mDp) The streamfunction contours (2Ds and mDs) are shown as solid lines, in black and magenta. As the same number of streamlines is chosen for both visualizations, the lines overlap and the magenta becomes visible only, where there are differences in the output. The magenta color below the black appears mainly in the vicinity of the fracture, where the lines exhibit higher curvature. There is only one exception of this general conclusion: the streamlines obtained by the approach mDs show differences compared to the other methods. This is especially visible in the flow patterns for the fracture that is aligned with the baseflow. In the upper right subplot biggest deviation between the streamlines can be observed directly above and below the fracture. In the lower right subplot, the differences are more pronounced: streamlines from mDs (magenta) show much less inclination towards the fracture. The deviation of results obtained by this approach is obviously bigger for low ratios of hydraulic conductivities.

SUMMARY AND DISCUSSION
Several numerical approaches were compared to simulate and visualize stationary flow fields in the vicinity of a thin structure within a 2D porous matrix. Aside from the full dimensional 2D model we considered techniques by which the thin structures are treated as entities of lower dimension, 1D in this case. Three numerical approaches were applied for hydraulic head: the fully 2D approach (2Dh), two approaches with a reduced dimension (1D) representation for the fractures, one with manual construction (mDh) and one using the fracture mode option of the software (mDp). The streamfunction was computed in 2D (2Ds) and in a mixed dimension model (mDs).
For two basic benchmark cases we scrutinized the results of the methods concerning the hydraulic head distributions and the streamlines, in dependence of the hydraulic conductivity ratio, varying between 10 and 1000. For streamline visualization we utilized the streamfunction or flow path tracing, alternatively. Numerical investigations were performed for thin structures of higher conductivity (cracks), but with re-interpretation concerning potential and streamfunction the results are equally valid for cases in which the thin structure is less conductive, as outlined in section 2. Concerning the hydraulic head distribution, no significant deviations between the different methods were found (see Figure 4.) Concerning the streamfunction both methods, the full dimensional and the mixed dimensional approaches, showed significant deviations within the chosen parameter range. In the results of the mixed dimension approach the jump at the fracture, a measure for the flow within the fracture, is highly underestimated. For high conductivity ratios, i.e. 1000 and higher, the results of the two approaches match. It was shown that solutions from both approaches converge towards the analytical solution that exists for an infinite conductivity ratio, as illustrated in Figure 5.
Flow path tracing either from the 2D model or the 2D/1D approaches deliver identical results (see Figure 6). The streamfunction approaches, 2D and mixed 2D/1D, produce differing output. After the previous discovery of differences in the streamfunction jump results that outcome was expected. However, for large ratios of hydraulic conductivities the mentioned differences become negligible. Calculations with a ratio of 3000 showed only marginal deviations in the flow patterns and the streamfunction jump for all approaches.

CONCLUSION
The reduced dimensional approach for fractures leading to mixed dimensional multi-physics models is a powerful tool that allows flow simulations with less computational requirements. Both 2D/1D implementations for the hydraulic potential led to identical results, which also coincide with the full-dimensional 2D model and compare well with the analytical solutions. For the streamfunction model the results are less favorable: the output of the 2D/1D approach showed significant deviations from the full 2D model outcome. These deviations were detected concerning the streamfunction jump at the fracture as well as the streamline pattern. For conductivity ratios above 1000 the results between all model approaches converged, approaching the analytical solution for infinite conductivity ratios.
Aside from the so far discussed accuracy the streamfunction approach has implementation issues. For multi-fracture networks the model domain has to be sub-divided into a multitude of sub-domains, following the principle sketched in Figure 3 on the right. For complex fracture networks, as the one shown in Figure 1 on the right, the manual construction of such a subdivision is practically not possible. The partition of the domain into sub-domains, as well as the construction of the corresponding extrusion operators would have to be automated, i.e. constructed by additional programming work. Otherwise the presented approach is feasible for a limited number of fractures only. The results shown in Figure 7 were obtained by a mixed dimension approach for the hydraulic head.
Concerning 3D, it has already been noted elsewhere that the streamfunction approach cannot be easily generalized for modeling [30]. Figure 8 visualizes a flow field in 3D. The geometry was already shown in Figure 1 on the left. In the model the fractures are represented as 2D objects with high conductivity. A head gradient is imposed by the boundary conditions. In the red color indicates plot high head values, blue low head values. Moreover, the figure depicts selected streamlines and an arrow field representing the velocity distribution; both in gray color. The computation and visualization were performed using COMSOL Multiphysics software, illustrating the applicability of the mixed dimension approach for 3D/2D networks of thin structures