Next Article in Journal / Special Issue
Full-Size Experimental Assessment of the Aerodynamic Sealing of Low Velocity Air Curtains
Previous Article in Journal
Numerical and Analytical Studies of Soret-Driven Convection Flow Inside an Annular Horizontal Porous Cavity
Previous Article in Special Issue
Evaluation of Vortex Generators in the Heat Transfer Improvement of Airflow through an In-Line Heated Tube Arrangement
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Breakthrough Investigation of Advective and Diffusive Transport in a Porous Matrix with a Crack

by
Ekkehard Holzbecher
Department of Applied Geosciences, German University of Technology in Oman (GUtech), Muscat 130, Oman
Fluids 2021, 6(10), 358; https://doi.org/10.3390/fluids6100358
Submission received: 12 August 2021 / Revised: 17 September 2021 / Accepted: 22 September 2021 / Published: 11 October 2021
(This article belongs to the Collection Challenges and Advances in Heat and Mass Transfer)

Abstract

:
Fluid flow and transport processes in fractured porous media are of particular interest for geologists and in the material sciences. Here a systematic investigation is presented, dealing with a generic geometric set-up of a porous matrix with a crack. In such a combined porous medium/free fluid system flow patterns have been examined frequently, while the resulting transport patterns have attracted less attention. Using numerical modeling with finite elements the problem is approached using a dimensionless formulation. With a reduced number of dimensionless parameter combinations (Darcy-, Peclet- and Reynolds-numbers) solution dependencies are examined in parametric sweeps. Breakthrough curves are fitted in comparison to those of 1D model approaches, yielding effective diffusivities and velocities. The computations reveal highest sensitivity concerning the angle between crack axis and flow direction, followed by the Peclet number and the crack axes ratio. As a dimensionless representation is used the results are scale independent. Thus, they deliver estimations concerning effective heat and solute transport parameters that can be relevant in all application fields.

1. Introduction

Fractures, cracks, vugs, and inclusions in a porous medium provide a focus for researchers in many application fields. They are of particular interest in the geo-sciences, and also receive increasing attention in the material sciences. In geology, such systems are most frequently met in the deep sub-surface where the material is consolidated or has been transformed by metamorphic processes. Fractures are of interest from different perspectives, one of which is fluid flow [1,2]. In general, studies deal with the problem as to how inhomogeneities affect flow [3]. Fractures can be conduits or barriers. We adopt the term crack for a structural discontinuity that is a conduit open to fluid flow, in contrast to an anti-crack that acts as a barrier [4]. Broberg [5] defines a crack as a material separation by opening or sliding, with the separation distance substantially smaller than the separation extent—the crack length. In material sciences, flow through fractured porous media is of special concern for the design of technical devices as fuel cells [6,7,8].
Several studies were interested in flow patterns that are induced by small-scale geometric features within the porous matrix. Some deal with a single small-scale entity, which also is the focus here, while others examine more elaborated configurations of different types, structured and unstructured, random or non-random networks etc. Berre et al. [9] gave a review of such studies of flow in fractured porous media. It is a classical task to derive effective parameters that are relevant on the larger scale depending on the parameters of the smaller (micro)-scale [10]—that also is an aim of this study.
While most of the studies deal with flow patterns, fewer are concerned with the effect on the transport processes. In his textbook, Sahimi [11] covered both flow and transport. Here the interest lies on the transport phenomena that are induced by the flow within matrix-crack systems. Studies published with the same aim cover different geometric constellations and mathematical approaches. Bodin et al. [12,13] examined details of transport in a single crack of constant thickness. For a similar set-up, Liu et. al. [14] derived an analytical solution for what they call velocity dispersion, i.e., dispersive effects caused by the velocity distributions. The topic is more recently been re-examined using domain-random walk numerical method [15]. Natarajan et al. [16] explored a sinusodial fracture geometry. More complex constellations with cracks of different geometries were also scrutinized [17].
Altogether the studies cover flow, heat and solute transport, reaction, and deformation. For example, Bagalkot et al. [18] specialized in colloid transport, including gravity effects. Zhu et al. [19] are interested in the effects of reaction rates between the rock matrix and the crack. Lonergan et al. [20] included mineralization in their study. More complex situations covering reactive transport also attracted interest recently [21]. While transport studies in general deal with advection, diffusion, dispersion, sorption, and reaction we restrict our systematic approach to the physical processes, i.e., to advection and dispersion on a larger scale, what can be termed velocity dispersion [14].
Concerning modeling there are various alternatives to deal with fractured media [22,23]. A multitude of random small-scale entities might be described by a continuum approach. Fracture networks are often treated using a multi-dimensional approach: 1D fractures in the 2D space, or 2D fractures in the 3D space [24,25]. This is justified, if higher-dimensional effects can be excluded. The aim here is to examine situations where such a simplification is not possible, so that matrix and crack have to be discretized in the same dimension. Here 2D transport is modeled with clear 2D characteristics (see the illustrating example in Section 5.2).

2. Crack-Matrix Set-Up

