A multiscale modeling study for the convective mass transfer in a subsurface aquifer
Abstract
Quantitative and realistic computer simulations of mass transfer associated with CO disposal in subsurface aquifers is a challenging endeavor. This article has proposed a novel and efficient multiscale modeling framework, and has examined its potential to study the penetrative mass transfer in a CO plume that migrates in an aquifer. Numerical simulations indicate that the migration of the injected CO enhances the vorticity generation, and the dissolution of CO has a strong effect on the natural convection mass transfer. The vorticity decays with the increase of the porosity. The time scale of the vertical migration of a CO plume is strongly dependent on the rate of CO dissolution. Comparisons confirm the near optimal performance of the proposed multiscale model. These primary results with an idealized computational model of the CO migration in an aquifer brings the potential of the proposed multiscale model to the field of heat and mass transfer in the geoscience.
keywords:mass transfer multiscale modeling porous media CO plume subsurface flow aquifer wavelet method numerical modeling.
List of symbols  

Concentration of CO  
Reference concentration of CO  
Diffusion coefficient ()  
Microscopic length scale ()  
Acceleration due to gravity ()  
Vertical length scale of the aquifer ()  
Permeability ()  
Horizontal length scale of the aquifer ()  
Molar mass of CO ()  
Number of nodes on a regular mesh  
Number of nodes on an adaptive mesh  
Buoyancy frequency ()  
Darcy number,  
Reynolds number at macroscale,  
Reynolds number at pore scale,  
Schmidt number,  
Spacetime mean Darcian velocity ()  
Spacetime mean intrinsic velocity ()  
Deviation from in a REV ()  
Deviation from in a REV ()  
Dimensionless spacetime mean intrinsic velocity  
Dimensionless spacetime mean intrinsic concentration  
Cartesian coordinate in the th direction  
Greek symbols  
Ratio of molecular viscosity to effective viscosity  
Solutal expansion coefficient  
Rate of CO dissolution ()  
Representative elementary volume (REV)  
The fraction of REV occupied with fluid  
time step  
Local step size in the direction  
Minimum of for  
Porosity  
Reference density  
Characteristic length scale of a REV  
Error tolerance  
Molecular viscosity ()  
Kinematic viscosity () 
1 Introduction
The natural and mixed convection heat and mass transfer is an important topic of scientific scrutiny in a subsurface flow problem that investigates the disposal of anthropogenic CO from the atmosphere into saline aquifers [33]. Example of such projects include the Utsira sand at Sleipner, Norway [12], the Mt. Simon aquifer in the Illinois basin [26, 4], saline aquifers in the Alberta basin, Canada [38], and the CarrizoWillcox aquifer in Texas (CWT) [29]. Seismic data from the Sleipner shows a marked increase in the CO flux near the reservoir top (see Fig 1 of [12]), which suggests further investigation of the rate of vertical migration and progressive development of CO plumes (e.g. [12]). Chadwick et al [12] hypothesized that this marked variability may be due to the multiscale natural convection mass transfer associated with the nonDarcian plume migration through numerous pathways in the aquifer. In the same vein, multiscale processes associated to a plume migration in the CarrizoWillcox aquifer was investigated numerically by Pruess and Nordbotten [37] using a classical macroscopic Darcian model along with a subgrid scale parameterization scheme (albeit the exact form of the scheme was not outlined with full details in [37]). Clearly, a complete understanding of the multiscale mass transfer mechanism in aquifers remains an active research area in the field of geoscience and reservoir engineering [29, 37, 27, 40, 16, 34, 15]. Hence, there are increasing interests on extending the adaptive mesh and multiscale finite volume/element methods (AMR based methods) for high performance numerical simulations of flow and transport in saline aquifers (e.g., see, [16, 20, 24, 30, 31, 23, 34, 15]).
The present article has investigated the development of a novel adaptive wavelet multiscale modeling and simulation methodology for studying the nonDarcian flow and transport through aquifers. One objective of this article is to investigate an effective methodology to minimize computation work and to improve the accuracy for transport problems in aquifers that deal with multiscale phenomena, in comparison to the commonly used classical numerical methods. For example, some authors verified with a classical method that a between and is necessary to capture multiscale features (e.g., see [38, 37, 34]), for which, the number of the grid points is about in the domain (e.g. [38]), and hence, would be far beyond the limit of modern computers for such a faithful simulation in a vertical cross section () of the CWT [37]. In contrast, this article investigates a multiscale methodology to put the computational effort locally in the physical domain (see [44, 2, 3, 1, 47, 46] for a technical details), where it is necessary to adapt the computational work to fieldscale geological structures and complex physical processes in aquifers.
In the same vein, AMR based classical adaptive and mutiscale methods require to find an individual error monitoring scheme for each new applications (e.g. [23, 24, 34]). These methods were developed for steady elliptic problems, and their extension to multiscale flow and transport in aquifers is being investigated by a number of other authors (see the recent work of Künze et al [24]). Note also that an extremely small time step () is needed when these AMR based methods explicitly simulate natural and mixed convection heat and mass transfer in aquifers. In contrast, the present article investigates a multiscale methodology that captures multiscale physics with a more robust wavelet based computation, where the error is controlled a priori according to a prescribed tolerance () [47, 46], and the time step () can be chosen independent of the spatial grid () [2, 3]. In addition, the present research employs the volumetric mean of the timeaveraged conservation laws with respect to a representative elementary volume (REV) so that the nonDarcian multiscale features (e.g. [12]) in a flow may be resolved with a nonDarcian multiscale model (see [45, 14] and the refs therein). This approach provides a natural framework for adopting appropriate parameterization for phenomena which are not directly resolved with the mean conservation laws (e.g. [37]). Here, the effect of CO dissolution on the gravitational segregation of a CO plume has been parameterized. Most importantly, the size of the REV has been adapted dynamically to the local variation of the CO mass fraction. To the best of knowledge, the present development is a first time investigation for the multiscale natural convection mass transfer in a saline aquifer, and it complements the growing trend in computational modeling of the natural/mixed convection phenomena in subsurface flows as well as in other fields [23, 31, 34, 36, 47, 46]. Following are few remarks on the proposed multiscale methodology.
Remarks:

