Free convection in fractured porous media: a numerical study (2024)

1 Introduction

The study of flow in porous media finds application in many different areas ranging from industrial, to biomedical and environmental applications.In the context of geological porous media many relevant problems such as geothermal engineering and contaminant transport are nowadays addressed with the help of numerical simulations, due to the necessity of estimating flow rates and pathways.

In many cases, due to either temperature or the presence of solutes, fluids can present density variations in space and time.Temperature gradients are particularly important in geothermal applications,while solute concentration gradients are found for instance in the study of contaminant plume migration or seawater intrusion phenomena.

In situations of inverse density gradients (denser fluid on top) the fluid may spontaneously develop unstable convective plumes.It is particularly important to predict the possible onset of such instabilities due to the strong rates at which convection can transport solutes with respect to diffusion alone.

The first studies dealing with free convection in porous media [6, 4] examined simple scenarios of inverse density gradients. The possibility of convection is linked to the Rayleigh number, which, accounting for both fluid and porous medium parameters, quantifies the strength of convection against that of diffusion.For large values of the Rayleigh number, there is the possibility of free convection.

In recent years, different authors have tried to understand how the presence of fractures influences the possibility and strength of free convection [10, 9, 17, 16, 11].While high-density fracture networks can be mostly treated through averaging/upscaling of porous media properties, low-density fracture networks can manifest unexpectedly high convection strength due to the particular geometry of the network.Even for simple geometries, e.g. for regular horizontal or vertical fracture grids, the unstable nature of convective plumes makes prediction very hard if not impossible [10].

The main focus of this study is the development of a new method based on the eigenvalue analysis of equilibrium solutions for assessing the possibility of convection in fractured porous media. Leveraging this new method will enable us to better understand the particular mechanisms by which fractures enhance convection. This idea is at the core of the results provided by [6], where the analysis is carried out analytically for homogenoues media, and in [8] for layered porous media. Here, given the geometrical complexity of the domain, eigenvalues will be computed by a suitable numerical algorithm starting from the discretized problem.

We start in section 2 with a detailed description of the mathematical model both in homogeneous and fractured media.As the focus of this study is the effect of fracture geometry, simple constitutive models are preferred over more elaborate (and more accurate) models.The Darcy law is used as constitutive law for filtration, Fick’s law is used to model solute diffusion and density is modeled as a linear function of solute concentration.

Particular care is given to the averaging procedure used for the dimensional reduction of fractures, consistent with the mixed-dimensional approach presented in [7].

In section 3 we present the discretization of the continuous model using the finite volume method in space and the implicit Euler scheme in time.The particular finite volume method used in the implementation is the Multipoint Flux Approximation method [1] for both the fluid mass conservation and the solute transport problems.The MPFA method is particularly appropriate due to its mass conservation properties and consistency for general grids.Again, density dependence must be treated with care for the numerical implementation to be consistent [12].

Two methods will be described for assessing the possibility of convective cell formation.

Section 3.1 describes what we call the direct method of assessing stability: starting from a non-equilibrium initial condition and integrating the time-dependent flow equations until steady state may result in either convective motion or a diffusive equilibrium solution. The two regimes are easily distinguished by measuring the amount of solute transport through the domain.

Section 3.2 describes the eigenvalue method for assessing stability. By linearizing the problem and numerically studying the eigenvalues of the resulting discrete system, we can assess the stability of arbitrary perturbations to any equilibrium solution without the need of solving the time-dependent problem.

Sections 4.1 and 4.2 test the two methods respectively with the Elder and HRL problems relying on the results of [17] and [2]: both are well known benchmark cases for density driven flow.Section 4.4 validates and examines a three-dimensional generalization of the HRL problem, based on the studies in [11, 16].In the HRL scenarios, the complementarity of the two approaches described above will enable us to better understand and examine the peculiarities of free convection in the presence of fractures.

The variety of scenarios is also useful for understanding the advantages and limitations of the eigenvalue method with respect to the direct method. A comparison of their computational cost is presented in section 4.5.

Finally, section 5 will conclude this study with a critical evaluation of what it has accomplished and suggest different options for further developments.

2 Mathematical model