We examine scenarios in 2D of a single crack of elliptic shape in a porous medium (matrix). The constellation is sketched in Figure 1, which shows the entire model region in grey and the crack in blue. The crack is assumed to have the shape of an ellipse, a concept that can be found in many studies on the mechanics of fracture propagation and closure [26,27], but that has been adopted in the hydro-sciences [23,27]. The main axis of the ellipse forms an angle α with the flow direction. There are two opposite sides without normal fluxes, neither for fluid nor for solute or heat (in the sketch the horizontal sides). The other sides of the region are open, and the conditions are chosen to induce inflow on one side and outflow on the opposite one. The boundary conditions are chosen in a way that they would induce a 1D constant flow field in a set-up without the crack.
Our aim is a systematic examination of the described system to cover all solution characteristics and their sensitivities within the relevant parameter range. The two main geometric parameters that characterize the system sketched in Figure 1 are the tilt angle α between the main axis of the crack, and the flow direction and the crack aperture d. The numerical simulations are performed to find out how the transport through the system changes dependent upon these parameters.
To cover various situations concerning the physical parameters, i.e., the properties of fluid and porous medium, we examine the system in a dimensionless analytical description [28]. This has the advantage that results can be generalized [29]. In particular, they are scale independent, or scale invariant [30]. Another relevant advantage is that the number of parameters usually is much smaller than the original set [30]. In Section 3, we describe a formulation in which only three dimensionless parameter combinations (Re, Da and Pe) are used for the mathematical description of the equations.
The results are evaluated on the basis of the breakthrough curves, which describe the mean concentration across the outlet boundary as a function of time. The breakthrough curves of the 2D models are used in a fitting procedure with a 1D model approach. In this way, effective transport parameters are determined, allowing a check how far the input value of the 2D model differs from a 1D approach with effective parameters. The final results reveal that aside of the geometrical parameters only one of three dimensionless parameter combinations Re, Da and Pe is relevant for the characterization of the breakthrough curves.

3. Differential Equations

Let us first note the differential equations on which hydrodynamic flow and transport modeling is based. We rely on the equations for Darcy flow in the matrix and for Stokes flow in the open fracture, as well as on the advection-diffusion equation for transport. For the flow we assume stationary conditions, whereas transport phenomena are transient.
The equations for steady-state constant density flow in a porous medium are noted as [28]:
v + k μ p = 0 v = 0
with permeability k and fluid viscosity μ. The first equation is a formulation of Darcy’s law for the velocity vector v, with p denoting the hydraulically relevant pressure. The second equation is a formulation of the principle of mass conservation, under the assumption of constant density. The given formulation holds for stationary flow in porous media in 1D, 2D or 3D.
Inside the crack we assume that the flow is free, i.e., not hindered by any obstacles. Thus, the Navier–Stokes equations for laminar flow are assumed to be valid [25,28].
ρ v v + p μ v + v T = 0 v = 0
where ρ denotes density. denotes the nabla operator, and the subscript T the transpose of a matrix. Equation (2) hold for steady laminar flow, which can be assumed to be valid because flow in porous media systems is usually very slow. Often the first term that accounts for inertia effects is neglected and the equation reduces to the description of Stokes flow [31].
To obtain scale-independent results we construct the models in terms of normalized variables and parameters. Normalized to a reference length L and a reference velocity v 0 , the first equation in Equation (1) can be re-written as:
v 0 v + k L μ p = 0
with normalized pressure p / ρ v 0 2 one can reformulate Equation (3) using dimensionless parameters Da and Re:
v = Da Re p
with Darcy number Da = k / L 2 and Reynolds number Re = v 0 L ρ / μ .
Using the same variable transformations for the first of the Navier–Stokes Equation (2) yields:
v 0 2 L v v + 1 ρ L p v 0 L 2 μ ρ v + v T = 0
with normalized pressure one thus obtains [32]:
v v + p 1 Re v + v T = 0 v = 0
with the Reynolds number as only parameter.
Mass and heat transport can be described by the advection-diffusion equation that can be written as [33]:
θ t = D θ κ v θ
The variable θ represents concentration in case of solute transport and temperature in case of heat transport. D is the diffusivity. The values for the diffusivity D are of a very different scale for solute and thermal transport. For thermal transport the values are in the order of 10−6 m2/s, in mass transport in the order of 10−9 m2/s. The factor κ in the advection term is relevant only for heat transport, where it denotes the ratio of specific heat capacities κ = ( ρ C ) f / ( ρ C ) f s . The ( ρ C ) terms denote the heat capacities where the subscripts refer to the fluid (f) and the fluid-solid (fs) phases. For solute mass transport Equation (7) holds with κ = 1.
For the general treatment of the problem, we normalize the transport equation. With the variable transformations given above the transport equation is written as [28]:
θ t = 1 Pe 2 θ v θ
with Peclet number Pe defined by:
Pe = v 0 L D s o l     for   solute   transport κ v 0 L D t h e r m     for   heat   transport
where the subscripts on D denote the solute and thermal cases. Pe measures the relation between advection and diffusion. For Re > 1 one speaks of advection dominated transport, while for Pe < 1 diffusive processes dominate. Please note that with the normalization of length and velocity as described, the time is implicitly transformed by: t v 0 / L t .

4. Numerical Modeling

