Background theory
Governing equations
Numbat implements the Boussinesq approximation to model density-driven convective mixing in porous media. To reduce the computational burden, only a single fluid phase is considered. This is a simplification to the actual physical process, where a gas phase may be present. This simplification is often used in practice, see Emami-Meybodi et al. (2015) for a discussion about the use of this simplifying assumption.
The more complicated two-phase model can be implemented using the porous_flow module
The governing equations for density-driven flow in porous media are Darcy's law (1) where is the velocity vector, is permeability, is the fluid viscosity, is the fluid pressure, is the fluid density as a function of solute concentration , is gravity, and is the unit vector in the direction.
The fluid velocity must also satisfy the continuity equation (2) and the solute concentration is governed by the convection - diffusion equation
(3) where is the porosity, is time and is the diffusivity.
Darcy's law and the convection-diffusion equations are coupled through the fluid density, which is given by (4) where is the equilibrium concentration, and is the increase in density of the fluid at equilibrium concentration.
The boundary conditions are (5) which correspond to impermeable boundary conditions at the top and bottom boundaries, given by and , respectively, and a saturated condition at the top boundary.
Initially, there is no solute in the model (6)
Numbat solves Eq. 1 and Eq. 3 with density coupled to concentration as in Eq. 4.
Dimensionless formulation
The governing equations can also be solved using a streamfunction formulation in 2D and a vector potential formulation in 3D. As a result, we shall consider the two cases separately.
2D solution
If we consider an anisotropic model, with vertical and horizontal permeabilities given by and , respectively, we can non-dimensionalise the governing equations in 2D following Ennis-King and Paterson (2005). Defining the anisotropy ratio as (7) we scale the variables using (8) where refers to a dimensionless variable. The governing equations in dimensionless form are then (9) (10) (11) where we have dropped the hat on the dimensionless variables for brevity.
The dimensionless boundary conditions are (12) where is the Rayleigh number, defined as (13) In this form, the Rayleigh number only appears in the boundary conditions as the location of the lower boundary. Therefore, can be interpreted in this formalism as a dimensionless model height, and can be varied in simulations by simply changing the height of the mesh.
Finally, the dimensionless initial condition is (14) For isotropic models, where and hence , we recover the dimensionless equations given by Slim (2014).
The coupled governing equations must be solved numerically. To simplify the numerical analysis, we introduce the streamfunction such that (15) This definition satisfies the continuity equation, Eq. 10, immediately.
The pressure is removed from Eq. 9 by taking the curl of both sides and noting that for any , to give (16) where we have introduced the streamfunction using Eq. 15.
The convection-diffusion equation, Eq. 11 becomes (17) The boundary conditions become (18) while the initial condition is still given by Eq. 14.
3D solution
We now consider the case of a three-dimensional model. For simplicity, we consider the case where all lateral permeabilities are equal (). The governing equations for the 3D model are identical to the 2D model. In dimensionless form, they are given by Eq. 9 to Eq. 11, with boundary conditions given by Eq. 12, and initial condition given by Eq. 14.
To solve these governing equations in 3D, a different approach must be used as the streamfunction is not defined in three dimensions. Instead, we define a vector potential such that (19) It is important to note that the vector potential is only known up to the addition of the gradient of a scalar as (20) as for any scalar . This uncertainty is referred to as guage freedom, and is common in electrodynamics. Taking the curl of Eq. 9 and substituting Eq. 19, we have (21) where we have again used the fact that . If we choose to specify the guage condition, this simplifies to (22) As shown in E and Liu (1997), is satisfied throughout the domain if (23) The governing equations are then (24) (25) where the continuity is satisfied automatically because for any .
Finally, it is straightforward to show that in order to satisfy and , which means that the vector potential has only and components, (26) and therefore the fluid velocity is (27) Note that if there is no dependence, Eq. 24 and \eqref{eq:convdiff3d} reduce to (28) It is simple to show that and at are only satisfied if in the entire domain. In this case, the governing equations reduce to the two-dimensional formulation, as expected.
References
- W. E and J. G. Liu.
Finite difference methods for 3d viscous incompressible flows in the vorticity-vector potential formulation on nonstaggered grids.
J. Comp. Phys., 138:57–82, 1997.[BibTeX]
- Hamid Emami-Meybodi, Hassan Hassanzadeh, Christopher P. Green, and Jonathan Ennis-King.
Convective dissolution of CO$_2$ in saline aquifers: Progress in modeling and experiments.
International Journal of Greenhouse Gas Control, 40:238–266, 2015.[BibTeX]
- J. Ennis-King and L. Paterson.
Role of convective mixing in the long-term storage of carbon dioxide in deep saline aquifers.
SPE J., 10:349–356, 2005.[BibTeX]
- A.C. Slim.
Solutal-convection regimes in a two-dimensional porous medium.
J. Fluid Mech., 741:461–491, 2014.[BibTeX]