The problem of our interest stems from the coupling between solute transport and density driven flow in a porous medium.The mathematical model is based on the one described in the review paper by Diersch and Kolditz [2], which will be here extended to account for the presence of fractures in the domain, since our main goal is to understand the impact of fracture networks on the onset of free convection.

where $\rho\;[\mathrm{M}/\mathrm{L}_{v}^{3}]$ is the density of the fluid in which the solute is dissolved, $\omega\;[\mathrm{M_{s}}/\mathrm{M}]$ is the concentration of the solute, $\bm{u}\;[\mathrm{L}/\mathrm{T}]$ is the flow velocity, $\rho\phi\bm{i}\;[\frac{\mathrm{M_{s}}}{\mathrm{L}^{2}\mathrm{T}}]$ is the solute mass flux due to diffusion and $\phi\;[\mathrm{L_{v}^{3}}/\mathrm{L^{3}}]$ is the porosity of the medium.The subscript $v$ in the length dimension $L$ is used to differentiate void volumes from total volumes, while $s$ is used to differentiate solute mass from fluid mass.

where $\bm{g}=-g\bm{e}_{z}$ is the gravity acceleration vector, $k$ is the permeability tensor of the porous medium $k\;[\mathrm{L}_{v}^{2}]$, and $\mu\;[\mathrm{M}/\mathrm{L}_{v}\mathrm{T}]$ is the viscosity of the fluid.

The boundary condition for (4) can beof Dirichlet type where the pressure is specified: $p=p_{D}$or of Neumann type where the normal flow velocity is prescribed, $\bm{u}\cdot\bm{n}=u_{N}$.

A large part of the following discussion applies unchanged to flows driven by temperature gradients instead of solute concentration gradients.Indeed, apart from a reinterpretation of some of the quantities introduced above e.g. thermal conductivity instead of diffusivity,the same model can be applied.

2.2 Horton-Rogers-Lapwood (HRL) problem

The HRL problem is a simple scenario aiming to study the possible onset and strength of natural convection.The idea is to impose, through the boundary conditions, a layer of heavier fluid overlaid on lighter fluid in a two-dimensional vertical cross section of a homogeneous porous medium.The density gradient can be caused by temperature difference (as in the original description [6]), or by solute concentration difference.In a situation of inverse mass gradient, the fluid will form convective cells only if diffusivity is small enough to allow it.The contrast between convection and diffusion speed, respectively $v_{g}$ and $v_{d}$, is described by the Rayleigh number

