Buoyancy-driven convection in porous media plays a central role in a wide range of industrial and environmental settings and has received renewed attention recently because of the role of convection in geologic carbon dioxide (CO2) sequestration-- a viable solution to mitigate climate change by reducing the concentration of atmospheric CO2. The idea is to inject CO2 into brine saturated formations, whereby a gradual dissolution of CO2 into the underlying brine forms a mixture that is denser than either fluids. The unstable stratified flow eventually results in a Rayleigh--Benard type flow, in which the formation of sinking plumes acts to mix the CO2 more thoroughly into the aquifer and increases the security of storage. Consequently, accurate prediction and characterization of the mixing process is crucial in estimating and managing storage security against leakage risks. Computational modeling of CO2 storage in subsurface formations, however, is a complex multiscale transport process because of several competing flow paths, regimes, and displacement patterns accompanied by a series of geochemical reactions, across a hierarchy of length and time scales associated with multiphase flow in porous formations. This has motivated studies of simplified system where although various features of real formations are neglected, it provides a valuable framework to investigate the underlying key processes in detail. To this end, the present work aims to improve the understanding of buoyancy-driven convection in an idealized 2D porous layer by addressing two fundamental issues that have not been investigated in the past using multiple theoretical and high resolution numerical simulation: (i) convective mixing in a vertically-layered porous media; and (ii) convective mixing in a continuously perturbed porous media. We uncover new physics, both in the dynamics of convective flow in a layered porous media as well as natural convection in a system subjected to continuous forcing. These contributions can be used as a stepping stone for modeling geological scale systems.

Among the main contributions of this study is the finding that, when the porous medium is vertically-layered, thick permeability layers enhance instability compared to thin layers when heterogeneity is increased. In contrast, for thin layers the instability is weakened progressively with increasing heterogeneity to the extent that the corresponding homogeneous case, with the same density contrast, is more unstable. A resonant amplification of instability is observed within the linear regime when the dominant perturbation mode is equal to half the wavenumber of permeability variation. A weaker resonance also occurs when the dominant perturbation mode of the heterogeneous system coincides with the corresponding homogeneous system. On the other hand, substantial damping occurs when the perturbation mode is equal to the harmonic and sub-harmonic components of the permeability wavenumber. The phenomenon of such harmonic interactions influences both the onset of instability as well as the onset of convection.

Of particular physical importance is a multimodal horizontal perturbation structures, in contrast to the situation for vertical permeability variation. As a consequence, the standard eigenvalue analysis can not be used.

In the case of a continuously perturbed porous system, perturbations that are required to induce convection are introduced in the form a spatial variation of porosity in the system, a feature reflecting realistic geological settings. This form of perturbation results in an unconditionally unstable system for which the prescription of initial perturbation time and shape function are not needed. This is in contrast to a system which is perturbed in the conventional manner by introducing disturbances in the initial concentration. Using a reduced nonlinear method, the effect of harmonic variations of porosity in the transverse and streamwise direction on the onset time of convection and late time behavior is examined. It was found that the choice of perturbation method has a noticeable effect on the onset of convection and the subsequent nonlinear regime, in that the onset time of convection is reached more quickly in an impulsively perturbed system. Subsequently, an optimization procedure based on a Lagrange multiplier technique are utilized to find the optimal porosity structure that leads to the earliest onset time of convection. Scaling relationships for the optimal onset of convection and wavenumber are developed in terms of aquifer properties and initial perturbation magnitude.