The flow and transport at the reservoir scale is approximated with the classical Darcian approach. The multiscale mass transfer has been resolved with a nonDarcian approach. Subgrid scale parameterization schemes for both the momentum diffusion and the CO dissolution have been adopted.

A multilevel algorithm has been developed for the nonlinear coupling between the mass and momentum transfer processes in aquifers, which occurs at different temporal and spatial scales.

Spatial differential operators are discretized using significant wavelets, where each wavelet represents the change in scale near a grid point, and multiscale features are captured with a multiscale wavelet theory.

There are two important computational benefits. First, the number of grid points is significantly small compared to the number of grid points needed for a classical numerical method if the same level of accuracy is desired. Second, if increases, then the CPU time increases approximately linearly. Note that for a given tolerance , may increase due to change in gradients of solution.

The error for such a simulation is controlled in both space and time according to the prescribed tolerance () and the maximum allowable CFL number, which is verified with a large number() of numerical experiments.
Section 2 summarizes the new developments towards capturing the multiscale features in a subsurface plume migration. The adaptive multiscale methodology has been presented in section 3. Representative results from a series of numerical experiments have been presented in section 4. Finally, section 5 discusses the potential extension of the present development in the field of computational heat and mass transfer analysis.
2 The proposed multiscale modeling framework
This model aims to simulate the natural convection mass transfer phenomena associated with CO storage and migration in a typical aquifer, such as the CarrizoWilcox aquifer, Texas, US (a further details in section 4.2.1). In the following development, the aquifer height, , is the length scale, and a typical flow speed, , in the aquifer is the velocity scale. A vertical cross section of the aquifer is assumed to be simulated, which is a good approximation to represent overall multiscale features of the problem(see, [38, 37]). The boundary conditions assume impermeable caproks on the bottom and top boundaries as well as open lateral boundaries. The aquifer is assumed a homogeneous and isotropic porous medium.
The present model diagnoses the flow and transport at the reservoir scale from the classical Darcian approach. The multiscale mass transfer has been resolved with a nonDarcian approach, where the REVs are resized locally in order to adapt the effective numerical resolution to the localized features of mass transfer. Subgrid scale parameterization schemes for both the momentum diffusion and the CO dissolution have been adopted.
2.1 Classical upscaling approach for diagnosing macroscopic transfer
The pore scale (fine scale) heat, mass, and momentum transfer is replaced to the reservoir scale (macroscale) (e.g. [34, 6]), where
are the time and volume averages of the velocity field, respectively. The Darcian velocity is a spacetime mean with respect to the REV () and the time scale, , which contains both the solid and the fluid (see [6] for details). Therefore, the natural convection mass transfer in an aquifer is typically studied with the following upscaling representation (e.g., [34, 17, 42, 21, 39, 35])
(1) 
(2) 
In this model, the momentum transfer is in a steady state, both the advection and the deviatoric stress make a negligible contribution to momentum transfer, the drag between the porous medium and the fluid is assumed linearly proportional to the pressure gradient, and the Reynolds number, does not exceed a value about (e.g., see Chapter 5, [6]).
As discussed in the introduction, using the macroscopic upscaling (12), faithful simulations of the convective mass transfer in an aquifer often requires a high numerical resolution (e.g. [38, 37, 34]) because the characteristic length scale () of the REV is typically at the order of . As demonstrated schematically in Fig 1, needs to be sufficiently small in order for so that computations on grid points would represent the actual mass transfer along streamlines. In the following two sections, a multiscale representation of the fine scale mass and momentum transfer has been studied, using the most robust volumetric mean of the timeaveraged conservation laws [6, 10, 9, 48].
2.2 Multiscale intrinsic mass transfer model
The intrinsic spacetime average is
where is the volume of the fluid contained in a REV. Note the difference between and , where [22]. The time and volume averages are assumed to commute; i.e,. one may apply the temporal average into a volume averaged variable, and viceversa. Using the intrinsic average , the spacetime transmissivity in a REV may be resolved, which is neglected in the classical macroscopic model (12).
To illustrate the multiscale framework, we begin with the decomposition
at three scales, where the first component, , captures the mean discharge per unit area, but does not represent a details of the flow. The second component, , adds the missing details into such that and . The third component, , represents a further details that satisfies and . Clearly, the spacetime mean intrinsic velocity resolves an additional details with respect to the classical mean Darcy velocity . Similarly, the spacetime intrinsic mean of the concentration of CO can be obtained. In order to simplify the symbolic representation, is used for the dimensionless mean velocity , and for the dimensionless mean concentration in the rest of this article, where all quantities are assumed uniform in the direction. Note that this twodimensional assumption is an idealization for the radial symmetry of the dynamics of the CO plume (e.g., [31, 37]), and has been adopted to aid the investigation on the proposed multiscale model development.
The derivation of the present intrinsic multiscale model for the mass transfer is similar to what was detailed by Breugem and Rees [10] for the heat transfer, and by Lage et al [25] for the momentum transfer (see also [6, 9, 48, 22, 30]). However, in the heat transfer model, an independent diffusion equation was considered in [10] to account for thermal conduction through the solid phase. In the present mass transfer model, the dissolution of the invaded phase has been parameterized through the averaging process. At a depth below m (see, [7]), the geothermal effects may be compensated by the geopressure gradient (see, [27] and chapter 2 of [14]), which is further illustrated by Pruess and Nordbotten [37]. The pressure gradient terms in eqs. (45) accounts for the geothermal pressure gradient. Since studies indicate a weak dependence between the molar volume of dissolved CO and density of the binary mixture, the Boussinesq approximation is reasonable [5, 27, 14, 34]. Following the derivation of Lage et al [25], the macroscopic upscaling (12) has been replaced with the multiscale upscaling (36):
(3) 
(4) 
(5) 
(6) 
In this model, the Schmidt number is , the macroscopic Reynolds number and the Darcy number are defined by and , respectively, where and . Clearly, if , we have , which corresponds to a macroscopic model (12).
As described in Appendix A, the spacetime average of the nonlinear advection of momentum takes the form
which resulted into the additional term
of which, the first two components have been neglected in the present laminar mass transfer model, and the last term can be parameterized with the Brinkman model [11], which – in the dimensionless form – has taken the form of the second term on the right hand sides (rhs) of both (4) and (5). The third term on the rhs of (4, 5) has appeared due to the intrinsic average of the viscous and pressure stress, which models the density of the pressure drag and skin friction,
when the Reynolds number is small, i.e., [14].
Similarly, the nonlinear advection of the mass flux has resulted into the additional term
of which the first two components have been neglected and the last component has been parameterized to model the effect of dissolution with the first term on the right hand side of (6).
2.3 The dissolution of CO through the dispersive mass flux
Literature review does not indicate a common approach for modeling the effect CO dissolution into the resident saline, despite there are few attempts (e.g. see [37]). If a CO plume proceeds upward from an isolated source, the saline density at the interface increases by about  [34]. To a first order approximation, the global density of the background environment adopts a vertically decreasing profile when the plume migrates upward, albeit more specific observational data would confirm the actual density profile (e.g. see also the density distribution presented by Pruess and Nordbotten [37]).
The present work has proposed a simple model
for the mass dispersion associated with the microscopic spacetime variation of mass and momentum, where a steady state horizontally homogeneous background concentration has been assumed since the molecular mass of CO is larger than that of brine (see [40, 20, 34]). With this model, one can define a Buoyancy frequency by to characterize the effect of dissolution, where is the molar mass of CO and is a reference density. Here, corresponds to a situation with ; i.e., in this case, CO has been accumulated near the reservoir top, where the dissolution would result into gravitational fingers studied by Pau et al [34]. corresponds to ; i.e., CO dissolution is now associated with a vertically migrating plume [12]. Accordingly, a Froude number is given by
and using data from the CarrizoWilcox aquifer, Texas, e.g., and [37, 29], we estimate that corresponds to and . Clearly, results into , and the mass transfer analysis for varying exploits the effect of the CO dissolution.
2.4 The natural convection mass transfer as a function of
The dominant mechanism for the onset of background dissolution on the natural convection mass transfer during the migration of a CO plume has now been studied with a dimensional analysis, which explains the effect of the variation of the Froude number. To estimate the order of magnitude of each term in the mass and momentum conservation laws (see Breugem and Rees [10]), consider the solutal Grashof number,
the Schmidt number, , the horizontal length scale, , and the vertical length scale, . Note the large aspect ratio of typical aquifers (e.g. the CarrizoWilcox aquifer, Texas, US [29]) and the order of magnitude for the inertial term, . Using horizontal and vertical length scales, and , respectively (e.g. [20]), where , as well as a fixed , and in the limit of , we obtain the following dimensionless linear system of PDEs
The system is independent of the length and time scales, as well as of other dimensional parameters, and hence, exhibits a selfsimilar solution. The existence of a self similar solution indicates that the vertical length scale can be determined from the dimensional parameters , , , , and those appear in the dimensional system of equations. A dimensional reasoning can be used to define the vertical length scale as a function of the remaining other parameters; i.e.,
which gives an aspect ratio between the vertical scale and the horizontal scale:
Clearly, in the limit of for a fixed and , the vertical length scale extends to infinity; i.e., . This indicates that a nonzero gradient, , tends to reduce vertical migration of the CO plume. For too small a value of , the vertical length scale is also too small; the plume will not continue to migrate vertically upward for small . This dimensional analysis has a good agreement with the numerical simulations presented in section 4.3.2.
3 The adaptive wavelet multiscale simulation methodology
To outline the proposed numerical method, eqs (46) has been written in the following compact form:
(7) 
where represents , , or , and () represent the right hand sides of eqs (46). Since the velocity () and the concentration () depends on each other simultaneously in a natural convection mass transfer application, a fully implicit time integration scheme has been adopted, where , , and are computed simultaneously.
3.1 Time integration
A fractionalstep method, originally proposed by Chorin [13], has been applied to solve (7) for ,
and it takes the following symbolic form
(8) 
where superscripts and mean the present time step and the first fraction of the next time step, respectively.
In the second fraction of a time step, eq. (3) is satisfied, , such that the macroscopic model (1) is approximately diagnosed from
(9) 
where . Setting from the nonlinear algebraic system (8), is also diagnosed from the elliptic problem (9), where (8) and (9) must be solved with efficient iterative solvers at each time step.
Since (8) is a nonlinear system, the classical Newton’s method
such that
must evaluate the matrixvector product at every th iteration with operations, where is the Jacobian matrix and is the error to be found. In order to reduce this high cost to , let us consider the Frechet derivative
where is a small real number. The performance of this approach – known as the Jacobian free NewtonKrylov (JFNK) algorithm – was verified for multiphysics simulations. However, the improvement in the operation count is paid off by requiring a preconditioner, which is a serious drawback for extending the JFNK to the simulation of heat and mass transfer. In the present work, we study an alternative, where the multiscale wavelet method captures the multiscale physics, and the Newton’s method along with the Frechet derivative is used to reduce the residual of (8) by a certain fraction at each level of the present multiscale wavelet based solution methodology.
3.2 The multiscale wavelet methodology
The wavelet method (see [2, 28, 43]) captures the multiscale physics, using the best terms of the multiscale decomposition
(10) 
according to a prescribed error tolerance on the magnitude of , where the error for resolving mass and/or energy is . In (10), the first term represents a sampling of (, , or ) on a coarse grid of points, and the second term models the additional details of , which is not captured by the first term. In D, with the exception of east and north boundaries, each coarse grid sampling corresponds to three details data (for ). For example, if the REV is divided by a factor of in each direction, for the th vertex , we have three neighbors , , and , where the corresponding detailed information is stored. When a REV is not sufficient to resolve on the th grid point, at least one will have a magnitude larger than , thereby providing with a natural framework for adapting the size of the REV to the local physical property. Therefore, one would begin with REVs on a coarse grid, perform a wavelet analysis, and recursively resize only those REVs, where wavelet coefficients are large. The wavelet transform may be computed with the Wavelet toolbox of Matlab without knowing the explicit information of and . Note that the grid adaptation is automatic with the wavelet method [44].
To illustrate the benefits of the multiscale wavelet representation, consider a prescribed concentration that decays exponentially with respect to a circle of radius . Taking the wavelet transform of , and recursively adapting the grid until all new wavelet coefficients satisfy a tolerance , an adaptive wavelet grid has been obtained, which is presented in Fig 2. Note that an obstacle of size has been placed at the center of the domain, showing that the wavelet transform can be computed on complex domains.
4 Results and verification
4.1 Verification results
The shear driven or natural convection transport of CO in an aquifer is an idealized model, and is useful for computational verification, where CO moves horizontally just below the impermeable caprock, or vertically after it has been injected through an injection well. Adapted from [38, 32], a shear driven case and a natural convection case have been shown schematically in Fig 3. Note that the shear driven case has been chosen for the availability of reference data so that the numerical model can be quantified. In the next two simulations, the parameters , , , , have been adapted from [48].
4.1.1 A shear driven mass transfer model for
The shear driven mass transfer at is a benchmark example, where the momentum transfer is independent of the mass transfer. To study the performance of the proposed multiscale model, results of simulations of a shear driven mass transfer case with have now been summarized, where for each of the CFL numbers are , , , , , , the tolerance values are , , , , . In the present multiscale model, controls , as well as – the length scale of the smallest REV, where the grid points are adapted dynamically. To utilize the full advantage of this adaptivity, the maximum time step, , has also been adapted dynamically based on CFL = . Clearly, an increase of the CFL number (or ) would increase the global temporal (or spatial) error. These simulations provide an understanding of how the spatial and temporal errors are controlled globally in the present model. Note the use of a large CFL number (CFL=) is a distinct feature of the present methodology, and is an important contribution to the field of computational heat and mass transfer.
As demonstrated in Fig 4(), the variation of the spatial error tolerance and the CFL number, confirms the spacetime error control of the proposed model. This velocity profile has a good agreement with the one, which was presented by Yang et al [48] and Guo and Zhao [19]. Fig 4() shows that the computational degrees of freedom, , varies like , which means that does not increase linearly if is decreased. Note that a direct numerical simulation of heat and mass transfer phenomena with a classical finite difference (FD) model, would increase linear with a reduction of the error. The efficiency of the present model can be seen from the comparison in Fig 4(). Further more, with the error tolerance, , the present model needs , which is about and of the grid points required by the simulations of Ghia et al [18] and Botella and Peyret [8], respectively. Although the physical setting of the present simulation is quite different than that of refs [18, 8], the similarity of dimensionless parameters allows us for a brief comparison in order to provide with a hint to the potential of the present model.
The laminar boundary layer thickness just below the top of the aquifer has been compared with the theoretical estimate, , in Fig 4(). A good fit between and indicates that the laminar boundary layer has been resolved more accurately at higher . This idealized simulation at exploit the potential of the present model for simulating the best solution using the least number of the degrees of freedom, .
4.1.2 A natural convection case for
An idealization of the natural convection mass transfer in the aquifer has been simulated, where the boundary conditions, and on and , represents the injection and removal of CO across and , respectively. Other boundary conditions are on and , on and . The velocity is assumed zero on all other boundaries. Under similar conditions, Nordbotten and Celia [32] studied mixed convection mass transfer in an aquifer, and showed that the migration of CO occurs beneath the top boundary of the aquifer when natural convection dominates over the forced convection. Clearly, at large , when the natural convection is dominant, the numerical simulation of the mass transfer in an aquifer would require high spatial resolution in order to resolve the solutal boundary layer. Hence, an estimate for as a function of provides with a good understanding for the efficiency of the present multiscale model. Here, we summarize representative simulations with values of (, , , , , , and ), which show that increases approximately as . In Fig 5, this result is also compared with the estimate from the classical numerical simulation.
4.2 Simulation of a vertically migrating CO plume
4.2.1 The CarrizoWilcox aquifer in Texas (CWT)
Over a period of years, the injection of approximately supercritical CO per year into the central section of the CWT would store about one fifth of the CO emissions in Texas [29, 20]. The CWT formation is about deep, and extends for a length of about at an angle with the horizontal, reaching a depth of about beneath the earth’s surface. A 2D vertical section of the CWT beneath a horizontal caprock at a depth of , which is thick, and long, has been simulated. The length to width aspect ratio of the model region is , which is different than the actual aspect ratio of the CWT. The temperature of the aquifer is assumed at the initial value. A wide CO plume is assumed instantaneously at the center of the bottom boundary. The migration of CO takes place under the action of gravity. Initially, both the resident brine and the invaded CO are assumed at rest. The normal components of the gradient of mass and momentum fluxes are assumed zero on the lateral boundaries for all time.
In Fig 6, simulated vertical migration of a CO plume and the corresponding adapted grid points are presented at year. Clearly, the REVs have been adapted dynamically to resolve sharp gradients of the plume with respect to the error tolerance . The minimum size of a REV is given by and near the center of the plume, and the maximum size of the same is given by and away from the plume. The simulation needed an average at each time step.
As depicted in Fig 6(c), the plume remains approximately in a region that extends from to in the horizontal direction with respect to the center of the domain. However, only in this part of the domain, one might get the qualitatively equivalent results by solving the multiscale model (36) with a classical numerical method, using the uniform resolution of and . In that case, at least grid points would be needed. In comparison to this estimate, the present methodology (with ) is about efficient, albeit a direct comparison with a previous simulation would be realistic.
4.2.2 A comparison of the computational gain with a TOUGH2 simulation
The recorded CPU time for the present simulation is approximately per grid point on a Dell T7400 workstation. A useful reference on the computational effort that is needed for a similar simulation is given by Pruess and Nordbotten [37], where the reported CPU time is about per grid block using the parallel code, TOUGH2, with processors in a Dell 5400 workstation. This is equivalent to about per grid block for a single processor computation. The CPU time comparison with this particular TOUGH2 simulation exploits the advantage of computational complexity, which is one important development of the present multiscale model.
Since varies at each time step due to adaptivity, we have recorded the CPU time for advancing the solution each time step, and two representative results are presented in Fig 5() and Fig 7, showing that the CPU time is approximately proportional to . These results explore the promise of this development.
4.3 Analysis of natural convection mass transfer in an aquifer
4.3.1 Penetrative mass transfer
Ruckenstein [41] studied a generalized penetration theory for mass transfer in the vicinity of a fluid–fluid interface. In this section, we have briefly studied the multiscale nature of the penetrative mass transfer and the associated vorticity generation in a subsurface aquifer. In order to study the vorticity generation, we have simulated natural convection mass transfer at various values of and .
In the limit of , a first order estimate for the vorticity as a function of concentration is , which indicates that a counter clockwise vorticity is generated in the region that is to the immediate left of the plume when the plume migrates upward [40]. Similarly, a clockwise vorticity is expected in the region that is to the immediate right of the plume. These clockwise and counter clockwise vortices interact with the porous structure, accelerate horizontal migration, and are responsible for transferring mass and momentum from one scale to the other. In Fig 8, the plume, the vorticity, and the adapted grid are presented for the porosities, and , where . A complicated multiscale dynamics is observed. For the higher porosity, only of a REV is occupied by the fluid. A reduction of the porosity would decrease the volume of fluid faction in the REV, thereby strengthening the interaction between the fluid and the porous structure. For , the mass transfer is accompanied with local sharp spatial gradients at multiple length scales, where grid points are needed because REVs are resized to resolve sharp gradients of the plume. If the porosity is doubled to , the rate of mass transfer is enhanced by smoothing out the spatial gradients, and thus, the number of grid points is reduced to . A better understanding of the vorticity field can be quantified with simulations at various values of .
For the fixed , let us present the time evolution for the mean vorticity for a variation of the Darcy number, . Since the variation of is between and , these simulations explore the sensitivity of perturbations introduced by the porous structure into the penetrative mass transfer. The growth of the mean vorticity, and its time evolution has been presented in Fig 9. Clearly, the Darcy number affects both the maximum of the mean vorticity and its amplitude of the fluctuation. Note that the natural convection mass transfer enhances the vorticity at the lowest of these values of , which indicates that the penetrative mass transfer in a saline aquifer is accompanied with vorticity.
4.3.2 The effect of CO dissolution into the natural convection mass transfer
In the present model, the dissolution of CO into the resident saline has been characterized by the Froude number . In Fig 10, simulated distribution of CO has been presented at for three representative values of the Froude number, , , and . The plume migrates vertically until it reaches the impermeable caprock at the top of the reservoir when , when the brine is saturated with CO such that . When , the plume migrates through a background environment with . Comparing Fig 10(a) with Fig 10(c), we see that the plume has traveled less than half way in the vertical direction when .
This simulation indicates a marked variability in the vertical migration of CO. These twodimensional idealized simulations hint on the potential of the present model for analysing natural convection mass transfer and other related phenomena associated with CO storage in aquifers. Further analysis with three dimensional simulations would be useful for explaining the multiscale and complex phenomena associated with the carbon capture and storage program.
5 Summary and discussion
5.1 Conclusion
This article introduces a multiscale modeling and simulation approach for the natural convection mass transfer in an aquifer. In order to tackle the challenges of multiscale phenomena, we have adopted the classical volume averaging technique to model a full details of the multiscale transport. Taking average heat and mass transfer with respect to a REV leads to a general framework for parameterizing the effect of multiscale phenomena which is not resolved with the averaging process. A wavelet based multiscale simulation methodology has been studied to dynamically resize the REV so that localized mass transfer can be resolved efficiently. A brief summary of the present development and key findings have been outlined below.