Based on the dimensionless formulation, as given by Equations (4), (6) and (8) a numerical model for flow and transport was set up in order to examine the breakthrough behavior. Extensive parametric sweeps were run to study the breakthrough curves. As result of the normalization three dimensionless parameter combinations could possibly be of relevance: Re, Da and Pe.
The model region is sketched in Figure 1. The crack is of elliptic shape with a long axis of length L and a ratio of half-axes of d (smaller axis/larger axis). The region has a height (or width) H and a length that is selected as double of the crack length. Several ratios of geometric entities define the shown set-up. Most relevant are the characteristics of the crack, i.e., the half-axis ratio d and the angle α by which the crack is tilted to the main flow direction. For α = 0 the crack is aligned with the flow, for α = 90° it is perpendicular to the flow direction.
The boundary conditions of the model are gathered in Equation (10). x denotes the coordinate axis in flow direction, while y denotes the perpendicular direction. Because of Darcy’s Law (1) the requirements for the normal components of the velocity v n introduce a Neumann condition for pressure p.
v n ( x = 1 , y , t ) = v 0 ,   p ( x = 1 , y , t ) = 0 v n ( x , y = 0.5 , t ) = 0 ,   v n ( x , y = 0.5 , t ) = 0 θ ( x = 1 , y , t ) = 0 ,   θ / n ( x = 1 , y , t ) = 0 θ / n ( x , y = 0.5 , t ) = 0 ,   θ / n ( x , y = 0.5 , t ) = 0
On the inflow side a Neumann condition is prescribed for flow v n = v 0 , while on the outflow side there is a Dirichlet condition with p = 0. For transport there is a Dirichlet condition θ = 0 at the inflow boundary and the common Neumann condition θ / n = n at the outflow boundary. All other sides have no-flow boundary conditions for both flow and transport. The initial condition for the transient transport model is θ ( x , y , t = 0 ) = 1 .
For the normalization we chose the crack length L and the inflow velocity v 0 . Thus, the model set-up is different to a Darcy experiment, where a pressure gradient is prescribed. The advantage of this formulation is that it is possible to distinguish between effects on advection and dispersion. Details are outlined below.
For the entire set-up we end up with five dimensionless parameters whose relevance is to be examined: Da, Re, Pe, d and α. These parameters are varied systematically in parametric sweeps in order to obtain a general view of the system’s breakthrough behaviour. Reference values in physical dimensions and the variations of the non-dimensionless combinations are listed in Table 1.
The model was set-up using COMSOL Multiphysics 5.6. This software is very versatile for Finite Element numerical modeling of coupled processes in physics, chemistry, biology, geology, medicine etc. The software allows the use of various physical modes (formulations) for flow and transport and free flow and porous media. Most modes are implemented for use of variables and parameters with physical units. For this study, these formulations were transformed for the dimensionless form in the above given Equations (4), (6) and (8).
COMSOL Multiphysics has been applied for porous systems including cracks in previous studies. Using different numerical approaches Holzbecher et al. [34] examined flow patterns in porous media with thin structures. Using a similar 2D setting with a single crack Romano-Perez & Diaz-Viera [35] investigated flow results depending on different representations of the crack (explicit, embedded, domain decomposition), but not transport phenomena as provide the focus here. Perko et al. [36] applied COMSOL Multiphysics in a verification and validation study on a sample with a given crack network. Their model covered flow and transport, similar to our approach. They were able to visualize and validate the complex transport behavior by comparison with experimental measurements. ‘A model validation exercise based on a 2D laboratory tracer test gave convincing results in terms of reproducing plume spreading (qualitative comparison) and breakthrough curves (quantitative comparison)’. The authors used the software in another study concerning the transport of radionuclides in cracked cement [37].
In the Finite Element formulation, the flow variables, pressure and the velocity components were represented by linear elements (P1 + P1), the transport variable by 2-order elements. Meshing was always made with a considerable refinement near and within the crack. Figure 2 shows an example mesh. A typical mesh consists of 3308 nodes, 6520 elements and 9750 degrees of freedom (DOF) for flow and 12,825 DOF for the transport discretization. Checking on the P2 + P1 alternative for the flow discretization led to 2899 DOFs. With an average element quality of 0.88 the meshing is good.
As the study deals with advection-dominated flow with high Peclet-numbers it was necessary to use several stabilization schemes. We only used consistent stabilization, as the streamline and crosswind diffusion. COMSOL Multiphysics proved to be quite effective dealing with problems that may arise during modeling of advection dominated transport [38]. The temporal simulation is performed using time-stepping. The time steps are updated during the simulation, following the Courant-Friedrichs-Lewy (CFL) condition.

5. Results

5.1. Flow Pattterns

Flow in the combined matrix-crack system was modelled assuming Darcy flow in the matrix and Stokes flow in the crack. The relevance of the inertial term in the Navier–Stokes Equation (2), as discussed above, was scrutinized. A graphical representation of selected results is depicted in Figure 3, showing velocity magnitude as a colormap, as well as flow paths and contour lines of the hydraulically relevant pressure p.
The figure shows results for two angles α between crack orientation and the direction of matrix flow. For the large angle the maximum velocity magnitude that is observed in the crack doubles the mean flow velocity (that is normalized to v 0 = 1 ). For small angles the velocity within the crack reaches values that are more than 6 times higher than the mean. These differences already highlight the relevance of α that will be further scrutinized in the breakthrough curves. For large α pockets of very slow, almost stagnant flow appear near the tips of the crack. For the more upstream tip the low velocity pocket appears downstream; for the more downstream tip it appears upstream.
The figures show that the maximum velocity appears in the center of the crack, and that a boundary layer with reduced velocity is located near the interface between the porous matrix and the open crack. The results shown were obtained with consideration of the inertial term in the Navier–Stokes equations. Control model runs neglecting this term delivered output with marginal differences. As the computational costs of considering the inertial term are small, we left this option enabled for all following runs.
As outlined above the model is normalized to inflow velocity v 0 . The influence of the crack on the mean flow velocity cannot be measured directly, but can be evaluated indirectly by the pressure gradient, calculated by the model. Increased advective transport dependent on velocity is related to a decrease of the pressure gradient. As the pressure at the outlet is fixed at zero value, the computed average pressure at the inlet can be taken as representative for the overall pressure gradient. The ratio between the inlet pressures for systems without and with crack is the key measure for the increase of velocity that is induced by the crack. For a system with crack the higher velocity is represented by a small pressure gradient, so that this ratio delivers values greater than 1.
A change of mean velocity results in an equal change of advective transport through the system. Thus, we introduce an advection factor that quantifies the increase of advective transport due to the crack. Figure 4 depicts the advection factor as function of the crack angle and the crack opening. Advective transport is higher for the wider crack, as expected. It depends strongly on the angle α. If the crack is oriented perpendicular to the flow direction the factor is 1, i.e., there is no influence, while the factor increases with decreasing α. Advection may increase by 50% if the crack is aligned with the flow.

