Arash Andrea Roknian1Anna Scotti2Alessio Fumagalli
(MOX - Dipartimento di Matematica “F. Brioschi”, Politecnico diMilano, via Bonardi 9, 20133 Milan, Italy 1 arashandrea.roknian@mail.polimi.it 2 anna.scotti@polimi.it 3 alessio.fumagalli@polimi.itSeptember 5, 2024)
Abstract
The objective of this study is to better understand the influence of fractureson 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 enablingconvection for complex fracture networks, making stability criteria based on the fracture network alone somewhat limited in applicability.
Keywords: free convection, variable density flow, fractured porous media, linear stability analysis
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.
Let us begin by considering a generic advection-diffusion equation expressing solute mass conservation in a porous medium:
(1)
where is the density of the fluid in which the solute is dissolved, is the concentration of the solute, is the flow velocity, is the solute mass flux due to diffusion and is the porosity of the medium.The subscript in the length dimension is used to differentiate void volumes from total volumes, while is used to differentiate solute mass from fluid mass.
For the diffusive flux a Fick-type law is used: The diffusion tensor 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): .In this study only the part related to molecular diffusion has been considered, thus
(2)
In the following, it will be useful to indicate the combined advective and diffusive solute flux as:
(3)
The boundary conditions for equation (1) can beof Dirichlet type where we set the value of the solute concentration ,or of Neumann type where the total mass flux in the normal direction is specified , where 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:
(4)
and a constitutive law classically used in the context of porous media, the Darcy law:
(5)
where is the gravity acceleration vector, is the permeability tensor of the porous medium , and is the viscosity of the fluid.
The boundary condition for (4) can beof Dirichlet type where the pressure is specified: or of Neumann type where the normal flow velocity is prescribed, .
For the density constitutive law, we will assume a linear dependence on the concentration
(6)
where 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.
2.1 Oberbeck-Boussinesq approximation
The analysis and the solution of system (1) -(6) 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 nonlinearitiesin the system.In particular we will take everywhere except in the gravity term.
Using this approximation, and after some algebraic simplifications the system reduces to
(7)
with boundary conditions
where , .
Among the simplifications to reach the reduced system, we have replaced the unknown with the excess pressure by removing the hydrostatic component (defined for the reference density, in the absence of solute):
For readability, in what follows we will drop the subscript and use the variable to indicate the excess pressure .
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 and , is described by the Rayleigh number
(8)
where 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 : is the critical number i.e. a threshold such that, for 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 , defined as
(9)
where is the diffusive inflow surface (top boundary) and is the diffusive flow in the absence of convection. When convection is not present inside the domain, and thus . Note that the boundary conditions for the HRL problem being of zero fluid mass flow all around, the solute transport on the inflow surface 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 makes and .
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
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 is split between bulk medium and fracture , where denotes the interface between subdomains. We assume that can be expressed as
where we have assumed a planar midsurface as shown in figure 3. The lateral boundary of will be treated differently from the top and bottom boundaries. We define
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 :
(10)
where is the normal vector defined on , exiting from .At the end of dimensional reduction, the equations solved on the three-dimensional domain will be replaced by the solution of (different) equations on the two-dimensional domain .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.
where we have denoted by , the projection operators specific to the fracture: and, where is the normal vector to the fracture plane (as indicated in figure 3), and is the identity operator.We also denoted the in-plane gradient operator by , and by the integral of the tangential flow field:
Indeed, given these definitions,
To obtain a law for , we integrate the in-plane component of the Darcy law:
where , and . Differently from , 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:
We conclude by replacing with its center plane and extending : . Even though the fracture domain 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:
(11)
In writing (11), we have introduced the new variable and the jump operator .For defined on , we will consider to be defined on .Note that collapsed onto a lower-dimensional object (the boundary of the center plane ) 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 ) 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 .
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
(12)
We integrate the conservation equation over the fracture, splitting the fluxes into their normal and tangential parts:
As we did for the fluid mass conservation equation, we introduced fracture quantities and the integrated tangential solute mass flux . 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 into , where is constant across the fracture and is a null average fluctuation, and similarly for . We will neglect the product of fluctuations assuming that they are small, i.e.
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
Finally, we obtain a mixed-dimensional system for the transport of a solute in a fractured porous medium,
(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 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 -dimensional domain () with multiple, possiblyintersecting fractures. In this case, the procedure is conceptually carried out hierarchically: -dimensionalquantities are coupled to -dimensional quantities through -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 in a single domain , and denote by the interface between domains and . We also introduce mixed-dimensional variables , fluxes and in-plane gradient operator defined on , interface fluxes defined on .
With these quantities at hand, we can generalize systems (11, 13) to
(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 , such that convective motion is possible for .
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.
for .
With dimensional reduction, the system of equations to be solved is replaced by (13), which yields a piecewise linear concentration profile:
being the fracture aperture. Fractures which qualify as thin enough to be treated as lower-dimensional regions will always satisfy ,thus making the error in the concentration profile small: . 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 .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 will be denoted by , the set of edges of by .Even though the interfaces are co-located with the lower dimensional domains , we keep the two separate by defining interface meshes .In what follows, will denote a generic element of mesh , a generic element of the interface mesh and a generic face of .We will also indicate the normal pointing inside the fracture at element as .
The set , or equivalently the faces of any element , 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):.Also, thanks to the conforming mesh hypothesis, faces in can be identified with elements of the interface mesh , thus enabling us to legitimately write integrals such as where .
With the notation in place, we can start integrating the conservation equations in system (14) over a generic element :
Working on the integrals one by one we have
Now introducing the discrete variables and fluxes
(15)
we can write the discrete version of the conservation equations in (14):
(16)
The last step is the discretization of the constitutive laws for the fluxes.Integrating the constitutive laws gives
We rewrite each of these laws in terms of the discrete variables defined in (15):
(17)
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 [1], 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:
where is a characteristic cell diameter, is a characteristic length of the domain.
In all the the following numerical experiments, sufficiently fine grids will be computationally feasible for to be .The number is nonetheless numerically monitored in all the following simulations.
3.1 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 4.2 for details.
To discretize in time, define a set of timestepsand write our discretized system as:
(18)
where , , .
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 for the concentration increment.
3.2 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 (16) and (17) can be written abstractly as
(19)
where collects the degrees of freedom relative to the discrete variable , and collects the degrees of freedom relative to the discrete variables , and . and collect the linear and nonlinear discrete operators in (16).Any discrete equilibrium solution will satisfy the system
To assess whether the equilibrium solution is also asymptotically stable we perturb the time-dependent system:
and linearize, taking advantage of the fact that are small perturbations:
Renaming the partial derivatives to ease the notation
the system becomes
Now, relying on the invertibility of , we can eliminate to obtain a single evolution equation for the perturbation :
Substituting an exponential in time trial solution in the evolution equation yields an eigenvalue problem:
or, by defining the matrix ,
(20)
Note that describes the evolution of the perturbation in time, while vector its shape in space because each component represent the value in a grid cell.The equilibrium solution is then linearly stable if and only if all the eigenvalues associated to system (20) have negative real part. Conversely, if we can find one or more eigenvalues with positive real part the perturbation 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 (19) can yield the numerical value of the matrices 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 . 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.
3.3 Implementation details
The implementation of the numerical methods outlined above is based on the PorePy library [5], 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 (3.2), 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 .For the computation of few leading eigenpairs, the Krylov-Schur algorithm (outlined in [14]) has been combined with the dynamic restarting scheme described in [13].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 ARPACK mostly due to the possibility of monitoring convergence and as a possible starting point for devising more efficient algorithms.
4 Results
The model and its numerical approximation have been validated against three reference papers: [2] which treats the Elder problem, [17] and [16] 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.
4.1 Elder problem
The Elder problem was originally proposed in a paper by Elder [3], studying thermal convection in a Hele-Shaw cell.It was later reformulated into a solute convection problem by Voss and Souza [15],where the system of equations is similar to the homogeneous problem (7).
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 (7), the domain being . Boundary conditions are
where and .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:.
Initial conditions prescribe .Equations are integrated in time until .
All the other parameters of the problem are reported in table 1.
Permeability
Porosity
0.1
1
Viscosity
Freshwater density
Solute expansion coefficient
1
Maximum concentration
1
1
Gravitational acceleration
9.81
Diffusivity
In the absence of gravity, and for small enough , the solute will diffuse from the inlet until the solution reaches the diffusive steady state. The characteristic time of the evolution is , which, for the parameters listed above is.
The corresponding characteristic time associated to convection induced by density differences is much smaller:
Note that the Rayleigh number presented above as a ratio of velocities 8 can be equivalently interpreted as a ratio of these timescales:
As noted in [2], 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 [2] for corresponding levels of grid refinement.
We use quadrilateral grids identical to the ones reported in the reference paper: number of cells and fixed timestep . The qualitative comparison using the contours of the concentration profiles is reported in figure 10.
While the differences in the continuous model (which is not explicitly detailed in [2])and discretization method ([2] 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.
4.2 HRL problem
The HRL problem, already introduced in sections 2.2 and 2.5, is the test case analyzed in [17],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 11 and table 2.
Permeability
Porosity
0.1
1
Viscosity
Freshwater density
Solute expansion coefficient
1
Maximum concentration
0.1
1
Gravitational acceleration
9.81
Diffusivity
The system of equations solved is system (14), with boundary conditions given by
where and .As for the Elder problem, the model is supplemented by the additional constraint 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 perturbedand the solution is advanced in time, gradually increasing the timestep .Criteria for adapting the timestep include both the number of Newton iterations required for convergence andthe norm of the concentration difference .If 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 3.Concentration profiles for a few selected cases are also presented in figure 12.In all the test cases, there is agreement between our results and [17] 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 and .As for the significant differences in cases such as 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.
Scenario
1.00
1.00
0.00
1.36
1.38
0.01
1.56
1.77
0.14
1.75
1.79
0.02
1.00
1.00
0.00
1.13
1.18
0.04
1.49
1.19
-0.20
1.32
1.29
-0.03
1.00
1.00
0.00
1.17
1.02
-0.13
1.21
1.06
-0.13
1.21
1.13
-0.06
1.06
1.08
0.02
Scenario
1.08
1.08
0.00
1.00
1.00
0.00
1.00
1.00
0.00
1.01
1.02
0.01
1.07
1.06
-0.01
1.08
1.09
0.01
1.28
1.27
-0.01
1.15
1.16
0.01
1.15
1.16
0.01
1.45
1.45
0.00
1.37
1.37
0.00
1.39
1.39
0.00
2.69
2.92
0.09
The analysis of the different scenarios presented above can be complemented with the stability analysis based on eigenvalues outlined in section 3.2.We begin with a detailed analysis of scenario .The method can provide the eigenpairs corresponding to eigenvalues with largest real part, for reasonably small .Figure 13 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 , where is introduced in (20). Grid independence instead is assessed by computing the eigenvalues on a coarser grid and again computing a relative error: , where are the eigenvalues computed on a coarser grid.The corresponding errors for this scenario are reported in table 4. The errors , are overall showing good accuracy for almost every scenario.
61.75
0.0000
0.0014
-9.86
0.0001
0.0004
-15.43
0.0001
0.0007
-18.83
0.0001
0.0007
-37.26
0.0001
0.0014
-39.40
0.0004
0.0015
-41.48
0.0003
0.0016
-47.97
0.0008
0.0020
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 denote the scalar product between (normalized) eigenfunctions: we have that, if the value is about three orders of magnitude smaller than 1 for test case D11. We remark that the 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 and the eigenfunctions , let us define the scalar functions
(21)
representing the projections of on the eigenfunction basis.At steady state, different eigenfunctions give non-negligible contribution to the solution : 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 also provides a good estimate of its growth rate.
Let us now consider other scenarios. Table 5 provides the five eigenvalues with largest real part for each of the scenarios from [17] 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.
Scenario
1.00
-9.87
-10.88
-16.61
-25.55
-39.35
0.00919
0.00075
1.84
1.36
51.92
30.36
2.30
-9.87
-39.43
0.01402
0.01351
1.86
1.56
16.03
13.45
8.44
0.17
-6.44
0.03232
0.07508
1.93
1.75
11.08
10.99
0.69
-0.17
-9.87
0.00652
0.03490
1.47
1.00
-9.87
-10.76
-15.70
-24.87
-39.43
0.00500
0.00110
1.75
1.13
31.29
26.31
0.96
-9.87
-18.21
0.00265
0.04576
1.73
1.49
15.31
9.31
1.63
-4.05
-9.87
0.00658
0.04265
1.64
1.32
11.62
4.39
-0.54
-3.63
-9.87
0.02193
0.05261
1.39
1.00
-9.87
-10.59
-13.33
-26.06
-40.57
0.00827
0.02932
1.84
1.17
1.90
0.91
-4.11
-9.87
-39.45
0.00155
0.32562
1.59
1.21
1.02
-1.06
-8.18
-9.87
-23.91
0.01008
0.08807
1.53
1.21
1.84
-1.02
-6.77
-9.87
-22.38
0.00230
0.03980
1.61
1.08
41.75
-9.86
-14.06
-16.43
-39.39
0.00110
0.00526
1.71
1.00
-8.52
-9.86
-16.43
-20.71
-39.40
0.00012
0.00967
1.71
1.00
-1.99
-9.86
-16.44
-16.68
-39.39
0.00093
0.02219
1.71
1.02
17.12
-9.86
-13.38
-16.40
-39.39
0.00015
0.02364
1.68
1.06
35.55
-9.86
-13.60
-16.43
-38.18
0.00003
0.00642
1.69
1.09
48.39
-9.86
-13.02
-15.31
-39.40
0.00097
0.01294
1.71
1.27
305.11
-9.86
-11.95
-15.69
-37.38
0.00089
0.02138
1.74
1.16
47.51
47.21
-9.86
-11.93
-27.74
0.00002
0.00433
1.68
1.16
84.85
7.23
-9.86
-26.09
-26.26
0.00418
0.01027
1.68
1.45
160.60
-9.86
-13.42
-16.15
-17.35
0.00001
0.01819
1.70
1.37
61.75
-9.86
-15.43
-18.83
-37.26
0.00014
0.00139
1.77
1.39
73.01
22.01
-9.86
-12.47
-20.32
0.00009
0.01913
1.74
1.06
239.23
67.69
32.76
-9.87
-10.11
0.00707
0.05096
1.60
2.92
495.15
290.82
238.04
178.01
86.36
0.00016
0.01846
1.53
As for the connection between the spectrum and the evolution dynamics at early times, scenarios , and have been analyzed, after checking the (approximate) orthogonality of the eigenfunctions.In figure 15, the comparison between the actual time evolutions and the one predicted by the spectrum as are represented. Here, we have used for cases , and for .
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 16.The eigenvalue with largest real part thus also provides a numerical upper bound on the rate of growth.
4.3 The role of fracture circuits
In [17] 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 and in [17] 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 (see figure 17), 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 is purposefully built to highlight this phenomenon, we can however find similar configurations in more complex scenarios already introduced in [17].Let us consider for instance network , shown in figure 18: velocity vectors plotted over the concentration profile clearly indicate, unlike what is claimed in [17], one large convective cell instead of two, separate ones.As in scenario , 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 using the approximate eigenvector basis (following the apprach outlined in section 4.2).In figure 19 we see how the dominant modes , i.e. the modes corresponding to the largest magnitudes of , both involve fluid motion across the gap in the fracture network, through the porous matrix,while modes looping around large continuous fracture circuits contribute only marginally to the steady state solution.
Let us parametrize the geometry of a fracture circuit with a gap as in figure 20, to relate it to the possibility of convection with the analysis of eigenvalues sign. Results are reported in figure 21.As expected, for very small matrix permeability gaps stop flow, thus inhibiting convection.On the other hand, for larger matrix permeability , the corresponding Rayleigh number is , very near the critical 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 (segments and ), and the gap of width (denoted by ). Assuming that across the gap flow occurs only across an area , and that exchange with the porous matrix elsewhere re negligible, to ensure convection in the circuit we have to guarantee that the conductivity , which translates into
or, equivalently,
This relation correctly predicts the qualitative behaviour of the numerical simulations reported in figure 21: 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 22. Although the threshold does not appear to be so sharp, and its numerical value is closer to than to , the horizontal separation of positive leading eigenvalues from negative ones indicates that the model is capturing part of the physics of the phenomenon.
4.4 Three-dimensional HRL problem
Both the direct method (exposed in section 3.1) and the eigenvalue method (section 3.2) 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 [11], 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 23 – which cannot be captured by a two-dimensional analysis.
In [16], 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 [16] being practically identical to the geometry studied in [11], the two results stand in contradiction.Indeed, [16] acknowledges the contradiction,attributing it to both the matrix-fracture coupling conditions and to the Rayleigh averaging strategy used in the analysis of [11].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 [17] (as briefly discussed in section 2.5).
We analyze the 3D problem with the eigenvalue method for different scenarios presented in [16], and collect the results in table 6, where we report for different geometries and different fracture apertures the unstable modes, and compare them with the convection modes predicted by [16].
As in the two-dimensional case, the instability thresholds are identical to the ones obtained in [16], apart from a slight misprediction for scenario .
In all scenarios furthermore, the eigenfunctions clearly follow the classification proposed in [11], corresponding to either interfracture or intrafracture modes.
Case
Mode [16]
Eigenfunctions
—
—
—
—
—
—
—
—
—
—
—
—
—
In particular, for scenario 6 we observe that (i) the interfracture mode is first unstable one appearing at (ii) already at , different intrafracture modes are available and dominate over the interfracture mode:(iii) 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 and .
Although our results seem to confirm the results of [11], where the intrafracture mode is also identified as dominant,they do not contradict the results in [16].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 (figure 16), even modes with decaying behaviour near equilibrium may turn out to be dominant at steady state.
4.5 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 (20) (including LU factorization of )(c) computation of largest eigenvalues.Step (c) can itself be subdivided in the two costly operations (i) matrix vector products 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 4.2, 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 5).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 7 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 methodis dominated by the assembly of the problem Jacobian in the two-dimensional casesand 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 and , 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 of the Krylov basis in the Krylov-Schur algorithm:for every iteration of the algorithm, the cost of orthogonalization scales with while the cost of the matrix-vector products scales with .Tuning for the problem at hand may be a possible option for improving the efficiency of the method.
Scenario
direct (s)
init
assembly
linsolve
eig (s)
init
lu
matvec
ortho
8192
127.1
0.03
0.70
0.15
22.5
0.17
0.02
0.13
0.23
8192
186.7
0.02
0.74
0.16
23.0
0.24
0.03
0.11
0.18
10016
249.4
0.03
0.73
0.17
37.0
0.17
0.02
0.29
0.28
10422
260.0
0.03
0.72
0.20
53.0
0.17
0.02
0.33
0.34
5000
125.0
0.01
0.77
0.13
20.0
0.22
0.03
0.05
0.08
5000
87.9
0.03
0.75
0.12
13.9
0.21
0.02
0.09
0.15
11379
160.1
0.08
0.68
0.18
165.5
0.07
0.00
0.40
0.42
12073
137.3
0.11
0.59
0.15
107.5
0.21
0.02
0.36
0.28
38400
1053.9
0.02
0.21
0.74
172.2
0.18
0.11
0.52
0.13
38400
886.1
0.04
0.22
0.72
191.7
0.22
0.12
0.48
0.12
In figure 25 we illustrate how the number of matrix-vector products changes with the number of mesh elements .The results show that depends weakly () on the size of the problem (although a lot of variability remains).Note that however the cost of each matrix-vector product scales as : the cost of each matrix-vector product scales as : the presence of in (20) makes the matrix dense.The scaling estimate for the total cost of the matrix-vector products is thus
The estimate is consistent with the data reported in table 7, for which a least-squares estimate gives .
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 26 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.
5 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 problemwe 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 4.2 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 [17] 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 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 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 speedmay enable a more systematic (or even statistical) study of how global properties of fracture networks such