Free convection in fractured porous media: a numerical study
Arash Andrea Roknian Anna Scotti Alessio Fumagalli









Free convection in fractured porous media: a numerical study
    Arash Andrea Roknian^1 Anna Scotti^2 Alessio Fumagalli 3
    MOX - Dipartimento di Matematica “F. Brioschi”, Politecnico di
Milano, via Bonardi 9, 20133 Milan, Italy

    ^1 arashandrea.roknian@mail.polimi.it

    ^2 anna.scotti@polimi.it

    ^3 alessio.fumagalli@polimi.it
    September 9, 2024
===============================================================================================================================================================================================================================================



§ ABSTRACT 
The objective of this study is to better understand the influence of fractures
on the possibility of free convection in porous media. To this aim, we introduce a mathematical model for density driven flow in the presence of fractures, and the corresponding numerical approximation. In addition to the direct numerical solution of the problem we propose and implement a novel method for the assessment of convective stability through the eigenvalue analysis of the linearized numerical problem. The new method is shown to be in agreement with existing literature cases both in simple and complex fracture configurations. With respect to direct simulation in time, the results of the eigenvalue method lack information about the strength of convection and the steady state solution, they however provide detailed (quantitative) information about the behaviour of the solution near the initial equilibrium condition. Furthermore, not having to solve a time-dependent problem makes the method computationally very efficient. Finally, the question of how the porous matrix interacts with the fracture network to enable free convection is examined: the porous matrix is shown to be of key importance in enabling
convection for complex fracture networks, making stability criteria based on the fracture network alone somewhat limited in applicability.












§ 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 <cit.> 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 <cit.>.
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 <cit.>.

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 <cit.>, where the analysis is carried out analytically for homogenoues media, and in <cit.> 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 <ref> 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 <cit.>.

In section <ref> 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 <cit.> 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 <cit.>.

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

Section <ref> 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 <ref> 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 <ref> and <ref> test the two methods respectively with the Elder and HRL problems relying on the results of <cit.> and <cit.>: both are well known benchmark cases for density driven flow.
Section <ref> validates and examines a three-dimensional generalization of the HRL problem, based on the studies in <cit.>.
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 <ref>.

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



§ 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 <cit.>, 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.

Let us begin by considering a generic advection-diffusion equation expressing solute mass conservation in a porous medium:

    ∂_t (ρϕω) + ∇·(ρϕωu )+∇·(ρϕi )=0 ,