5.2. Distributions of Concentration or Temperature

The pattern of transport through the combined system of a crack located within a porous matrix may differ significantly from the homogeneous case. It is the intention here to study extreme cases, where the distribution of the transport variable clearly deviates from the regular behavior that would be observed in a porous medium without a crack.
For a homogenous porous system, the studied set-up is essentially 1D and would show a front moving along straight flow paths under the regimes of advection and diffusion. The more complex behavior induced by the crack that is of interest here and is illustrated in Figure 5. The three color-maps show the distribution of the transport variable (concentration or temperature) at time instants, when the front is reaching (a), passing through (b) and leaving the crack (c). The colors are chosen to resemble a cold front gradually replacing the warmer fluid that initially filled the entire system.
As soon as the front reaches the front end of the crack a rapid migration through the highly permeable structure can be observed (a). While the front is passing further downstream, a plume of the infiltrating fluid builds up all along beyond the downstream interface between the crack and the matrix (b). Further on, this plume extends in transversal and longitudinal direction, leading to a quick breakthrough in the center at the outflow boundary (c). Near the upper and lower boundaries however, pockets of the initial fluid remain for a long time, leading to tailing in breakthrough curves. The effective transport behavior is examined in more detail by scrutinizing the breakthrough at the outlet.
The described pattern identifies some pockets with a quick passing through the crack and thus an early arrival at downstream positions. On the other hand, other pockets are temporarily stuck in the low velocity zones and thus reach downstream positions with retardation. Altogether the breakthrough will thus show a dispersion effect, where we understand dispersion as generalized diffusive transport. The precise outlook of this dispersion is examined in the following.
In Figure 5 the flow field is illustrated by streamlines and isobars. Isobars connect the upper and lower margins of the plots. The streamlines originate at the inflow boundary on the left at equidistant positions. The flow pattern clearly indicates the deviation from the horizontal lines due to the presence of the crack. Regions with high density of the streamlines, indicating elevated velocities can be recognized at both tips of the crack. The even higher velocities within the crack are not visible because of its small aperture. Regions of low streamline density, i.e., lower flow velocity that led to pockets with retarded transport can be identified adjacent to the crack.

5.3. Breakthrough Curves

The dispersive behavior of the matrix-crack system is examined by breakthrough curves (btc). The btcs show the mean concentration at the outlet as a function of time. The outlet concentration is computed by integrating the concentration along the outlet boundary. As we deal with heat- or mass transport, the breakthrough may refer to temperature (in the case of heat transport) or to concentration (in the case of mass transport). Thus, in the visualizations of the btcs the term ‘transport variable’ is used in the label of the y-axis. In both cases the inflow value in the model is 1, i.e., the variable is normalized to the inflow value, and thus dimensionless, too.

5.3.1. Influence of the Tilt Angle

To study the influence of the tilt, angle several parametric sweeps were run. Figure 6 depicts the outflow btcs with α taking the values 0, 15, 30, 45, 60, 75 and 90°. For α = 0° the crack is aligned with the imposed pressure gradient and flow field. For α = 90° the crack is directed in a right angle with the flow.
If the crack is directed perpendicular to the flow direction, the velocity field differs only marginally from the result of a 1D plug flow model (for details see below). In that case fluxes transverse to the main flow are so small that they hardly make a difference, and the passing front remains steep. More recognizable is the slight shift of the btc towards upstream, which indicates a slightly earlier arrival time. With decreasing angle, the gradient of the btcs becomes less steep, as the figure shows. This may be conceived as an additional diffusion. As the term dispersion is usually used for generalized diffusive processes, we will use this term subsequently. For here, we can conclude that the dispersion increases, the more the crack is aligned with the flow direction.
The aim here is also to check whether the btcs from the model can be obtained by a 1D approach. For that aim we used the analytical solution of the 1D advection-diffusion equation, given by the often-cited Ogata-Banks formula [39]:
c ( x , t ) = 1 1 2 erfc x v t 2 D t + exp v x D erfc x + v t 2 D t
The diffusivity D and velocity v in Equation (10) are fitted to the btcs of the 2D model, thus obtaining values for effective parameters Deff and veff. For this task we used the ‘transportfit’ simulator [40], a MATLAB program that visualizes solutions of the 1D transport equation for advection, diffusion, decay, and retardation, based on analytical solutions (11). The program has the additional option to fit parameters to given btcs—an option that was used here.
Figure 7 depicts a typical graphical output of ‘transportfit’ for one btc that is representative for all btcs that are in question here. The given btc is represented by circles, while the colored curves provide btcs for intermediate observation points in the 1D approach. The last of the colored btcs fits perfectly with the given data, shown by circles. In most cases the effective value veff and Deff of the 1D approach deviate from the values of the original 2D system. For these cases the 1D approach can be used to reproduce the output of the more complex 2D model. This is not always the case, as will be shown below.
Table 2 lists the effective parameters obtained from the described fitting procedure. The effective Peclet number Peeff is calculated from the effective velocity and effective diffusivity. As expected, the effective velocity changes only slightly. The major change of velocity and thus of advection is included in the advection factor that was outlined above. That is due to the normalization of the model to the inflow velocity, which in the physical reality is varying.
The change of diffusivity could clearly be observed in Figure 6. The outlined procedure provides numerical values for this change. The values listed in Table 2 show that the diffusivity changes by almost one order of magnitude dependent upon the tilt angle.

5.3.2. Influence of the Crack Opening