Using the spacetime average of the first principle conservation laws with respect to a representative elementary volume and a representative elementary time scale, this work has studied the governing equations so that heat, mass, and momentum transfer can be captured in a range of multiple length and time scales.

Results from a set of representative numerical experiments at and show that the spacetime error can be controlled with a priori prescribed error tolerance and CFL number. This a priori error control is an important contribution that would benefit future developments in the field of computational heat and mass transfer analysis.

Simulations with an idealized natural convection mass transfer observes that the number of computational degrees of freedom varies like . Clearly, exploiting multiscale physics in the computational model has the potential to compute the most significant proportion of the flow using a near optimal computational effort. Further more, all simulations have indicated that the CPU time remains approximately proportional to , which means that the amount of computational work remains proportional to the amount of actual physical change in the system. To the author’s opinion, this property of the present model is a distinct feature with respect to classical approaches those are commonly used in the computational heat and mass transfer analysis.

Present simulations with varying show a marked variation in the vertical migration of a CO plume, i.e. the associated time scale, which suggests that the effect of dissolution of CO in saline has a dominant role on the mass transfer mechanism. This would help to explain the dynamics of the subsurface CO plume in various storage facilities.

Present simulations observes that the vertical migration of a plume enhances the vorticity generation, where the dissolution of CO has a strong influence on the vertical time scale of a plume.