where ρ [M/L_v^3] is the density of the fluid in which the solute is dissolved, ω [M_s/M] is the concentration of the solute, u  [L/T] is the flow velocity, ρϕi  [M_s/L^2 T] is the solute mass flux due to diffusion and ϕ  [L_v^3/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.

For the diffusive flux a Fick-type law is used: i=  -D ∇ω .
The diffusion tensor D  [L^2/T] is usually modeled as the sum of two components: 
a part related to molecular diffusion and a part related to mechanical dispersion (function of the fluid velocity): D= D_d (u)+ D_m I.
In this study only the part related to molecular diffusion has been considered, thus 

    i=-D_m ∇ω .

In the following, it will be useful to indicate the combined advective and diffusive solute flux as
q [M_s/L^2 T]: 

    q = ωu + i .

The boundary conditions for equation (<ref>) can be 
of Dirichlet type where we set the value of the solute concentration ω = ω_D, 
or of Neumann type where the total mass flux in the normal direction is specified q·n = q_N, where n is the unit normal (by convention pointing outwards from the domain).

To close the system we still need an expression for the fluid velocity and a constitutive equation for density as a function of primary variables.
The first can be determined as the solution of the Darcy problem, which is the model used to describe filtration in porous media. In particular we have a mass conservation law for the fluid:

    ∂_t (ρϕ)+∇·(ρϕu )=0,

and a constitutive law classically used in the context of porous media, the Darcy law:

    u = k/ϕμ (- ∇ p+ρg) .


where g=-ge_z is the gravity acceleration vector, k is the permeability tensor of the porous medium k  [L_v^2], and μ [M/L_vT] is the viscosity of the fluid.

The boundary condition for (<ref>) can be 
of Dirichlet type where the pressure is specified: p = p_D 
or of Neumann type where the normal flow velocity is prescribed, u·n = u_N.

For the density constitutive law, we will assume a linear dependence on the concentration ω

    ρ = ρ(p, ω) = ρ_0 (1 + αω) ,

where ρ_0 is a reference value and α a dilation coefficient.

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.





 §.§ Oberbeck-Boussinesq approximation


The analysis and the solution of system (<ref>) -(<ref>)  can be substantially simplified by assuming the solution satisfies what is known as the Oberbeck-Boussinesq (OB) approximation.

The OB approximation consists in neglecting all but the most important among the different density related nonlinearities 
in the system.
In particular we will take ρ = ρ_0 (1 + αω) ≈ρ_0 everywhere except in the gravity term. 

Using this approximation, and after some algebraic simplifications the system reduces to 

    {   ∇·u = 0        in Ω, 
       ∂_t ω + ∇·q = 0         in Ω, 
       u= k/ϕμ (- ∇ p_e + ρ_0 αωg)       in Ω, 
       q= - D ∇ω + ωu      in Ω, 
    .

with boundary conditions

    {   u·n = u_N         on ∂Ω_N^p, 
        p = p_D       on ∂Ω_D^p,
    .    {   q·n = q_N       on ∂Ω_N^ω, 
       ω = ω_D       on ∂Ω_D^ω, 
    .

where ∂Ω = ∂Ω_N^p ∪∂Ω_D^p, ∂Ω = ∂Ω_N^ω∪∂Ω_D^ω.

Among the simplifications to reach the reduced system, we have replaced the unknown p with the excess pressure 
p_e by removing the hydrostatic component p_h= ρ_0 g (y_0 - y) (defined for the reference density, in the absence of solute):

    p  = p_e + p_h              = p_e + ρ_0 g (y_0 - y) , 
    
        -∇ p + ρg   = -∇ p_e - ρ_0 g +ρ_0 g + ρ_0 αωg
    

For readability, in what follows we will drop the subscript e and use the variable p to indicate the excess pressure p_e.





 §.§ 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 <cit.>), 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

    = v_g/v_d = k/ϕμρ_0 αω_max g /D/H ,


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 <cit.>, the possibility of convective motion is linked to the value of : _c = 4 π^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 <ref> 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 , defined as

    = ∫_A i·n/∫_A i_0 ·n =∫_A i·n/D/H A ,

where A is the diffusive inflow surface (top boundary) and i_0 is the diffusive flow in the absence of convection. When convection is not present inside the domain, i = i_0 and thus =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 i > i_0 and > 1.



 §.§ 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 


    k = k(x) = {    k_b       x∈Ω_b, 
        k_f               x∈Ω_f,
    
      
    .

where Ω_f denotes the fracture, Ω_b the bulk medium and Ω = Ω_b ∪Ω_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 <ref>. The domain Ω is split between bulk medium and fracture Ω = Ω_b ∪Ω_f ∪Γ, where Γ = Ω_b∩Ω_f denotes the interface between subdomains. We assume that Ω_f can be expressed as

    Ω_f = {x∈Ω:x=γ + αn, γ∈Γ^0, α∈ ( -b/2, b/2 ) } , 
       Γ^0 = {x∈Ω:x=x_0 + s_1(x_1 - x_0) + s_2(x_2 - x_0),  s_i ∈ (0,1) }  ,

where we have assumed a planar midsurface Γ^0 as shown in figure <ref>. The lateral boundary of Ω_f will be treated differently from the top and bottom boundaries. We define

    Γ^± = {x∈∂Ω_f : x = γ±b/2n, γ∈Γ^0 } , 
       Σ     = {x∈∂Ω_f : x = σ + αn,  σ∈∂Γ^0,  α∈ ( -b/2, b/2 ) }  ,

such that Γ = Γ^±∪Σ. 
Superscripts + and - will also be used to denote quantities evaluated on Γ^±.



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 Γ:


    {   u_b = 0       in Ω_b ,
       u_b=  (- ∇ p_b + ρ_0 αω_b g)           in Ω_b ,
       u_f = 0       in Ω_f ,
       u_f=  (- ∇ p_f + ρ_0 αω_f g)       in Ω_f ,
        p_b =  p_f       on Γ ,
       u_b ·n_b = u_f·n_b       on Γ ,
    .

where n_b is the normal vector defined on Γ, exiting from Ω_b.
At the end of dimensional reduction, the equations solved on the three-dimensional domain Ω_f will be replaced by the solution of (different) equations on the two-dimensional domain Γ^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 = ∫_-b/2^b/2u_f = ∫_-b/2^b/2 (Tu_f + Nu_f) 
    
      = ∇_τ·u_γ - (u_b ·n_b)^+ - (u_b ·n_b)^- ,

where we have denoted by T, N the projection operators specific to the fracture: 
Nv = (v·n) n and 
Tv = (I - N) v, where n is the normal vector to the fracture plane (as indicated in figure <ref>), and I is the identity operator.
We also denoted the in-plane gradient operator by ∇_τ = T∇, and by u_γ the integral of the tangential flow field:

    u_γ = ∫_-b/2^b/2Tu_f,    [L^2/T].

Indeed, given these definitions,

    ∫_-b/2^b/2 (Tu_f) = 
    ∫_-b/2^b/2∇_τ· (Tu_f) = 
    ∇_τ·∫_-b/2^b/2Tu_f =
    ∇_τ·u_γ.




To obtain a law for u_γ, we integrate the in-plane component of the Darcy law:

    u_γ 
    = ∫_-b/2^b/2Tu_f 
    = k_f/ϕμ∫_-b/2^b/2T (- ∇ p_f + ρ_0 αω_f g)
    = b k_f/ϕμ (-∇_τ p_γ + ρ_0 αω_γg_τ) ,

where p_γ = 1/b∫_-b/2^b/2 p_f, ω_γ = 1/b∫_-b/2^b/2ω_f and g_τ = Tg. Differently from u_γ, 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:

    (u_b ·n_b)^±≈b k_f/ϕμ (p_b^± - p_γ/b/2 + ρ_0 αω_b^± g_n).


We conclude by replacing Ω_f with its center plane Γ^0 and extending Ω_b:  Ω = Ω_b ∪Ω_f ∪Γ≈ (Ω∖Γ^0) ∪Γ^0 ∪Γ^±. Even though the fracture domain Ω_f and the interface Γ 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:

    {   u_b = 0 
          in Ω_b ,
       u_b=  (- ∇ p_b + ρ_0 αω_b g)
          in Ω_b ,
       ∇_τ·u_γ = [λ] 
          in Γ^0 ,
       u_γ = b (-∇_τ p_γ + ρ_0 αω_γg_τ) 
          in Γ^0 ,
       u_b ·n_b = λ      on Γ^±, 
       λ =  (- p_γ -  p_b/b/2 + ρ_0 αω_b g_n)           on Γ^±.
    .


In writing (<ref>), we have introduced the new variable λ [L/T] and the jump operator [v] = v^+ + v^-. 
For v defined on Γ^±, we will consider [v] to be defined on Γ^0.
Note that Σ collapsed onto a lower-dimensional object (the boundary of the center plane ∂Γ^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 λ 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 ∇ p. 



 §.§ 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

    {   ∂_t ω_b + ∇·q_b = 0 
          in Ω_b , 
       q_b= - D ∇ω_b + ω_b u_b
          in Ω_b ,  
       ∂_t ω_f + ∇·q_f = 0 
          in Ω_f ,  
       q_f= - D ∇ω_f + ω_f u_f            in Ω_f ,  
       ω_b = ω_f 
          on Γ ,  
       q_b ·n_b = q_f·n_b 
          on Γ . 
    .

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

    0
    = ∫_-b/2^b/2∂_t ω_f + ∇·q_f
    = b ∂_t ω_γ + ∇_τ·q_γ - (q_b ·n_b)^+ - (q_b ·n_b)^- = 0 .

As we did for the fluid mass conservation equation, we introduced fracture quantities ω_γ = 1/b∫_-b/2^b/2ω_f and the integrated tangential solute mass flux q_γ = ∫_-b/2^b/2Tq_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 ω_f into ω_γ + ω̃_f, where ω_γ is constant across the fracture and ω̃_f is a null average fluctuation, and similarly for u_f. We will neglect the product of fluctuations assuming that they are small, i.e. 

    ∫_-b/2^b/2 (ω_f u_f)
    ≈ω_γu_γ.


For fractures of non-negligible aperture however, internal motions (see e.g. intrafractrure mode 2A in figure <ref>) 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

    (q_b ·n_b)^± = D (ω_b^± - ω_γ/b/2) + ω_b^±λ^±


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


    {   ∂_t ω_b + ∇·q_b = 0                                                                            
          in Ω_b ,
       q_b = -D ∇ω_b + ω_b u_b                                                                           
          in Ω_b ,
       b ∂_t ω_γ + ∇_τ·q_γ = [θ]
          in Γ^0 ,
       q_γ = -b D ∇_τω_γ + ω_γu_γ      in Γ^0 ,
       q_b ·n_b = θ      on Γ^±, 
       θ = D (ω_b - ω_γ/b/2) + ω_b λ      on Γ^± ,
    .


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 (<ref>) and (<ref>), are entirely decoupled apart from their interaction through the interface variables λ and θ. 
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, possibly
intersecting fractures. In this case, the procedure is conceptually carried out hierarchically: n-dimensional
quantities are coupled to (n-1)-dimensional quantities through (n-1)-dimensional interface fluxes (figure <ref>).




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 ∈{0, 1, 2, 3} in a single domain Ω_d, and denote by Γ_d the interface between domains Ω_d and Ω_d+1. We also introduce mixed-dimensional variables ω_d, p_d, fluxes u_d, i_d and in-plane gradient operator ∇_d defined on Ω_d, interface fluxes λ_d, θ_d defined on Γ_d.

With these quantities at hand, we can generalize systems (<ref>, <ref>) to 

    {   ∇_d ·u_d = [λ_d] 
          in Ω_d ,
        b^3-d∂_t ω_d + ∇_d ·q_d = [θ_d]
          in Ω_d ,
       u_d = b^3-d k_d/ϕμ (-∇ p_d + ρ_0 αω_d g_d^τ) 
          in Ω_d ,
       q_d = - b^3-d D ∇ω_d + ω_d u_d 
          in Ω_d ,
       u_d ·n_d = λ_d-1      on Γ^d-1, 
       q_b ·n_b = θ_d-1      on Γ^d-1, 
       λ_d = b^2-d k_d/ϕμ ( p_d+1 - p_d /b/2 + ρ_0 αω_d g_d^n )
                  on Γ_d ,
       θ_d = b^2-d D ω_d+1 - ω_d /b/2 + ω_d λ_d 
          on Γ_d .
    .


with suitable boundary conditions.



 §.§ 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 _c, such that convective motion is possible for ≥_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 <cit.>, 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 <ref>, 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 <cit.>, 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. 





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

    ω_j = ω_maxy/H,  i_j = -D ∇ω_j = -D ω_max/He_y

for j ∈{ B, f }.

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

    ω_b = 
        {   (ω_max - δω)  y/H           y < H/2 , 
       (ω_max - δω)  y/H + δω         y > H/2 ,
        .
    
       ω_f = ω_max/2 
       δω = ω_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 ≪ H, 
thus making the error in the concentration profile small: δω≪ω_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 <ref> 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 <ref> on the right). 








§ NUMERICAL DISCRETIZATION



This section is dedicated to the discretization of system (<ref>).
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 Ω_d,  d ∈{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 Ω_d will be denoted by 𝒯_d, the set of edges of 𝒯_d by ℰ_d.
Even though the interfaces Γ_d are co-located with the lower dimensional domains Ω_d, we keep the two separate by defining interface meshes 𝒢_d.
In what follows, K will denote a generic element of mesh 𝒯_d, γ a generic element of the interface mesh 𝒢_d and σ a generic face of ℰ_d.
We will also indicate the normal pointing inside the fracture at element γ as n_γ.

The set ℰ_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): 
∂ K = ∂ K_i^p ∪∂ K_N^p ∪∂ K_f = ∂ K_i^ω∪∂ K_N^ω∪∂ K_f.
Also, thanks to the conforming mesh hypothesis, faces in ∂ K_f can be identified with elements of the interface mesh 𝒢_d, thus enabling us to legitimately write integrals such as ∫_γλ_d where γ∈∂ K_f.



With the notation in place, we can start integrating the conservation equations in system (<ref>) over a generic element K ∈𝒯_d:

    ∫_K ∇_d ·u_d  x^d = ∫_K [λ]  x^d , 
       ∫_K b^3-d∂ω_d/∂ t x^d + ∫_K ∇_d ·q_d  x^d = ∫_K [θ]  x^d .

Working on the integrals one by one we have

   ∫_K ∇_d ·u_d x^d 
= ∫_∂K u_d ·n x^d-1
= ∑_∂K_i ∫_σ u_d ·n x^d-1
+ ∑_∂K_f ∫_σ λ_d-1 x^d-1
+ ∑_∂K_N^p ∫_σ u_N x^d-1 , 

    ∫_K [λ_d] x^d = ∫_K^± λ_d x^d , 

    ∫_K b^3-d _t ω_d x^d = b^3-d d/dt ∫_K ω_d x^d , 

   ∫_K ∇_d ·q_d x^d 
= ∫_∂K q_d ·n x^d-1
= ∑_∂K_i ∫_σ q_d ·n x^d-1
+ ∑_∂K_f ∫_σ θ_d-1 x^d-1
+ ∑_∂K_N^ω ∫_σ q_N x^d-1 , 

    ∫_K [θ_d] x^d = ∫_K^± θ_d x^d . 



Now introducing the discrete variables and fluxes

    P_K = 1/|K|∫_K p_d  x^d ,
            U_Kσ≈∫_σu_d ·n_K  x^d-1 ,
           Λ_K^± = ∫_K^±λ_d  x^d , 
       W_K = 1/|K|∫_K ω_d  x^d ,
           Q_Kσ≈∫_σq_d ·n_K  x^d-1 ,
          Θ_K^± = ∫_K^±θ_d  x^d .

we can write the discrete version of the conservation equations in (<ref>):

    ∑_∂ K_i U_Kσ
    + ∑_∂ K_fΛ_γ
    + ∑_∂ K_N^p U_N,σ = Λ_K^± , 
        b^3-d |K| d/dt W_K
    + ∑_∂ K_i Q_Kσ 
    + ∑_∂ K_fΘ_γ 
    + ∑_∂ K_N^ω Q_N,σ 
    = Θ_K^± .

The last step is the discretization of the constitutive laws for the fluxes.
Integrating the constitutive laws gives

    ∫_σu_d ·n x^d-1 
    = b^3-d( 
        ∫_σ -∂ p_d/∂ n x^d-1 
        + ρ_0 α g·n_σ∫_σω_d  x^d-1) , 
       ∫_γλ_d   x^d 
    = b^2-d∫_γ(  p_d+1 - p_d/b/2 + ρ_0 αω_d+1 g·n_γ)  x^d , 
       ∫_σq_d ·n x^d-1 
    = b^3-d(
        D ∫_σ -∂ω_d/∂ n x^d-1
        + ∫_σω_d u_d·n x^d-1) , 
       ∫_γθ_d  x^d 
    = ∫_γ( b^2-d D ω_d+1 - ω_d/b/2  + ω_d λ_d )  x^d .

We rewrite each of these laws in terms of the discrete variables defined in (<ref>):

    U_Kσ =       b^3-d |σ| (   ∇ P_σ           + ρ_0 α g·n_σ  W_σ) , 
       Λ_K^± =   b^2-d |K|      (   P_K^± - P_K/b/2 + ρ_0 α g·n_K^±  W_K^±) , 
       Q_Kσ =       b^3-d |σ| D  ∇ W_σ           + W_σ U_Kσ , 
       Θ_K^± =    b^2-d |K|      D W_K^± - W_K/b/2 + W_K^±Λ_K^± ,


where the quantities ∇ϕ_σ and ϕ_σ, i.e. gradients and face values, depend on the particular finite volume scheme. We choose the Multipoint Flux Approximation scheme described in <cit.>, which computes the gradient on a face by considering values of all cells sharing a node with the face.
As for all finite volume schemes, fluid mass and solute mass conservation is guaranteed.
In contrast to a two-point scheme (TPFA) however, the MPFA scheme is consistent on general grids.

Note that we make use of a centered scheme for the advective term in the concentration equation. This choice, unlike upwind, is known to produce numerical oscillations for convection dominated flows. 
In all our numerical experiments however, convection is mild enough for our centered scheme to be numerically stable.
We can easily relate the already introduced Rayleigh number to the Peclet number:

    𝑃𝑒 = u h/D = u H/Dh/H = h/H

where h is a characteristic cell diameter, H is a characteristic length of the domain.

In all the the following numerical experiments, sufficiently fine grids will be computationally feasible for 𝑃𝑒 to be O(1).
The 𝑃𝑒 number is nonetheless numerically monitored in all the following simulations.



 §.§ Time discretization and direct solution method



We denote as "direct solution method" the integration of the model equations forward in time starting from a zero-solute or equilibrium initial condition, to assess the possible onset of convection.
Once a steady state has been reached, we must verify whether convective motion is present.
We will use the implicit Euler method for advancing in time for its unconditional stability. An adaptive time-stepping is also used to reduce the amount of computation necessary to reach steady state, see section <ref> for details.

To discretize in time, define a set of timesteps 
{t^n}_n=0… N
and write our discretized system as:

    ∑_∂ K_i U_Kσ^n+1
    + ∑_∂ K_fΛ_γ^n+1
    + ∑_∂ K_N^p U_N,σ^n+1 = Λ_K^±^n+1 , 
        b^3-d |K| W_K^n+1-W_K^n/Δ t^n
    + ∑_∂ K_i Q_Kσ^n+1 
    + ∑_∂ K_fΘ_γ^n+1
    + ∑_∂ K_N^ω Q_N,σ^n+1 
    = Θ_K^±^n+1,


where K∈𝒯_d,   d∈{0,1,2,3}, Δ t^n = t^n+1-t^n, n∈{0, …, N-1}.

If the timestep is long enough and the solution stops changing (according to the norm of the difference between two timesteps of the solution), we declare the system to have reached steady state.
Note that since the system of equations is nonlinear, each timestep requires the solution of a nonlinear problem. 
In our case we use Newton iterations by leveraging the automatic differentiation capabilities of the implementation. 
In the test cases, the tolerance for the Newton procedure is set to 1e-8 for the concentration increment.



 §.§ Eigenvalue analysis



The previously outlined method for assessing stability, while effective at predicting stability, has some shortcomings.
Among the disadvantages are its reliance on the choice of perturbation for the hydrostatic solution (if not starting from zero solute everywhere) and having to reach and assess the steadiness of the solution.
The computational cost of reaching the steady state may become an issue: advancing a nonlinear equation in time with an implicit scheme requires the solution of multiple linear systems for each timestep advancement. 

An alternative method, presented in detail below, relies on inspecting the nonlinear system of equations linearized around the equilibrium solution. 
The nonlinear discrete system (<ref>) and (<ref>) can be written abstractly as


    M d W/d t + F(W, Y)=0 , 
        G(W, Y)=0 ,

where  W ∈ℝ^n collects the degrees of freedom relative to the discrete variable W, and Y ∈ℝ^N-n collects the degrees of freedom relative to the discrete variables P, Λ and Θ. 
F: ℝ^n ×ℝ^N-n→ℝ^n and G: ℝ^n ×ℝ^N-n→ℝ^N-n collect the linear and nonlinear discrete operators in (<ref>).
Any discrete equilibrium solution (W_s, Y_s) will satisfy the system

    F(W_s, Y_s)=0 , 
       G(W_s, Y_s)=0 .

To assess  whether the equilibrium solution is also asymptotically stable we perturb the time-dependent system:


    M (d W_s/d t + d δ W/d t )+    F(W_s+ δ W, Y_s + δ Y)=0 , 
        G(W_s+ δ W, Y_s + δ Y)=0 ,

and linearize, taking advantage of the fact that δ W, δ Y are small perturbations:

    M d W_s/d t + M d δ W/d t +   0 F(W_s, Y_s) + ∂ F/∂ W|_W_s, Y_sδ W + ∂ F/∂ Y|_W_s, Y_sδ Y=0 , 
       0 G(W_s, Y_s) + ∂ G/∂ W|_W_s, Y_sδ W + ∂ G/∂ Y|_W_s, Y_sδ Y=0 .

Renaming the partial derivatives to ease the notation 

    A_ww = ∂F/∂W |_W_s, Y_s ,            A_wy = ∂F/∂Y |_W_s, Y_s , 

    A_yw = ∂G/∂W |_W_s, Y_s ,              A_yy = ∂G/∂Y |_W_s, Y_s , 


the system becomes

M ∂δW/∂t + 
   A_ww δW       + A_wy δY       = 0 , 
 
   A_yw δW       + A_yy δY       = 0 .

Now, relying on the invertibility of A_yy , we can eliminate δ Y to obtain a single  evolution equation for the perturbation δ W:

    M ∂δ W/∂ t = (A_wy A_yy^-1 A_yw - A_ww) δ W.

Substituting an exponential in time trial solution δ W(t) = w   e^λ t in the evolution equation yields an eigenvalue problem:

    λ M w = (A_wy A_yy^-1 A_yw - A_ww) w,

or, by defining the matrix S= M^-1 (A_wy A_yy^-1 A_yw - A_ww),

    S w = λ w .

Note that λ describes the evolution of the perturbation in time, while vector w its shape in space because each component represent the value in a grid cell. 
The equilibrium solution W_s is then linearly stable if and only if all the eigenvalues associated to system (<ref>) have negative real part. Conversely, if we can find one or more eigenvalues with positive real part the perturbation w will grow and be sustained in time.

Three computational considerations: (i) to assess stability there is no need to compute the entire spectrum, it is enough to compute the eigenvalue of largest real part (ii) automatic differentiation of the discrete system (<ref>) can yield the numerical value of the matrices A_ww, A_wy, A_yw, A_yy without having to write explicit expressions for them (iii) iterative methods for computation of eigenvalue spectra are available e.g. power iterations which do not require explicit expression for the matrix under study, only the ability of computing matrix-vector products. In our case, this removes the need of explicitly inverting the matrix A_yy. Furthermore, since the matrix only depends on the equilibrium solution, we can factorize the matrix only once using e.g. LU decomposition for fast matrix-vector products during the computation of the spectrum.



 §.§ Implementation details


The implementation of the numerical methods outlined above is based on the  library <cit.>, which provides the necessary tools for meshing fractured domains and assembling the discrete mixed-dimensional operators. The framework also implements forward automatic differentiation, providing the numerical Jacobians used for performing Newton iterations.

For the eigenvalue analysis (<ref>), the automatic differentiation part provides the Jacobian matrix.
Once the different matrix blocks are identified, we can easily define the matrix-vector product procedure yielding S v.
For the computation of few leading eigenpairs, the Krylov-Schur algorithm (outlined in <cit.>) has been combined with the dynamic restarting scheme described in <cit.>. 
Using the inner product defined by the mass matrix for the orthogonalization part of the algorithm has been particularly beneficial in accelerating convergence. 
The use of a custom procedure for computing eigenvalues has been preferred to a default implementation such as  mostly due to the possibility of monitoring convergence and as a possible starting point for devising more efficient algorithms.




§ RESULTS


The model and its numerical approximation have been validated against three reference papers: <cit.> which treats the Elder problem, <cit.> and <cit.> that treat the HRL scenario respectively in two and three dimensions. Both problems have been extensively used in the literature as benchmarks in the context of density driven flows.



 §.§ Elder problem



The Elder problem was originally proposed in a paper by Elder <cit.>, studying thermal convection in a Hele-Shaw cell.
It was later reformulated into a solute convection problem by Voss and Souza <cit.>,
where the system of equations is similar to the homogeneous problem (<ref>).

What this benchmark case aims to highlight and validate is the possibility of flow driven purely by density differences: no pressure gradient is being enforced by the boundary conditions.

The problem we want to solve is (<ref>), the domain being Ω = [0, 600] × [0,150]. Boundary conditions are

    {   u·n = 0         on ∂Ω ,
       ω = ω_max      on ∂Ω_i ,
       ω = 0       on ∂Ω_o ,
       q·n = 0         on ∂Ω∖ (∂Ω_i ∪∂Ω_o) ,
    .

where ∂Ω_i = (150, 450) × 150 and ∂Ω_o = (0, 600) × 0.
Since the boundary conditions for the Darcy problem are of Neumann type over the whole boundary we have an ill-posed problem; 
we can however restore the well-posedness by adding an additional constraint, such as imposing zero mean pressure over the whole domain:
∫_Ω p = 0.

Initial conditions prescribe ω(x, 0) = 0 ,  x ∈Ω.
Equations are integrated in time until T = 20yr.



All the other parameters of the problem are reported in table <ref>.




In the absence of gravity, and for small enough , the solute will diffuse from the inlet ∂Ω_i until the solution reaches the diffusive steady state. The characteristic time of the evolution is T_diff = H^2/D, which, for the parameters listed above is 
T = (150)^2/3.565e-6^2≈200yr.

The corresponding characteristic time associated to convection induced by density differences is much smaller: 

    T_adv = Hϕμ/k ρ_0 αω_max g≈0.5yr.

Note that the Rayleigh number presented above as a ratio of velocities <ref> can be equivalently interpreted as a ratio of these timescales:

    = T_diff/T_adv≈ 400 .


As noted in <cit.>, the solutions of the Elder problem present grid dependence even for moderately fine meshes, also due to the possible non-uniqueness of the solution.
For this reason, instead of trying to achieve grid independence, we directly validate the implementation by comparing the solutions with <cit.> for corresponding levels of grid refinement.

We use quadrilateral grids identical to the ones reported in the reference paper: number of cells n =  2^2l + 1 ,  l ∈{ 4, 5, 6 } and fixed timestep Δ t = [parse-numbers=false]1/12yr. The qualitative comparison using the contours of the concentration profiles is reported in figure <ref>.



While the differences in the continuous model (which is not explicitly detailed in <cit.>)
and discretization method (<cit.> uses adaptive time stepping and Galerkin-FEM) cause different solutions for low levels of grid refinement,
for higher levels of grid refinement the solutions are in good agreement.



 §.§ HRL problem



The HRL problem, already introduced in sections <ref> and <ref>, is the test case analyzed in <cit.>, 
which we will use as reference solutions for validating both the direct method and the eigenvalue method. 
The geometry of the domain and the parameters common to the different simulations are reported in figure <ref> and table <ref>.





The system of equations solved is system (<ref>), with boundary conditions given by

    {   u·n = 0         on ∂Ω ,
       ω = ω_max      on ∂Ω_i ,
       ω = 0       on ∂Ω_o ,
       q·n = 0         on ∂Ω∖ (∂Ω_i ∪∂Ω_o) ,
    .

where ∂Ω_i = (0, 20) × 10 and ∂Ω_o = (0, 20) × 0.
As for the Elder problem, the model is supplemented by the additional constraint ∫_Ω p = 0 to obtain a well-posed problem.
The solution strategy is inspired by the one outlined in the reference paper: the initial diffusive steady state concentration is perturbed
and the solution is advanced in time, gradually increasing the timestep Δ t.
Criteria for adapting the timestep Δ t include both the number of Newton iterations required for convergence and 
the norm of the concentration difference ‖ω^n+1 - ω^n‖.
If Δ t is large enough with respect to characteristic time scales, and the concentration difference between timesteps is small enough, we consider the solution to have reached steady state.

A quantitative comparison between our results and the one presented in the reference paper is presented in table <ref>.
Concentration profiles for a few selected cases are also presented in figure <ref>.
In all the test cases, there is agreement between our results and <cit.> on whether convective motion is present at steady state.
In the majority of cases, there is both qualitative agreement in the concentration profiles and quantitative agreement on the strength of convection, as measured by the Sherwood number, even in test cases with complex fracture configurations such as 𝖤9𝖺 and 𝖤9𝖻. 
As for the significant differences in cases such as 𝖡3, 𝖢2, 𝖢3
they might be due to important differences in the model solved in this study and the one in the reference paper:
we use the Boussinesq approximation, neglect dispersivity, and use a different numerical method. 
These differences however seem to have a modest overall impact.





The analysis of the different scenarios presented above can be complemented with the stability analysis based on eigenvalues outlined in section <ref>.
We begin with a detailed analysis of scenario 𝖣11.
The method can provide the eigenpairs corresponding to k eigenvalues with largest real part, for reasonably small k.
Figure <ref> shows some of the computed eigenpairs, and, as we can see, one of them is positive, indicating the presence of natural convection. 



The computed eigenvalues are consistent with the direct simulation: indeed, the instability of the diffusive steady state is confirmed by the presence of one eigenvalue of positive real part.
The numerical error due to the iterative nature of the eigenvalue computation is estimated by the formula ϵ = ‖ S x - λ x ‖/λ x, where S is introduced in (<ref>). Grid independence instead is assessed by computing the eigenvalues on a coarser grid and again computing a relative error: ϵ_g = | λ - λ_g |/λ, where λ_g are the eigenvalues computed on a coarser grid. 
The corresponding errors for this scenario are reported in table <ref>. The errors ϵ, ϵ_g are overall showing good accuracy for almost every scenario.



Apart from predicting the possibility of convection, the eigenfunctions and eigenvalues can give insight into the evolution dynamics in the vicinity of the diffusive initial condition.
We start by noticing that not only all of the computed eigenvalues have imaginary part equal to zero, but the related eigenfunctions furthermore are all mutually approximately orthogonal. Let ⟨ e_i, e_j⟩ denote the scalar product between (normalized) eigenfunctions: we have that, if i≠ j the value is about three orders of magnitude smaller than 1 for test case D11. We remark that the S matrix is not symmetric for reasons due to the implementation of boundary conditions, so there is no obvious reason to expect these results.

The orthogonality in particular suggests the possibility of studying the time evolution of the concentration in the space spanned by these eigenfunctions: given the solution ω(t) and the eigenfunctions { e_1, …, e_k }, let us define the scalar functions 

    α_i(t) = ⟨ω(t) - ω_0, e_i ⟩, i=1,…,k ,


representing the projections of ω(t)-ω_0 on the eigenfunction basis.
At steady state, different eigenfunctions give non-negligible contribution to the solution ω(t): no obvious relationship exists between the eigenpairs and the steady state solution. This should come as no surprise given that the eigenvalue analysis is localized around the initial condition.
In the early stages of the simulation however, not only the eigenvalue analysis correctly predicts the shape of the growing perturbation (identical to the only eigenfunction associated to a positive eigenvalue), but the eigenvalue λ_1 also provides a good estimate of its growth rate.



Let us now consider other scenarios.  Table <ref> provides the five eigenvalues with largest real part for each of the scenarios from <cit.> presented above. With no exception, a Sherwood number greater than 1 is associated to an eigenvalue with positive real part, thus correctly predicting the possibility or impossibility of convection. Moreover, larger  are (weakly) associated to greater number of eigenvalues or eigenvalues greater in magnitude.



As for the connection between the spectrum and the evolution dynamics at early times, scenarios 𝖠2, 𝖡2 and 𝖤9𝖻 have been analyzed, after checking the (approximate) orthogonality of the eigenfunctions.
In figure <ref>, the comparison between the actual time evolutions and the one predicted by the spectrum as ∑_i=0^N exp(λ_i t) e_i are represented. Here, we have used N=8 for cases 𝖠2, 𝖡2 and N=12 for 𝖤9𝖻.

The presence of multiple positive eigenvalues precludes the possibility of pinning down uniquely the shape and rate of growth of the instability. In all cases however, the growing perturbation has the form of one particular eigenfunction and its rate of growth is given by the associated eigenvalue, see figure <ref>.
The eigenvalue with largest real part thus also provides a numerical upper bound on the rate of growth.
 










 §.§ The role of fracture circuits



In <cit.> the onset and strength of convection is linked to the existence of continuous fracture circuits, and different numerical experiments confirm its importance.
In this section we argue that quasi-continuous fracture circuits can also enable convective motions, for a surrounding porous matrix of sufficiently high permeability.

Cases 𝖣1 and 𝖣2 in <cit.> are compared to show the necessity of circuit continuity for convection to occur.
However, slightly changing the fracture geometry (parameters remain unchanged), as in scenario 𝖣2^* (see figure <ref>), we see that, although the strength of convection is significantly reduced, it is still possible.
The permeability of the surrounding medium is indeed large enough for the medium to be part of the convective circuit, 
though not large enough for convection to occur without the presence of fractures.



Scenario 𝖣2^* is purposefully built to highlight this phenomenon, we can however find similar configurations in more complex scenarios already introduced in <cit.>.
Let us consider for instance network 𝖤9𝖻, shown in figure <ref>: velocity vectors plotted over the concentration profile clearly indicate, unlike what is claimed in <cit.>, one large convective cell instead of two, separate ones. 
As in scenario 𝖣2^*, the convective loop is not limited to the fracture network but crosses the porous matrix as well.



As further confirmation, we can decompose the solution of 𝖤9𝖻 using the approximate eigenvector basis (following the apprach outlined in section <ref>).
In figure <ref> we see how the dominant modes e_6, e_8, i.e. the modes corresponding to the largest magnitudes of α_i, both involve fluid motion across the gap in the fracture network, through the porous matrix,
while modes looping around large continuous fracture circuits e_1, e_2, e_4, e_7 contribute only marginally to the steady state solution.



Let us parametrize the geometry of a fracture circuit with a gap as in figure <ref>, to relate it to the possibility of convection with the analysis of eigenvalues sign. Results are reported in figure <ref>.
As expected, for very small matrix permeability gaps stop flow, thus inhibiting convection.
On the other hand, for larger matrix permeability k_m = 3e-16, the corresponding Rayleigh number is ≈ 20, very near the critical _c = 4 π^2 for which the matrix can exhibit convection even without the aid of fractures; convection is thus possible across very large gaps. 
For matrix permeability between these two extremes, the geometry of the gap determines whether convection across is possible or not.





This behavior can be explained, at least qualitatively, focusing on the top (broken) edge of the fracture circuit composed by the two fracture segments of length Δ x (segments A and C), and the gap of width ϵ (denoted by B). Assuming that across the gap flow occurs only across an area A, and that exchange with the porous matrix elsewhere re negligible, to ensure convection in the circuit we have to guarantee that the conductivity G_B≥ G_A=G_C, which translates into


    A k_mϵ≥ k_f bΔ x,


or, equivalently, 


    k_m Ak_f bΔ xϵ≥ 1.


This relation correctly predicts the qualitative behaviour of the numerical simulations reported in figure <ref>: small gaps with large surface areas favour convection across the gap. 
The results of a more quantitative evaluation, checking the sign of eigenvalues for different combinations of the parameters, are illustrated in figure <ref>. Although the threshold does not appear to be so sharp, and its numerical value is closer to 10^-1 than to 1, the horizontal separation of positive leading eigenvalues from negative ones indicates that the model is capturing part of the physics of the phenomenon.






 §.§ Three-dimensional HRL problem



Both the direct method (exposed in section <ref>) and the eigenvalue method (section <ref>) make no assumptions on the dimensionality of the ambient space.
In this section, we are going to apply the latter to simple three-dimensional scenarios.

In <cit.>, a classification of possible modes of convection in fractured porous media is proposed.
The paper further explains how the assumption of a two-dimensional domain severely limits the applicability to experimental cases.
This is due to the three-dimensional character of the dominant mode of convection – mode 2B in figure <ref> – which cannot be captured by a two-dimensional analysis.

In <cit.>, different three-dimensional scenarios are analyzed and numerically simulated. 
The numerical results from the regular three-dimensional fracture circuit (here denoted as scenario 6), show that when convection is available (for large enough fracture aperture), the dominant mode is the interfracture mode 1.

The geometry simulated in <cit.> being practically identical to the geometry studied in <cit.>, the two results stand in contradiction.
Indeed, <cit.> acknowledges the contradiction, 
attributing it to both the matrix-fracture coupling conditions and to the Rayleigh averaging strategy used in the analysis of <cit.>.
In particular, the use of a Rayleigh number based on an averaged permeability was shown to be ineffective in the case of low density fracture networks in <cit.> (as briefly discussed in section <ref>).



We analyze the 3D problem with the eigenvalue method for different scenarios presented in <cit.>, and collect the results in table <ref>, where we report for different geometries and different fracture apertures the unstable modes, and compare them with the convection modes predicted by <cit.>. 

As in the two-dimensional case, the instability thresholds are identical to the ones obtained in <cit.>, apart from a slight misprediction for scenario 9𝖺.

In all scenarios furthermore, the eigenfunctions clearly follow the classification proposed in <cit.>, corresponding to either interfracture □ or intrafracture ↻ modes.



In particular, for scenario 6 we observe that (1) the interfracture mode is first unstable one appearing at b ≈1.4e-5
(2) already at b ≈1.6e-5, different intrafracture modes are available and dominate over the interfracture mode:
λ^↻_1 = 21.0, λ^↻_2 = 20.1, λ^□ = 17.6
(3) intrafracture convection modes remain dominant for larger fracture apertures.

Thus, excluding apertures for which the setup is very near to the instability threshold, intrafracture convection modes dominate over interfracture convection modes when both modes are available.
This conclusion also holds for scenarios 9𝖺 and 9𝖻.

Although our results seem to confirm the results of <cit.>, where the intrafracture mode is also identified as dominant,
they do not contradict the results in <cit.>. 
Indeed the eigenvalue analysis is strictly localized around the diffusive equilibrium solution. In particular, the eigenvalues give no indication on which modes will be present at steady state. 
As seen in e.g. scenario 𝖠2 (figure <ref>), even modes with decaying behaviour near equilibrium may turn out to be dominant at steady state.






 §.§ Computational considerations



We want to discuss some of the computational aspects of the eigenvalue method, comparing it in particular to the direct method of assessing stability in terms of computational cost.

The computationally expensive steps of the eigenvalue method can be divided in 
(a) grid construction, initial operator discretization and Newton iterations to reach steady state (needed as initial condition)
(b) construction of the S operator defined in (<ref>) (including LU factorization of A_yy)
(c) computation of k largest eigenvalues.
Step (c) can itself be subdivided in the two costly operations (i) matrix vector products S v and (ii) orthogonalization.
We will initially compare the computational cost of assessing the possibility of convection using the eigenvalue method against the direct method.
Next, we will discuss how the cost changes with number of degrees of freedom and number of sought eigenvalues.

In section <ref>, we saw how the direct method and the eigenvalue method can be used in a complementary fashion,
each of them giving answers inaccessible to the other.
We can however restrict the question to the possibility of convection since both methods are a valid way of finding an answer (e.g. see table <ref>).
For a fairer comparison, the direct forward simulations will be as soon as the Sherwood number exceeds 1, which would already indicate the possibility of convection.
Table <ref> reports the timing comparison and the relative time spent in each of the most costly operations for different test cases.

The cost of the direct method 
is dominated by the assembly of the problem Jacobian in the two-dimensional cases
and by the linear solves in the three-dimensional cases (which possess a much larger number of degrees of freedom).
Both of these operations are performed for each Newton iteration. 
The number of Newton iterations per timestep advancement mostly varies between 2 and 3 for all the simulations.
Relaxing the convergence tolerance of the Newton method is one of the options to speed up the method, at the cost of lower accuracy of the solution.
This might be a valid option considering that what we are interested in is whether convection is possible, not in the detailed behaviour of the solution.

The eigenvalue method is much more efficient with respect to the direct method in low fracture density cases.
For test cases with complex fracture configurations such as 𝖤9𝖺 and 𝖤9𝖻, the performance of the two methods (direct and eigenvalue) approximately matches.

In the different test cases, the eigenvalue method dedicates approximately equal times to orthogonalization and matrix-vector products.
Their relative weights can actually be adjusted by varying the size m of the Krylov basis in the Krylov-Schur algorithm:
for every iteration of the algorithm, the cost of orthogonalization scales with m^2 while the cost of the matrix-vector products scales with m.
Tuning m for the problem at hand may be a possible option for improving the efficiency of the method.



In figure <ref> we illustrate how the number of matrix-vector products N_mv changes with the number of mesh elements .
The results show that N_mv depends weakly (N_mv∼^0.5) on the size of the problem (although a lot of variability remains).
Note that however the cost of each matrix-vector product S v scales as C_mv∼^2: the cost of each matrix-vector product S v scales as C_mv∼ N_mv^2: the presence of A_yy^-1 in (<ref>) makes the S matrix dense.
The scaling estimate for the total cost of the matrix-vector products 𝑇𝐶_mv is thus

    𝑇𝐶_mv = N_mv C_mv∼^0.5^2 = ^2.5 .

The estimate is consistent with the data reported in table <ref>, for which a least-squares estimate gives 𝑇𝐶_mv∼^2.2 ± 0.4.



Finally, the convergence of the  Krylov-Schur algorithms is represented for four different test cases during the earch for multiple eigenvalues. As shown in figure <ref> in all the test cases, once the algorithms pins down an eigenvalue, the corresponding error starts decreasing at exponential rate. Depending on the test case, we can have scenarios in which the eigenvalues all practically converge together, as opposed to cases where computing each successive eigenvalue requires considerably more computation.






§ CONCLUSION



The aim of this study was to better understand the influence of fractures on the possibility of free convection in porous media.
To this aim, we have described a mathematical model for density driven flow in the presence of fractures,
and the corresponding numerical approximation. In addition to the direct ”forward” numerical solution of the problem
we have proposed and implemented a novel method for the assessment of convective stability through the eigenvalue analysis of the linearized numerical problem.

The new method is shown to be in agreement with existing literature cases both in simple and complex fracture configurations.
With respect to direct simulation in time, its results provide further information on the possibility of convection.
In particular, as shown in section <ref> the sign of the leading eigenvalue correctly predicts the onset of convection, and its magnitude provides quantitative estimates of the rate of growth or decay of perturbations.
Furthermore the computational cost of the method has proven to be in the worst cases equal, and in the best cases up to an order of magnitude faster than the direct solution method.

The fact that the eigenvalue method closely mimics the analytical method of investigating stability clarifies what inferences can be made from the the results of a linear stability analysis. In particular, analyzing stability through the study of the system around equilibrium solutions (both numerically and analytically) is restrictive in that it cannot predict whether, far from the equilibrium solution, transient convective motion is possible and which convective modes will be dominant at steady state.

Moreover, an in depth study of the particular scenario E9b in <cit.> has further complicated the question of analyzing free convection in the presence of fractures: for realistic sets of problem parameters, the porous matrix is indeed able to participate in the convective motion. 
Thus, stability criteria based on the fracture network alone, e.g. the presence of large continuous fracture circuits as a trigger for convection, are shown to be somewhat limited in applicability.

Given the results of this study, we could expand the work in different directions.

On the numeric side, the method used to compute the eigenvalues at the moment does not take into account the structure of the particular problem.
However, the S matrix is built up from submatrices of the Jacobian associated to the linearized problem.
Therefore, a close study of the Jacobian may suggest ways of speeding up the computation of the eigenvalues by exploiting the structure of the S matrix.

Finally, though the eigenvalue method has proven to be faster than the direct method of assessing convective stability, the method could be further sped up if we are willing to sacrifice accuracy, since, for this purpose, we are only interested in the sign of the most positive eigenvalue. If during the computation of the eigenvalues we find one to be positive, with error reasonably smaller than its magnitude, we can already predict the possibility of convective motion.
In cases where the method is applicable, such as the fractured HRL problem, the additional speed 
may enable a more systematic (or even statistical) study of how global properties of fracture networks such as fracture density, connectivity and characteristic sizes are related to the possibility of convective motion.



§ ACKNOWLEDGEMENTS

The authors acknowledge the support by MUR, grant Dipartimento di Eccellenza 2023–2027.


10

aavatsmark02
Ivar Aavatsmark.
An Introduction to Multipoint Flux Approximations for
  Quadrilateral Grids.
Computational Geosciences, 6(3/4):405–432, 2002.

diersch02
H.-J.G. Diersch and O. Kolditz.
Variable-density flow and transport in porous media: approaches and
  challenges.
Advances in Water Resources, 25(8-12):899–944, August 2002.

elder67
J. W. Elder.
Steady free convection in a porous medium heated from below.
Journal of Fluid Mechanics, 27(1):29–48, January 1967.

hortonrogers
C. W. Horton and F. T. Rogers.
Convection Currents in a Porous Medium.
Journal of Applied Physics, 16(6):367–370, June 1945.

porepy
Eirik Keilegavlen, Runar Berge, Alessio Fumagalli, Michele Starnoni, Ivar
  Stefansson, Jhabriel Varela, and Inga Berre.
PorePy: an open-source software for simulation of multiphysics
  processes in fractured porous media.
Computational Geosciences, 25(1):243–265, February 2021.

lapwood_1948
E. R. Lapwood.
Convection of a fluid in a porous medium.
Mathematical Proceedings of the Cambridge Philosophical
  Society, 44(4):508–521, October 1948.

martin05
Vincent Martin, Jérôme Jaffré, and Jean E. Roberts.
Modeling Fractures and Barriers as Interfaces for Flow in
  Porous Media.
SIAM Journal on Scientific Computing, 26(5):1667–1691, January
  2005.

rees
D.A.S. Rees and L. Storesletten.
The onset of convection in a two-layered porous medium with
  anisotropic permeability.
Transport in Porous Media, 128, 2019.

shafabakhsh19
Paiman Shafabakhsh, Marwan Fahs, Behzad Ataie-Ashtiani, and Craig T. Simmons.
Unstable Density-Driven Flow in Fractured Porous Media:
  The Fractured Elder Problem.
Fluids, 4(3):168, September 2019.

shikaze98
Steven G Shikaze, E.A Sudicky, and F.W Schwartz.
Density-dependent solute transport in discretely-fractured geologic
  media: is prediction possible?
Journal of Contaminant Hydrology, 34(3):273–291, October 1998.

nield
C. T. Simmons, J. M. Sharp, and D. A. Nield.
Modes of free convection in fractured low‐permeability media.
Water Resources Research, 44(3):2007WR006551, March 2008.

starnoni19
M. Starnoni, I. Berre, E. Keilegavlen, and J. M. Nordbotten.
Consistent MPFA Discretization for Flow in the Presence of
  Gravity.
Water Resources Research, 55(12):10105–10118, December 2019.

stathopoulos98
Andreas Stathopoulos, Yousef Saad, and Kesheng Wu.
Dynamic Thick Restarting of the Davidson, and the Implicitly
  Restarted Arnoldi Methods.
SIAM Journal on Scientific Computing, 19(1):227–245, January
  1998.

stewart02
G. W. Stewart.
A Krylov–Schur Algorithm for Large Eigenproblems.
SIAM Journal on Matrix Analysis and Applications,
  23(3):601–614, January 2002.

voss1987variable
Clifford I. Voss and William R. Souza.
Variable density flow and solute transport simulation of regional
  aquifers containing a narrow freshwater‐saltwater transition zone.
Water Resources Research, 23(10):1851–1866, October 1987.

vg15
Katharina Vujević and Thomas Graf.
Combined inter- and intra-fracture free convection in fracture
  networks embedded in a low-permeability matrix.
Advances in Water Resources, 84:52–63, October 2015.

vg14
Katharina Vujević, Thomas Graf, Craig T. Simmons, and Adrian D. Werner.
Impact of fracture network geometry on free convective flow patterns.
Advances in Water Resources, 71:65–80, September 2014.