The effect of crack opening is demonstrated by comparing btcs from a parametric sweep over d, the ratio between the long and the short axes of the elliptic crack, using values between a minimum of 0.01 and 0.2. Figure 8 depicts the resulting btcs.
The effects of the aspect ratio on the btcs are relatively small, compared to the influence of α observed above. The effects become more pronounced for highly advection-dominated flow. Figure 8 depicts the results of parametric sweeps of the crack’s axes aspect ratio for Pe = 33.3 and 333, which clearly shows the additional effect of Pe.
There is an effect on the arrival time of the mean concentration or temperature value, i.e., on advection. For wider cracks the velocities in the crack increase and lead to an earlier breakthrough. With increasing d the btcs show higher gradients in early stages, but smaller gradients in the tailing. The pockets with almost stagnant water (Figure 3) can explain the longer tailing.

5.3.3. Influence of the Peclet Number

As a strong dependency on the tilt angle α was observed, as described above, we focus here on the influence of the Pe-number on solutions for the different values of α. We also include the cases of different crack aperture, for which a weaker influence was found above. Figure 9 depicts btcs to illustrate the effects.
The Pe-number is a measure for advection dominance. Higher gradients can be expected for higher Pe. Comparing both plots in Figure 9 with different degrees of advection-dominance confirms the expectation. The plots also show that there is an effect on the arrival of the mean concentration, which is an indicator for advection. A shift to earlier arrival can be observed, which indicates an additional effect on advection. While this shift is very small for cracks that are perpendicular to flow (α = 90°), it becomes more pronounced with decreasing angles. This shift depends on α and d, but less on Pe.
For a further study of the additional effect on advection and diffusion the btcs were examined with respect to the 1D advection-diffusion model. The parameter fitting option in the ‘transportfit’ program was used for this aim. Using the btcs of the 2D model as input, diffusivity and velocity were fitted to obtain effective parameters Deff and veff. With the results an effective Pe-number Peeff can be computed using Equation (9). Results of this procedure for different values of α and the 2D model Pe are listed in Table 3.
Concerning Deff the table shows a different Pe dependence for different α. For low α (<45°) Deff increases with Pe until Pe = 16, then decreases. For high α (≥45°) Deff decreases with advection dominance. For α = 90° the effective values coincide with the 2D input values which are listed in the second row of the table. A view on the dependence on α for given Pe reveals that there is no effect for low advection-dominated flow, but a strong decrease for high Pe. The latter is in line with the former observation that the flow and transport regime is almost unaltered from 1D plug flow.
The comparison of veff with the reference value of 1, with which the 2D models were run, shows a simpler behavior: the angle dependence is marginal and concerning the Pe influence there is a monotonous decrease from the maximum of 1.33 to values ≈1. The additional effect on advection that is not reflected in the pressure increase is stronger if the advection dominance is weaker.
It is most convenient to compare the Peclet numbers Pe that are input in the 2D model with the effective one Peeff: obtained from the btcs as output from the 2D model. Similar to the observation for Deff there is a different Pe dependence depending on α. For low α (<45°) Peeff increases with Pe until Pe = 16, then decreases. For high α (≥45°) Peeff increases with Pe. The extremes are observed for high advection dominance. If the crack is aligned with flow the input and effective Pe-numbers differ by a factor of 200/2.7 = 74. If the crack is directed 90° to the flow the Pe-numbers are almost the same. Remarkably Peeff > Pe holds for slight advection-dominance. This can be attributed to the higher values of veff in comparison to the reference value v 0 = 1 .
Almost all btcs could be represented nicely by the 1D model using effective parameters. This finding is not valid for Pe > 100. For an illustration, Figure 10 depicts btcs for high Pe-numbers (Pe = 100 and 1000). While the btcs for Pe = 100 show the common shape of 1D btcs, this is not true anymore if the Pe is increased further. This is especially visible in the btc for α = 0 and Pe = 1000, which contains two regions with increased slope: for an early stage (t < 1.5) and a late stage (t > 3).

5.3.4. Further Dependencies

Both Reynolds number Re and Darcy number Da were varied over several orders of magnitude, as noted in Table 1. It is assumed that the flow within the crack is laminar. Thus, Re was not varied to values above 2000, which marks the transition to turbulence [30]. According to its definition Da depends strongly on the scale of the model, i.e., the parameter L. For L = 1 m the parameter variation covers a range of permeable porous formations. Due to very low flow velocities low permeable media are of less physical interest. The parametric sweeps with Re and Da did not reveal any influence on the results.

6. Conclusions

The transport of components in a porous matrix containing a crack was investigated systematically by numerical modelling. Converting to dimensionless parameters the analytical description could be reduced to five dimensionless parameters: the geometry parameters α and d, as well as the parameter combinations Re, Pe, and Da. Using appropriate boundary conditions, it is possible to distinguish between effects of altered flow or transport. For a wide parameter range it is possible to fit the outflow breakthrough to the btc of a 1D model approach. In that way, effective parameters could be determined for the studied setting.
Most relevant results from parameter variation runs are:
  • Strong dependence on α; the larger the angle between crack orientation and the general flow direction, the larger the effective diffusivity, the smaller the effective Peclet number.
  • There is a much smaller influence from the crack axes ratio.
  • An advection factor was introduced as a measure for the increase of velocity and advection; the maximum value of the advection factor of ≈1.5 appears when the crack is aligned with the flow direction.
  • There is an additional increase of advection that is not related to the flow field, but to the transport; it is measured by veff, which takes its highest value of 1.33 for slightly advection dominated flow.
  • For highly advection dominated flow (Pe > 100) the btc cannot always be reproduced by a 1D approach.
  • The dimensionless parameters Da and Re are not relevant.