To the best of knowledge, this article, for the first time, has extended the second generation wavelet based adaptive technique to the field of subsurface flow and transport modeling. Hence, the present works adds a novel technology to the growing interests of multiscale modeling in the field of computational transport in aquifers.
5.2 Discussions
There are several possibilities to extend the present development. This includes the simulation of heat and mass transfer in a large aquifer – such as long CWT or any other similar sites – where the phenomena is fully threedimensional and gravitationally unstable. Future simulations, where the permeability, porosity, and other reservoir characteristics would have been obtained with advanced field technologies, may provide with an understanding of the degree to the proportion of actual multiscale features that can be resolved by the proposed model within the constraint of modern powerful computing resources. In a more realistic D simulation, if is still constrained by the computer power, then the smallest scale of a REV may need to be increased, where further improvements of the parameterization scheme would be necessary. Therefore, the present work leaves potential open questions those may be addressed in the future. However, the present research indicates that D simulation will be benefited greatly if the present modeling approach is extended. Furthermore, the effect of variable permeability, temperature perturbation, and pathway transmissivities have not been examined in the present article. These works are currently underway.
Appendix A
The spacetime double decomposition methodology
This section outlines the proposed model of resolving multiscale physics of a natural convection mass transfer when a plume migrates after CO has been injected into an aquifer. The methodology aims to capture the average flow and transport with respect to a REV and a representative time scale for the REV, as well as to parameterize the effect of the unresolved flow. Researchers with interest in further mathematical details of the adopted double decomposition methodology may find the works of De Lemos [14] and Lage et al [25] useful. This section provides a brief outline for the decomposition of the inertia term into a resolved proportion and an unresolved proportion, where the intrinsic space average () is applied to time average (), and viceversa (see also, section 2.2).
If we take the time average of a spatially averaged quantity, and use to denote temporal fluctuation of spatial averages, then a spatial mean may be decomposed as
where is the deviation of with respect to the time average. The double decomposition is obtained by taking a spatial decomposition, , which is followed by a temporal decomposition; i.e.,
The following assumptions have been adopted to parameterize the effect of unresolved multiscale physics.