$$\mathit{Ra}=\frac{v_{g}}{v_{d}}=\frac{\frac{k}{\phi\mu}\rho_{0}\alpha\omega_{$$

(8)

where $H$ is the height of the domain, and other quantities have been introduced in the previous section.In the original presentation of this problem [6], the possibility of convective motion is linked to the value of $\mathit{Ra}$: $\mathit{Ra}_{c}=4\pi^{2}$ is the critical number i.e. a threshold such that, for $Ra>Ra_{c}$ convection occurs.Why is free convection only possible for such values of the Rayleigh number?For high enough diffusion, convective cells cannot sustain the concentration/temperature difference and thus they simply decay.Figure 1 illustrates how the diffusive flow tries to restore concentration imbalances caused by the convective motion.

Note that convective motion in general enhances solute transport.For this reason, an indicator of the presence of convection is the Sherwood number $\mathit{Sh}$, defined as

$$\mathit{Sh}=\frac{\int_{A}\bm{i}\cdot\bm{n}}{\int_{A}\bm{i}_{0}\cdot\bm{n}}=$$

(9)

where $A$ is the diffusive inflow surface (top boundary) and $\bm{i}_{0}$ is the diffusive flow in the absence of convection. When convection is not present inside the domain, $\bm{i}=\bm{i_{0}}$ and thus $\mathit{Sh}=1$. Note that the boundary conditions for the HRL problem being of zero fluid mass flow all around, the solute transport on the inflow surface $A$ is entirely due to diffusion, regardless of whether convection is present inside the domain. However, when convection is present, the (overall) larger concentration gradient at $A$ makes $\bm{i}>\bm{i_{0}}$ and $\mathit{Sh}>1$.

2.3 Dimensional reduction

Porous media often present heterogeneities in their material properties.A particularly strong kind of heterogeneity are fractures: regions of different material properties, with negligible aperture (thickness) with respect to both their length and the characteristic lengths of the medium.Very different material properties in fractures, such as permeability, can compensate for their small dimensions so that, overall, fractures can strongly influence the behaviour of flow in the medium. Let us define

See Also

Improving Dense Retriever Training with Corrector Networks¿Qué Significa Acento Diacrítico?Learning in Hybrid Active Inference Models

$$k=k(\bm{x})=\left\{\begin{aligned} &k_{b}\enspace&&\bm{x}\in\Omega_{b},\\&k_{f}&&\bm{x}\in\Omega_{f},\end{aligned}\right.$$

where $\Omega_{f}$ denotes the fracture, $\Omega_{b}$ the bulk medium and $\Omega=\Omega_{b}\cup\Omega_{f}$.

The approach known as discrete fracture modeling aims to solve this problem by treating fractures as separate domains, often of reduced dimensionality, coupled to the bulk medium by coupling conditions at their interface.

We will now rewrite the fluid mass conservation equation and associated Darcy law as a mixed-dimensional equation.To keep our derivation simple, we will treat a single two-dimensional fracture as illustrated in figure 2. The domain $\Omega$ is split between bulk medium and fracture $\Omega=\Omega_{b}\cup\Omega_{f}\cup\Gamma$, where $\Gamma=\overline{\Omega_{b}}\cap\overline{\Omega_{f}}$ denotes the interface between subdomains. We assume that $\Omega_{f}$ can be expressed as

	$\displaystyle\Omega_{f}=\{\bm{x}\in\Omega:\bm{x}=\bm{\gamma}+\alpha\bm{n},\>$
	$\displaystyle\Gamma^{0}=\{\bm{x}\in\Omega:\bm{x}=\bm{x}_{0}+s_{1}(\bm{x}_{1}-$

where we have assumed a planar midsurface $\Gamma^{0}$ as shown in figure 3. The lateral boundary of $\Omega_{f}$ will be treated differently from the top and bottom boundaries. We define

	$\displaystyle\Gamma^{\pm}=\{\bm{x}\in\partial\Omega_{f}:\bm{x}=\bm{\gamma}\pm$
	$\displaystyle\Sigma=\{\bm{x}\in\partial\Omega_{f}:\bm{x}=\bm{\sigma}+\alpha\bm$

such that $\Gamma=\Gamma^{\pm}\cup\Sigma$.Superscripts + and - will also be used to denote quantities evaluated on $\Gamma^{\pm}$.

Starting from the continuity equation and the associated Darcy law in the two domains, we enforce pressure continuity and conservation of mass across the interface $\Gamma$:

$$\left\{\begin{aligned} &\nabla\cdot\bm{u}_{b}=0&&\text{in }\Omega_{b},\\&\bm{u}_{b}=k_{b}/(\phi\mu)\;(-\nabla p_{b}+\rho_{0}\alpha\omega_{b}\bm{g})&$$

(10)

where $\bm{n}_{b}$ is the normal vector defined on $\Gamma$, exiting from $\Omega_{b}$.At the end of dimensional reduction, the equations solved on the three-dimensional domain $\Omega_{f}$ will be replaced by the solution of (different) equations on the two-dimensional domain $\Gamma^{0}$.For the sake of simplicity, we will suppose that fractures differ from the bulk only in their permeability: all the other material parameters will be common to both.

We proceed by integrating the mass conservation equation in the fracture across the aperture.

$$0=\int_{-b/2}^{b/2}\nabla\cdot\bm{u_{f}}=\int_{-b/2}^{b/2}\nabla\cdot(\mathrm{$$

where we have denoted by $\mathrm{T}$, $\mathrm{N}$ the projection operators specific to the fracture:$\mathrm{N}\bm{v}=(\bm{v}\cdot\bm{n})\bm{n}$ and$\mathrm{T}\bm{v}=(\mathrm{I}-\mathrm{N})\bm{v}$, where $\bm{n}$ is the normal vector to the fracture plane (as indicated in figure 3), and $\mathrm{I}$ is the identity operator.We also denoted the in-plane gradient operator by $\nabla_{\tau}=\mathrm{T}\nabla$, and by $\bm{u}_{\gamma}$ the integral of the tangential flow field:

$$\bm{u}_{\gamma}=\int_{-b/2}^{b/2}\mathrm{T}\bm{u_{f}},\quad[\mathrm{L}^{2}/$$

Indeed, given these definitions,

$$\int_{-b/2}^{b/2}\nabla\cdot(\mathrm{T}\bm{u_{f}})=\int_{-b/2}^{b/2}\nabla_{$$

To obtain a law for $\bm{u}_{\gamma}$, we integrate the in-plane component of the Darcy law:

$\displaystyle\bm{u_{\gamma}}=\int_{-b/2}^{b/2}\mathrm{T}\bm{u_{f}}=\frac{k_{f}$

where $p_{\gamma}=\frac{1}{b}\int_{-b/2}^{b/2}p_{f}$, $\omega_{\gamma}=\frac{1}{b}\int_{-b/2}^{b/2}\omega_{f}$ and $\bm{g_{\tau}}=\mathrm{T}\bm{g}$. Differently from $\bm{u}_{\gamma}$, scalar quantities in the fracture are averaged across the fractures thus maintaining the same dimensions as the corresponding bulk quantities.Coupling conditions at the interface must also be expressed in terms of the averaged variables:

$$(\operatorname{tr}\bm{u}_{b}\cdot\bm{n}_{b})^{\pm}\approx\frac{bk_{f}}{\phi\mu$$

We conclude by replacing $\Omega_{f}$ with its center plane $\Gamma^{0}$ and extending $\Omega_{b}$: $\Omega=\Omega_{b}\cup\Omega_{f}\cup\Gamma\approx(\Omega\setminus\Gamma^{0})$. Even though the fracture domain $\Omega_{f}$ and the interface $\Gamma$ have collapsed on one another geometrically, the two play distinct roles and must be kept conceptually separate.

Collecting the last steps, we can write the Darcy part of the mixed-dimensional system:

$$\left\{\begin{aligned} &\nabla\cdot\bm{u}_{b}=0&&\text{in }\Omega_{b},\\&\bm{u}_{b}=k_{b}/(\phi\mu)\;(-\nabla p_{b}+\rho_{0}\alpha\omega_{b}\bm{g})&&$$

(11)

In writing (11), we have introduced the new variable $\lambda\,\mathrm{[L/T]}$ and the jump operator $[v]=v^{+}+v^{-}$.For $v$ defined on $\Gamma^{\pm}$, we will consider $[v]$ to be defined on $\Gamma^{0}$.Note that $\Sigma$ collapsed onto a lower-dimensional object (the boundary of the center plane $\partial\Gamma^{0}$) onto which we will have to impose boundary conditions. We will ignore its contribution to the mass exchange between fracture and bulk by setting zero normal flux (note that this flux scales linearly with the fracture aperture $b$) on immersed fracture boundaries (or tips) whereas fractures will inherit boundary conditions from the bulk if they touch the boundary.

Note that two sources of modeling error associated with dimensional reduction: (i) by collapsing the thin dimension of the fracture we have reassigned part of the domain which previously belonged to the fracture to the bulk medium and (ii) the flux $\lambda$ exchanged between the bulk medium and the fracture is a first order approximation to the true flux due to the presence of the normal gradient of the pressure $\nabla p$.

2.4 Mixed-dimensional transport

We can now re-apply the reduction procedure to the transport problem. We start by splitting the equation in the two domains and prescribing compatibility conditions at the interface

$$\left\{\begin{aligned} &\partial_{t}\omega_{b}+\nabla\cdot\bm{q_{b}}=0&&\text{$$

(12)

We integrate the conservation equation over the fracture, splitting the fluxes into their normal and tangential parts:

$\displaystyle 0=\int_{-b/2}^{b/2}\partial_{t}\omega_{f}+\nabla\cdot\bm{q_{f}}=$

As we did for the fluid mass conservation equation, we introduced fracture quantities $\omega_{\gamma}=\frac{1}{b}\int_{-b/2}^{b/2}\omega_{f}$ and the integrated tangential solute mass flux $q_{\gamma}=\int_{-b/2}^{b/2}\mathrm{T}\bm{q}_{f}$. As done before, scalar quantities are averaged while vector quantities are integrated.

Unlike the continuity equation, a new approximation is needed in the averaging step to deal with the nonlinearity of the advective term: let us consider a splitting of $\omega_{f}$ into $\omega_{\gamma}+\tilde{\omega}_{f}$, where $\omega_{\gamma}$ is constant across the fracture and $\tilde{\omega}_{f}$ is a null average fluctuation, and similarly for $\bm{u}_{f}$. We will neglect the product of fluctuations assuming that they are small, i.e.

$\displaystyle\int_{-b/2}^{b/2}\nabla\cdot(\omega_{f}\bm{u}_{f})\approx\omega_{$

For fractures of non-negligible aperture however, internal motions (see e.g. intrafractrure mode 2A in figure 23) could make this approximation inappropriate.

Just as for Darcy’s law, Fick’s law will look identical when projected in the fracture plane, while for the normal projection, we will have to resort to a first order approximation

$$(\operatorname{tr}\bm{q}_{b}\cdot\bm{n}_{b})^{\pm}=D(\frac{\omega_{b}^{\pm}-$$

Finally, we obtain a mixed-dimensional system for the transport of a solute in a fractured porous medium,

$$\left\{\begin{aligned} &\partial_{t}\omega_{b}+\nabla\cdot\bm{q}_{b}=0&&\text{$$

(13)

to be complemented with boundary and initial conditions.It is important to notice that the equations in the bulk and the equations in the fracture, both in (11) and (13), are entirely decoupled apart from their interaction through the interface variables $\lambda$ and $\theta$.The way these variables appear in the equations reveals how the domains are coupled across different dimensions: in the higher-dimensional domain interface variables appear as Neumann boundary conditions, while in the lower-dimensional fracture they appear as sources inside the domain.Although the averaging procedure has been carried out for a three-dimensional domain with a single fracture,the exact same procedure works for a generic $n$-dimensional domain ($n=1,2,3$) with multiple, possiblyintersecting fractures. In this case, the procedure is conceptually carried out hierarchically: $n$-dimensionalquantities are coupled to $(n-1)$-dimensional quantities through $(n-1)$-dimensional interface fluxes (figure 4).

In view of dealing with the general case of multiple fractures of different dimensions, we want to write equations which hold for each dimension.Apart from making the system of equations more compact, it will also make the mechanism by which different domains interact more clear.

We collect all domains of equal dimension $d\in\{0,1,2,3\}$ in a single domain $\Omega_{d}$, and denote by $\Gamma_{d}$ the interface between domains $\Omega_{d}$ and $\Omega_{d+1}$. We also introduce mixed-dimensional variables $\omega_{d},p_{d}$, fluxes $\bm{u}_{d},\bm{i}_{d}$ and in-plane gradient operator $\nabla_{d}$ defined on $\Omega_{d}$, interface fluxes $\lambda_{d},\theta_{d}$ defined on $\Gamma_{d}$.

With these quantities at hand, we can generalize systems (11, 13) to

$$\left\{\begin{aligned} &\nabla_{d}\cdot\bm{u}_{d}=[\lambda_{d}]&&\text{in }$$

(14)

with suitable boundary conditions.

2.5 The impact of fractures on convection onset

The original discussion of the HRL problem addressed the question of the possible onset of convection.This was expressed as a critical Rayleigh number $\mathit{Ra}_{c}$, such that convective motion is possible for $\mathit{Ra}\geq\mathit{Ra}_{c}$.

In recent years, different studies have tried to understand in what way fractures influence the possibility of convection.It is clear that highly permeable fractures constitute preferential paths for flow, thus enabling convection or enhancing its strength.As shown in [17], for large fracture density, calculating an average Rayleigh number based on the average (upscaled) permeability (neglecting the specific fracture configuration) can be effective at predicting the onset and strength of convection.

For lower fracture densities however, this approach is not adequate.As shown in figure 5, a continuous fracture loop barely modifying the permeability, such that the average Rayleigh number is well below the critical Rayleigh number for homogeneous media, still exhibits convective motion.In [17], the key factor for enabling convection in the case of low-density fracture configurations is shown to be the presence of continuous fracture circuits, around which convection cells can form.Simple scenarios (also used as validation cases in this work) tried to relate the location, aspect ratio and size of fracture circuits to the possibility and strength of convection.

A more systematic study remains to be done, with the aim of uncovering better quantitative relations that can be extended to more complex fracture configurations. Moreover, we have to consider that in applicative scenarios the particular geometry of the underground fracture network is mostly unknown, or in the best cases described by statistical parameters such as fracture density, mean lengths and orientations.In these cases, being able to relate statistical parameters such as the ones mentioned to a quantitative estimate on the strength of convection and its uncertainty would be both useful from an applicative standpoint and interesting in its own right.

Remark 1.

Let us consider the diffusion of a solute in a domain cut by a horizontal fracture.The continuous equidimensional problem (12) admits the linear concentration profile solution for the boundary conditions prescribed by the HRL problem.

$$\omega_{j}=\omega_{max}\frac{y}{H},\quad\bm{i}_{j}=-D\nabla\omega_{j}=-D\frac{$$

for $j\in\{B,f\}$.

With dimensional reduction, the system of equations to be solved is replaced by (13), which yields a piecewise linear concentration profile:

	$\displaystyle\omega_{b}=\left\{\begin{aligned} &(\omega_{max}-\delta\omega)\,{$
	$\displaystyle\omega_{f}=\omega_{max}/2$
	$\displaystyle\delta\omega=\frac{\omega_{max}}{1+H/b},$

$b$ being the fracture aperture. Fractures which qualify as thin enough to be treated as lower-dimensional regions will always satisfy $b\ll H$,thus making the error in the concentration profile small: $\delta\omega\ll\omega_{max}$. This small model reduction error however can manifest itself in the form of (small) artificial fluxes around fracture tips: the concentration gradient that arises from matching the two solutions creates a circulating diffusive flux, as illustrated in figure 6 on the left. The exact same reasoning can then be applied for the fluid mass conservation equation: small artificial fluid mass fluxes may appear around fracture tips due to the discontinuity of pressure across the fracture (see figure 6 on the right).

3 Numerical discretization

This section is dedicated to the discretization of system (14).The method chosen for the spatial discretization of the system is the finite volume method.In the fractured problem, we have a sequence of domains $\Omega_{d},\,d\in\{0,1,2,3\}$.We start by defining a mesh on each of them.While the different meshes in principle can be completely independent, the mathematical formulation is more easily expressed using conforming meshes: lower dimensional meshes are implicitly defined by their higher dimensional neighbour.

The sequence of meshes for each of the domains $\Omega_{d}$ will be denoted by $\mathcal{T}_{d}$, the set of edges of $\mathcal{T}_{d}$ by $\mathcal{E}_{d}$.Even though the interfaces $\Gamma_{d}$ are co-located with the lower dimensional domains $\Omega_{d}$, we keep the two separate by defining interface meshes $\mathcal{G}_{d}$.In what follows, $K$ will denote a generic element of mesh $\mathcal{T}_{d}$, $\gamma$ a generic element of the interface mesh $\mathcal{G}_{d}$ and $\sigma$ a generic face of $\mathcal{E}_{d}$.We will also indicate the normal pointing inside the fracture at element $\gamma$ as $\bm{n}_{\gamma}$.

The set $\mathcal{E}_{d}$, or equivalently the faces of any element $K$, can be partitioned into three sets:(i) internal or belonging to the Dirichlet boundary,(ii) belonging to the Neumann boundary,(iii) adjacent to a lower-dimensional domain.Notice that we get two different partitions based on which boundary conditions we use to split the boundary (boundary conditions related to the flow problem or to the transport problem):$\partial K=\partial K_{i}^{p}\cup\partial K_{N}^{p}\cup\partial K_{f}=\partialK$.Also, thanks to the conforming mesh hypothesis, faces in $\partial K_{f}$ can be identified with elements of the interface mesh $\mathcal{G}_{d}$, thus enabling us to legitimately write integrals such as $\int_{\gamma}\lambda_{d}$ where $\gamma\in\partial K_{f}$.