As a dimensionless approach was used in this study the results are scale independent. However, the presented values should not be taken as valid for any constellations of a matrix-crack system. The outcome of this study shows which parameters and parameter combinations are most relevant for the changed advective and diffusive behavior, and which are not relevant. The results provide expected values, to which degree diffusion and advection change for a given parameter constellation.

Funding

This research received no external funding.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Liu, R.; Lou, S.; Jinag, Y. Recent advances in fluid flow in fractured porous media. Processes 2019, 7, 255. [Google Scholar] [CrossRef] [Green Version]
  2. National Academy of Sciences. Rock Fractures and Fluid Flow; National Academy Press: Washington, DC, USA, 1996. [Google Scholar]
  3. Dimmen, V.; Rotevatn, A.; Nixon, C.W. The relationship between fluid flow, structures, and depositional architecture in sedimentary rocks: An example-based overview. Geofluids 2020, 2020, 3506743. [Google Scholar] [CrossRef]
  4. Schultz, R. Geologic Fracture Mechanics; Cambridge University Press: Cambridge, UK, 2019. [Google Scholar]
  5. Broberg, K. Cracks and Fracture; Academic Press: San Diego, CA, USA, 1999. [Google Scholar]
  6. Shao, Q.; Fernández-González, R.; Mikdam, A.; Bouhala, L.; Younes, A.; Núñez, P.; Belouettar, S.; Makradi, A. Influence of heat transfer and fluid flow on crack growth in multilayered porous/dense materials using XFEM: Application to solid oxide fuel cell-like material design. Int. J. Solids Struct. 2014, 51, 3557–3569. [Google Scholar] [CrossRef] [Green Version]
  7. Tsushima, S.; Hirai, S. An overview of cracks and interfacial voids in membrane electrode assemblies in polymer electrolyte fuel cells. J. Therm. Sci. Technol. 2015, 10, JTST0002. [Google Scholar] [CrossRef] [Green Version]
  8. Kumano, N.; Kudo, K.; Suda, A.; Akimoto, Y.; Ishii, M.; Nakamura, H. Controlling cracking formation in fuel cell catalyst layers. J. Power Sources 2019, 419, 219–228. [Google Scholar] [CrossRef]
  9. Berre, I.; Doster, F.; Keilegavlen, E. Flow in fractured porous media: A review of conceptual models and discretization approaches. Transp. Porous Media 2019, 130, 215–236. [Google Scholar] [CrossRef] [Green Version]
  10. Pouya, A.; Ghabezloo, S. Flow around a crack in a porous matrix and related problems. Transp. Porous Media 2010, 84, 511–532. [Google Scholar] [CrossRef]
  11. Sahimi, M. Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches; Wiley & Sons: Hoboken, NJ, USA, 2011. [Google Scholar]
  12. Bodin, J.; Delay, F.; de Marsily, G. Solute transport in a single fracture with negligible matrix permeability: 1. Fundamental mechanisms. Hydrogeol. J. 2003, 11, 418–433. [Google Scholar] [CrossRef]
  13. Bodin, J.; Delay, F.; de Marsily, G. Solute transport in a single fracture with negligible matrix permeability: 2. Mathematical formalism. Hydrogeol. J. 2003, 11, 434–454. [Google Scholar] [CrossRef]
  14. Liu, L.; Neretnieks, I.; Shahkarami, P.; Meng, S. Solute transport along a single fracture in a porous rock: A simple analytical solution and its extension for modeling velocity dispersion. Hydrogeol. J. 2018, 26, 297–320. [Google Scholar] [CrossRef]
  15. Kuva, J.; Voutilainen, M.; Mattila, K. Modeling mass transfer in fracture flows with the time domain-random walk method. Comput. Geosci. 2019, 23, 953–967. [Google Scholar] [CrossRef] [Green Version]
  16. Natarajan, N.; Suresh Kumar, G. Numerical modeling of solute transport in a coupled sinusoidal fracture matrix system in the presence of fracture skin. Int. J. Energy Environ. 2010, 4, 99–104. [Google Scholar]
  17. Mehmani, A.; Kelly, S.; Torres-Verdin, C.; Balhoff, M. Quantification of fracture-matrix fluid transport in unconventional rocks using two-scale microfluidic chips. In Proceedings of the Unconventional Resources Technology Conference, Austin, TX, USA, 24–27 July 2017. [Google Scholar]
  18. Bagalkot, N.; Suresh Kumar, G. Colloid transport in a single fracture–matrix system: Gravity effects, influence of colloid size and density. Water 2018, 10, 1531. [Google Scholar] [CrossRef] [Green Version]
  19. Zhu, Y.; Zhan, H.; Jin, M. Analytical solutions of solute transport in a fracture–matrix system with different reaction rates for fracture and matrix. J. Hydrol. 2016, 539, 447–456. [Google Scholar] [CrossRef]
  20. Lonergan, L.; Wolkinson, J.; McCaffrey, K. Fractures, Fluid Flow and Mineralization; Special Publ. 155; Geological Society: London, UK, 1999. [Google Scholar]
  21. Yuan, T.; Wei, C.; Zhang, C.-S.; Qin, G. A numerical simulator for modeling the coupling processes of subsurface fluid flow and reactive transport processes in fractured carbonate rocks. Water 2019, 11, 1957. [Google Scholar] [CrossRef] [Green Version]
  22. Diodato, D.M. A Compendium of Fracture Flow Models; Argonne National Laboratory, Center for Environmental Restoration Systems, Energy Systems Division: Chicago, IL, USA, 1994. [Google Scholar]
  23. Adler, P.M.; Thovert, J.-F.; Mourzenko, V.V. Fractured Porous Media; Oxford University Press: Oxford, UK, 2013. [Google Scholar]
  24. Schwenck, N.; Flemisch, B.; Helmig, R.; Wohlmuth, B.L. Dimensionally reduced flow models in fractured porous media: Crossings and boundaries. Comput. Geosci. 2015, 19, 1219–1230. [Google Scholar] [CrossRef]
  25. Holzbecher, E. Benchmarking model approaches for thin structures in a porous matrix. Int. J. Multiphysics 2020, 14, 273–288. [Google Scholar]
  26. Wang, H.; Yi, S.; Mukul Sharma, M. A computationally efficient approach to modeling contact problems and fracture closure using superposition method. Theor. Appl. Fract. Mech. 2018, 93, 276–287. [Google Scholar] [CrossRef] [Green Version]
  27. De Borst, R. Fluid flow in fractured and fracturing porous media: A unified view. Mech. Res. Commun. 2017, 80, 47–57. [Google Scholar] [CrossRef]
  28. Langtangen, H.-P.; Pedersen, G.K. Scaling of Differential Equations; Springer International Publishing: Berlin/Heidelberg, Germany, 2016. [Google Scholar]
  29. Yarin, L.P. The Pi-Theorem; Springer Publishing: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
  30. Sonin, A.A. The Physical Basis of Dimensional Analysis; MIT, Department of Mechanical Engineering: Cambridge, UK, 2001. [Google Scholar]
  31. Guyon, E.; Hulin, J.-P.; Petit, L. Hydrodynamik; Vieweg: Braunschweig, Germany, 1997. [Google Scholar]
  32. Zimmerman, R.W.; Bodvarsson, G.S. Hydraulic conductivity of rock fractures. Transp. Porous Media 1996, 23, 1–30. [Google Scholar] [CrossRef] [Green Version]
  33. Holzbecher, E. Environmental Modeling; Springer Publishing: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
  34. Holzbecher, E.; Wong, L.; Litz, M.-S. Modelling flow through fractures in porous media. In Proceedings of the COMSOL Conference, Paris, France, 17–19 November 2010. [Google Scholar]
  35. Romano-Perez, C.A.; Diaz-Viera, M.A. A comparison of discrete fracture models for single phase flow in porous media by COMSOL Multiphysics® Software. In Proceedings of the 2015 COMSOL Conference, Boston, MA, USA, 7–9 October 2015. [Google Scholar]
  36. Perko, J.; Seetharam, S.; Mallants, D. Verification and validation of flow and transport in cracked saturated porous media. In Proceedings of the COMSOL Conference, Stuttgart, Germany, 26–28 October 2011. [Google Scholar]
  37. Perko, J.; Seetharam, S.; Diederik, J.; Mallants, D.; Cool, W.; Vermarien, E. Influence of cracks in cementitious engineering barriers in a near-surface disposal system: Assessment analysis of the Belgian case 2013. In Proceedings of the 15th International Conference on Radioactive Waste Management and Environmental Remediation, Brussels, Belgium, 8–12 September 2013. [Google Scholar]
  38. Bruining, H.; Darwish, M.; Rijnks, A. Computation of the longitudinal and transverse dispersion coefficient in an adsorbing porous medium using homogenization. Transp. Porous Media 2012, 91, 833–859. [Google Scholar] [CrossRef] [Green Version]
  39. Ogata, A.; Banks, R.B. A Solution of the Differential Equation of Longitudinal Dispersion in Porous Media; Professional Paper, No. 411-A; United States Geological Survey: Reston, VA, USA, 1961.
  40. Holzbecher, E. Transportfit.m, 2021, Github. Available online: https://github.com/eholzbe/1D-transport-visualization (accessed on 21 September 2021).