corresponds to momentum flux perturbation associated with the temporal fluctuations of the spatially averaged quantities , which has been neglected. This term needs to be parameterized if a transition to turbulence is important; i.e., if .

corresponds to momentum diffusion associated with the temporal mean of spatial fluctuation (same as the spatial average of temporal fluctuation). In order to allow a smooth transition of flow and transport through pores of a porous medium, Brinkman [11] suggested to incorporate the effect of momentum diffusion. In the present work, this term has been parameterized to incorporate the effect of momentum diffusion.

represents dispersion of momentum due to both time and spatial fluctuation, and can be neglected unless turbulence intermittency is important.
The above approximations, along with the conservation of mass and the Boussinesq approximation, leads to , and hence, the spatiotemporal average of the divergence of the momentum flux (divided by the density) can be derived recursively. First, let us apply the spatial decomposition, , and write
Second, both the mean, and the fluctuation, , are decomposed into a temporal mean and fluctuation, and as a result,
Third, the decomposition of the spatial fluctuation, into a temporal mean, , and a temporal fluctuation, , (i.e. ) results into the final form
Finally, using the average mass conservation,
Note that the above decomposition tells us what physical features of the problem has been modelled and the assumptions made to neglect some of the other physics. Similar approach has been used to derive Eq (6).
acknowledgements
The author acknowledges financial support from the National Science and Research Council (NSERC), Canada. Useful discussions with Prof. Jan Martin Nordbotten (jan.nordbotten@math.uib.no) is also greatly acknowledged. Many thanks to two anonymous reviewers for very useful comments and suggestions. The computational work was done partially with a Dell T7400 Workstation funded by the Industrial Research and Innovation Fund (IRIF), Govt of Newfoundland and Labrador, and partially on the ACEnet (www.acenet.ca) high performance computing cluster.
References
 Alam and Islam [2014] Alam J, Islam MR (2014) A multiscale eddy simulation methodology for the atmospheric ekman boundary layer. Geophysical & Astrophysical Fluid Dynamics 0(0):1–20, DOI 10.1080/03091929.2014.975127, in press
 Alam et al [2012] Alam JM, Kevlahan NKR, Vasilyev OV, Hossain Z (2012) A multiresolution model for the simulation of transient heat and mass transfer. Numerical Heat Transfer, Part B 61:1–24
 Alam et al [2014] Alam JM, Walsh RP, Alamgir Hossain M, Rose AM (2014) A computational methodology for twodimensional fluid flows. International Journal for Numerical Methods in Fluids 75(12):835–859
 Bachu et al [1994] Bachu S, Gunter W, Perkins E (1994) Aquifer disposal of CO: Hydrodynamic and mineral trapping. Energy Conversion and Management 35(4):269 – 279
 Barbero et al [1983] Barbero J, Hepler L, McCurdy K, Tremaine P (1983) Thermodynamics of aqueous carbon dioxide and sulfur dioxide: Heat capacities, volumes, and the temperature dependence of ionization. Can J Chem 61(11):2509–2519
 Bear [1972] Bear J (1972) Dynamics of fluids in porous media. Elsevier (New York)
 Benson and Cook [2005] Benson SM, Cook P (2005) Underground Geological Storage, IPCC Special Report on Carbon Dioxide Capture and Storage, Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, U.K., chap 5, pp 195–276
 Botella and Peyret [1998] Botella O, Peyret R (1998) Benchmark spectral results on the liddriven cavity flow. Computers & Fluids 27(4):421 – 433
 Braga and de Lemos [2009] Braga EJ, de Lemos MJ (2009) Laminar and turbulent free convection in a composite enclosure. International Journal of Heat and Mass Transfer 52(3 4):588 – 596
 Breugem and Rees [2006] Breugem W, Rees D (2006) A derivation of the volumeaveraged boussinesq equations for flow in porous media with viscous dissipation. Transport in Porous Media 63(1):1–12
 Brinkman [1949] Brinkman H (1949) A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Applied Scientific Research 1(1):27–34
 Chadwick et al [2009] Chadwick R, Noy D, Arts R, Eiken O (2009) Latest timelapse seismic data from sleipner yield new insights into CO2 plume development. Energy Procedia 1(1):2103 – 2110
 Chorin [1968] Chorin A (1968) Numerical solution of navierstokes equation. Math Commp 22:745–762
 De Lemos [2006] De Lemos M (2006) Turbulence in Porous Media: Modeling And Applications. Chemical, Petrochemical & Process, Elsevier Science Limited
 Elenius and Johannsen [2012] Elenius M, Johannsen K (2012) On the time scales of nonlinear instability in miscible displacement porous media flow. Computational Geosciences 16:901–911
 Farajzadeh et al [2007] Farajzadeh R, Salimi H, Zitha PL, Bruining H (2007) Numerical simulation of densitydriven natural convection in porous media with application for CO2 injection projects. International Journal of Heat and Mass Transfer 50(25 26):5054 – 5064
 Farajzadeh et al [2009] Farajzadeh R, Zitha PLJ, Bruining J (2009) Enhanced mass transfer of CO2 into water: Experiment and modeling. Industrial & Engineering Chemistry Research 48(13):6423–6431
 Ghia et al [1982] Ghia U, Ghia K, Shin C (1982) HighRe solutions for incompressible flow using the navierstokes equations and a multigrid method. Journal of Computational Physics 48(3):387–411
 Guo and Zhao [2002] Guo Z, Zhao TS (2002) Lattice Boltzmann model for incompressible flows through porous media. Physical Review E 66(3):036,304–9
 Hesse et al [2008] Hesse MA, Orr FM, Tchelepi HA (2008) Gravity currents with residual trapping. Journal of Fluid Mechanics 611:35–60
 Homsy [1987] Homsy GM (1987) Viscous fingering in porous media. Annual Review of Fluid Mechanics 19(1):271–311
 Hsu and Cheng [1990] Hsu C, Cheng P (1990) Thermal dispersion in a porous medium. International Journal of Heat and Mass Transfer 33(8):1587 – 1597
 Jenny et al [2006] Jenny P, Lee S, Tchelepi H (2006) Adaptive fully implicit multiscale finitevolume method for multiphase flow and transport in heterogeneous porous media. Journal of Computational Physics 217(2):627 – 641
 Künze et al [2013] Künze R, Lunati I, Lee SH (2013) A multilevel multiscale finitevolume method. Journal of Computational Physics 255(0):502 – 520
 Lage et al [2002] Lage JL, de Lemos MJS, Nield DA (2002) Modeling turbulence in porous media. Transport phenomena in porous media 2:198–230
 Leetaru et al [2009] Leetaru HE, Frailey SM, Damico J, Mehnert E, Birkholzer J, Zhou Q, Jordan PD (2009) Understanding CO plume behavior and basinscale pressure changes during sequestration projects through the use of reservoir fluid modeling. Energy Procedia 1(1):1799 – 1806
 Lindeberg and WesselBerg [1997] Lindeberg E, WesselBerg D (1997) Vertical convection in an aquifer column under a gas cap of CO2. Energy Conversion and Management 38, Supplement(0):S229 – S234
 Mallat [1999] Mallat S (1999) A wavelet tour of signal processing. Academic Press, New York
 Nicot [2008] Nicot JP (2008) Evaluation of largescale CO2 storage on freshwater sections of aquifers: An example from the Texas Gulf coast basin. International Journal of Greenhouse Gas Control 2(4):582 – 593
 Nithiarasu et al [1999] Nithiarasu P, Seetharamu K, Sundararajan T (1999) Numerical investigation of buoyancy driven flow in a fluid saturated nondarcian porous medium. International Journal of Heat and Mass Transfer 42(7):1205 – 1215
 Nordbotten [2009] Nordbotten J (2009) Adaptive variational multiscale methods for multiphase flow in porous media. Multiscale Modeling & Simulation 7(3):1455–1473
 Nordbotten and Celia [2006] Nordbotten JM, Celia MA (2006) Similarity solutions for fluid injection into confined aquifers. Journal of Fluid Mechanics 561:307–327
 Nordbotten et al [2005] Nordbotten JM, Celia MA, Bachu S (2005) Injection and storage of CO in deep saline aquifers: Analytical solution for CO plume evolution during injection. Transport in Porous Media 58:339–360
 Pau et al [2010] Pau GS, Bell JB, Pruess K, Almgren AS, Lijewski MJ, Zhang K (2010) Highresolution simulation and characterization of densitydriven flow in CO2 storage in saline aquifers. Advances in Water Resources 33(4):443 – 455
 Peaceman and Rachford [1962] Peaceman DW, Rachford HH (1962) Numerical calculation of the multidimensional miscible displacement. Society of Petroleum Engineering Journal 2:327–339
 Popiolek and Awruch [2009] Popiolek T, Awruch A (2009) An adaptive mesh strategy for transient flows simulations. In: Computational Modeling (MCSUL), 2009 Third Southern Conference on, pp 71 –76
 Pruess and Nordbotten [2011] Pruess K, Nordbotten J (2011) Numerical simulation studies of the longterm evolution of a CO plume in a saline aquifer with a sloping caprock. Transport in Porous Media 90:135–151
 Pruess and Zhang [2008] Pruess K, Zhang K (2008) Numerical modeling studies of the dissolution diffusion convection process during CO2 storage in saline aquifers. Tech. Rep. LBNL1243E, Lawrence Berkeley National Laboratory, California
 Riaz and Meiburg [2003] Riaz A, Meiburg E (2003) Threedimensional miscible displacement simulations in homogeneous porous media with gravity override. Journal of Fluid Mechanics 494:95–117
 Riaz et al [2006] Riaz A, Hesse M, Tchelepi HA, Orr FM (2006) Onset of convection in a gravitationally unstable diffusive boundary layer in porous media. Journal of Fluid Mechanics 548:87–111
 Ruckenstein [1968] Ruckenstein E (1968) A generalized penetration theory for unsteady convective mass transfer. Chemical Engineering Science 23(4):363 – 371
 Straughan [2008] Straughan B (2008) Convection in porous media. In: Stability and Wave Motion in Porous Media, Applied Mathematical Sciences, vol 165, Springer New York, pp 1–45
 Sweldens [1997] Sweldens W (1997) The lifting scheme: A construction of second generation wavelets. SIAM J Math Anal 29(2):511–546
 Vasilyev and Kevlahan [2005] Vasilyev OV, Kevlahan NR (2005) An adaptive multilevel wavelet collocation method for elliptic problems. J Comput Phys 206:412–431
 Whitaker [1999] Whitaker S (1999) The Method of Volume Averaging. Springer
 Wirasaet and Paolucci [2005] Wirasaet D, Paolucci S (2005) Application of an adaptive wavelet method to naturalconvection flow in a differentially heated cavity. ASME Conference Proceedings 2005(47330):499–511
 Wirasaet and Paolucci [2006] Wirasaet D, Paolucci S (2006) The application of an adaptive wavelet method to the 3d naturalconvection flow in a differentially heated cavity. ASME Conference Proceedings 2006(47861):581–592
 Yang et al [2012] Yang D, Xue Z, Mathias SA (2012) Analysis of momentum transfer in a liddriven cavity containing a BrinkmanForchheimer medium. Transport in Porous Media 92:101–118