Figure 1. Geometry of a region with a crack of elliptic shape in a porous medium, here the crack is tilted by an angle of 45° to the flow direction.
Figure 1. Geometry of a region with a crack of elliptic shape in a porous medium, here the crack is tilted by an angle of 45° to the flow direction.
Fluids 06 00358 g001
Figure 2. Finite element mesh with elliptic fracture; (a): complete model region, (b): detail.
Figure 2. Finite element mesh with elliptic fracture; (a): complete model region, (b): detail.
Fluids 06 00358 g002
Figure 3. Flow patterns in porous matrix and crack for tilting angles α = 15 (a) and 75° (b); color plot of velocity magnitude, flow paths (white with arrows) and isobars (white without arrows) for parameter values d = 0.1, Re = 10, Da = 10−11.
Figure 3. Flow patterns in porous matrix and crack for tilting angles α = 15 (a) and 75° (b); color plot of velocity magnitude, flow paths (white with arrows) and isobars (white without arrows) for parameter values d = 0.1, Re = 10, Da = 10−11.
Fluids 06 00358 g003
Figure 4. Advection factor as function of the crack angle for d = 0.1 and d = 0.01; Re = 10, Da = 10−11.
Figure 4. Advection factor as function of the crack angle for d = 0.1 and d = 0.01; Re = 10, Da = 10−11.
Fluids 06 00358 g004
Figure 5. Temperature distributions at three-time instances for a cold front entering from the left, replacing hot fluid; showing the distortion of the front while reaching (a), passing through (b) and leaving the crack (c) (see text for details).
Figure 5. Temperature distributions at three-time instances for a cold front entering from the left, replacing hot fluid; showing the distortion of the front while reaching (a), passing through (b) and leaving the crack (c) (see text for details).
Fluids 06 00358 g005aFluids 06 00358 g005b
Figure 6. Breakthrough curves for different tilting angles α, d = 0.05, Pe = 33.3, Re = 1, Da = 10−14.
Figure 6. Breakthrough curves for different tilting angles α, d = 0.05, Pe = 33.3, Re = 1, Da = 10−14.
Fluids 06 00358 g006
Figure 7. Example of fitting breakthrough curves from 2D models by 1D advection-diffusion model, d = 0.05, α = 0°, Pe = 33, Re = 1, Da = 10−14 (see text for details).
Figure 7. Example of fitting breakthrough curves from 2D models by 1D advection-diffusion model, d = 0.05, α = 0°, Pe = 33, Re = 1, Da = 10−14 (see text for details).
Fluids 06 00358 g007
Figure 8. Breakthrough curves as function of aspect ratio d; for cracks with α = 45°, Re = 1, Da = 10−14; (a) for Pe = 33.3, (b) for Pe = 333.
Figure 8. Breakthrough curves as function of aspect ratio d; for cracks with α = 45°, Re = 1, Da = 10−14; (a) for Pe = 33.3, (b) for Pe = 333.
Fluids 06 00358 g008aFluids 06 00358 g008b
Figure 9. Breakthrough curves for varying α and d, Pe = 100 (a) and Pe = 333 (b), Re = 1, Da = 10−14.
Figure 9. Breakthrough curves for varying α and d, Pe = 100 (a) and Pe = 333 (b), Re = 1, Da = 10−14.
Fluids 06 00358 g009aFluids 06 00358 g009b
Figure 10. Breakthrough curves for varying α, d = 0.05, Pe = 100 and 1000, Re = 1, Da = 10−14.
Figure 10. Breakthrough curves for varying α, d = 0.05, Pe = 100 and 1000, Re = 1, Da = 10−14.
Fluids 06 00358 g010
Table 1. Parameters, their reference values, and variations in parametric sweeps.
Table 1. Parameters, their reference values, and variations in parametric sweeps.
ParameterUnitSymbolReference ValueVariation
Model height/widthmH1/
Model lengthm 2/
Crack lengthmL1/
Crack, ratio of half-axes/d0.050.01, 0.05, 0.1
Crack tilting angle/α45°0, 7.5, 15, 22.5, 30, 45, 60, 67.5, 75, 82.5, 90°
Matrix permeabilitym2k10−14/
Fluid viscosityPa sμ0.001/
Fluid densitykg/m3ρ1000/
Inflow velocitym/s v 0 10−6/
Heat capacity ratio/κ1/
Diffusivitym2/sD3 × 10−8/
Reynolds number/Re10.1, 1, 10, 100, 1000
Darcy number/Da10−1410−14, 10−13, 10−12, 10−11, 10−10
Peclet number/Pe33.32, 3.3, 4, 8, 10, 16, 33.3, 64, 100, 200, 333, 1000
Table 2. Effective diffusivities dependent upon the tilt angle α for d = 0.05, Pe = 33.3, Re = 1, Da = 10−14.
Table 2. Effective diffusivities dependent upon the tilt angle α for d = 0.05, Pe = 33.3, Re = 1, Da = 10−14.
α15°30°45°60°75°90°
Deff0.2260.1980.1370.07770.03390.03120.0297
veff1.0051.0051.0061.0101.0141.0141.014
Peeff4.4505.0767.34313.029.9132.534.1
Table 3. Effective parameters for 2D-model Pe-numbers and selected values of α, d = 0.05, Re = 1, Da = 10−14.
Table 3. Effective parameters for 2D-model Pe-numbers and selected values of α, d = 0.05, Re = 1, Da = 10−14.
Pe24810163364100200
D0.50.250.1250.10.06250.03030.01560.010.005
Effective diffusivity Deff
α = 0°0.5170.3020.2110.20.1950.2260.2710.3040.357
α = 22.5°0.5120.2910.1910.1750.1610.1690.1890.2040.225
α = 45°0.50.2680.1520.1310.10.07770.0690.0650.0614
α = 67.5°0.4940.2530.1290.1040.0670.0350.0200.01460.00923
α = 90°0.4930.2490.1230.10.0620.02970.01560.010.00525
Pe24810163364100200
Effective velocity veff
α = 0°1.331.1461.0641.0491.0271.0050.9890.9790.964
α = 22.5°1.331.1441.0631.0481.0251.0050.9940.990.986
α = 45°1.331.1431.0641.0491.0291.011.0041.0031.006
α = 67.5°1.331.1431.0661.0521.0311.0141.0081.0051.005
α = 90°1.321.1441.0661.0521.0311.0141.0061.0031.002
Pe24810163364100200
Effective Peclet number Peeff
α = 0°2.573.795.045.255.274.453.653.222.7
α = 22.5°2.63.935.576.06.375.955.264.854.38
α = 45°2.664.267.08.0110.2913.014.615.4316.4
α = 67.5°2.694.528.2610.1215.3929.050.468.84109
α = 90°2.684.598.6710.516.6334.164.5100.3191
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Holzbecher, E. Breakthrough Investigation of Advective and Diffusive Transport in a Porous Matrix with a Crack. Fluids 2021, 6, 358. https://doi.org/10.3390/fluids6100358

AMA Style

Holzbecher E. Breakthrough Investigation of Advective and Diffusive Transport in a Porous Matrix with a Crack. Fluids. 2021; 6(10):358. https://doi.org/10.3390/fluids6100358

Chicago/Turabian Style

Holzbecher, Ekkehard. 2021. "Breakthrough Investigation of Advective and Diffusive Transport in a Porous Matrix with a Crack" Fluids 6, no. 10: 358. https://doi.org/10.3390/fluids6100358

APA Style

Holzbecher, E. (2021). Breakthrough Investigation of Advective and Diffusive Transport in a Porous Matrix with a Crack. Fluids, 6(10), 358. https://doi.org/10.3390/fluids6100358

Article Metrics

Back to TopTop