Stability of miscible Rayleigh-Taylor fingers in porous media with non-monotonic density profiles
Satyajit Pramanik Manoranjan Mishra
physics.flu-dyn















[]satyajit.math16@gmail.com
Nordita, Royal Institute of Technology and Stockholm University, Stockholm, Sweden
Department of Mathematics, Indian Institute of Technology Ropar, Punjab, India




[]manoranjan.mishra@gmail.com
Department of Mathematics, Indian Institute of Technology Ropar, Punjab, India
Department of Chemical Engineering, Indian Institute of Technology Ropar, Punjab, India











We study miscible Rayleigh-Taylor (RT) fingering instability in two-dimensional homogeneous porous media, in which the fluid density varies non-monotonically as a function of the solute concentration such that the maximum density lies in the interior of the domain, thus creating a stable and an unstable density gradients that evolve in time. With the help of linear stability analysis (LSA) as well as nonlinear simulations the effects of density gradients on the stability of RT fingers are investigated. As diffusion relaxes the concentration gradient, a non-monotonic density profile emerges in time. Our simple mathematical treatment addresses the importance of density gradients on the stability of miscible Rayleigh-Taylor fingering. In this process we identify that RT fingering instabilities are better understood combining Rayleigh number (Ra) with the density gradients–hence defining a new dimensionless group, gradient Rayleigh number (Ra_g). We show that controlling the stable and unstable density gradients, miscible Rayleigh-Taylor fingers are controllable. 



64.70.Dv, 68.15.+e, 68.45.Gd, 82.65.Dp





Stability of miscible Rayleigh-Taylor fingers in porous media with non-monotonic density profiles
    Manoranjan Mishra
    December 30, 2023
=================================================================================================





§ INTRODUCTION 

There are several situations in which hydrodynamic instabilities occur due to an adverse density difference in the gravity field when a heavy fluid overlies a light fluid. The study of Rayleigh-Taylor (RT) fingering in immiscible fluids dates back to the seminal work by Rayleigh <cit.>. In RT, a convective motion generates at the interface between two fluids. Understanding RT fingering is important in various geological and industrial applications such as carbon-dioxide (CO_2) sequestration, planetary dynamics and formation of stars <cit.>. The possibility to store anthropogenic CO_2 in underground reservoir has renewed the interest to understand miscible RT fingers in porous media <cit.>. 

Notwithstanding the intensive theoretical, numerical and experimental researches of this classical hydrodynamic instability in porous media, several questioned remained unanswered. Recently, Gopalakrishnan et al. <cit.> have experimentally and numerically verified the relative role of convective and diffusive mixing in the miscible RT fingering in porous media. Through systematic analysis these authors concluded that in the fingered zone both diffusion and convection have equal contributions. The other important aspect that has not been understood yet is the role of the “transient density gradients" on the miscible RT fingering. In miscible fluids, relaxation of the concentration gradient changes the density gradient within the diffusive layer, unlike its immiscible counterparts. 

Despite the fact that the density gradient continuously changes, viscous RT fingering in miscible fluids in a porous medium is generally characterized in terms of the Rayleigh-Darcy number <cit.>, 


    Ra = κ g l_c Δρ/ϕμ D ,

where Δρ is the density difference between the two fluids, κ is the permeability of the medium, μ is the dynamic viscosity of the fluid, D is the diffusivity of the solute, l_c is a characteristic length scale, ϕ is the fluid volume fraction (porosity), and g is the gravitational acceleration. For a steady-state solution of the base solute profile, the critical Rayleigh-Darcy number for the occurrence of the instability is Ra_c = 4 π^2 <cit.>, which is violated while considering a variation in the dynamic viscosity of the fluid <cit.>. 

The scenario is completely different when the base-state solution of the solute profile is time dependent–the critical Rayleigh-Darcy number depends on the onset of instability t_c, i.e., Ra_c(t_c). On the other hand, in reactive systems a non-monotonic density profile develops in time <cit.>. Similar non-monotonic profile is also evident in double diffusive systems <cit.> or as an influence of non-ideal mixing. Water density maximizes at temperature ≈ 4^∘C. This results penetrative convection in thermal Rayleigh-Bénard convection <cit.>. It is evident from equation (<ref>) that Ra considers the density difference, not its gradient. Therefore, we define gradient Rayleigh-Darcy number  


    Ra_g(t) = κ g l_c^2 [ ∂ρ (t)/∂ x ] /ϕμ D .

Note that for thermal Rayleigh-Bénard convection, the temperature field controlling the fluid density has a constant gradient. In that case the partial derivative of the local density gradient becomes an ordinary derivative and we can write dρ/ d x  = Δρ/ l_c for l_c equal to the domain depth, and hence, as described in equation (<ref>) the gradient Rayleigh-Darcy number (Ra_g) equals the Rayleigh-Darcy number (Ra). 

In this paper we show that the miscible RT fingers are better explained in terms of local density gradients, which is best described through Ra_g(t), not Ra. We model the fluid density as piecewise smooth functions of a solute concentration and perform a linear stability analysis based on the widely used method of quasi-steady state approximation (QSSA) <cit.> and nonlinear simulations of the full problem. Numerical simulations for both monotonic and non-monotonic density profiles confirm linear stability predictions. The remainder of the paper is organized as follows. In Section <ref>, we briefly discuss a general mathematical formulation for the macroscopic description of the fingering instabilities in porous media. Results from QSSA and nonlinear simulations are discussed in section <ref>. Finally, discussion and conclusions are presented in <ref>. 



§ MATHEMATICAL MODEL AND STABILITY ANALYSIS 
 

Consider the fluid motion in a two-dimensional, vertical, homogeneous porous medium. The dynamics viscosity of the fluid μ is a constant, while the density (ρ̃) depends on a solute concentration (C) such that ρ̃(C = C_0) = ρ̃_1 and ρ̃(C = 0) = ρ̃_2, which correspond to the density of the penetrating and defending fluids, respectively. Furthermore, we assume that these two fluids are initially separated by a sharp interface (C = C_0/2) and it expands in time, as concentration diffuses. 



 §.§ Equations of motion 
 

A general macroscopic description of single species miscible RT fingering instabilities in two dimensional porous media is given by the following coupled nonlinear partial differential equations: <cit.> 



    ∇_X ·U( X ) = 0, 
       ∇_X  P = -μ/κU( X ) + ρ̃ ( X, τ ) g, 
       ϕ∂_τ C(X, τ ) + ∇_X ·[ U( X ) C( X, τ ) ] 
                         = ∇_X ·[ ϕ D ∇_X  C( X, τ ) ],


in ( X, τ ) ∈Ω× (0, ∞) wherein ∇_X  = ( ∂/∂ X , ∂/∂ Y ). We assume that the gravity acts in the positive X-direction. Here U is the two-dimensional Darcy velocity vector, p is the fluid pressure, μ is the dynamic viscosity and ρ̃ is the density of the fluid, ϕ and κ are the porosity and permeability of the medium, and D is a constant dispersion coefficient of the solute C. 



 §.§ Boundary and initial conditions 
 

The domain of definition of our problem is Ω = [-H/2, H/2] × (-∞, ∞), where H is the depth of the porous media. This system of equations are associated with the following boundary and initial conditions. Here, we are interested in time scales smaller than the time at which the equilibrium state is obtained for the species. This allows us to consider the cases for which the species concentration diffuses very slowly so that the diffusive fronts do not reach the boundaries. Therefore, we have the following boundary and initial conditions. At the inlet boundary, we have 


    U = (0, 0),  ∂_X C = 0,     [X = -H/2].

At the outlet boundary, we have 


    U = (0, 0),  ∂_X C = 0,     [X = H/2].

We assume periodic boundary conditions in the transverse direction. The initial condition is zero velocity and a step-like profile for the solute concentration. That is, 



    U = (0, 0), 
        C(X, τ = 0) = C_0(X) =
    {[  C_0, X ≤ 0;    0, X > 0;       ]., ∀ Y.





 §.§ Non-dimensionalization 
 

This system inherits (a) a diffusive time scale, τ_d = D/U^2, and (b) a convective time scale, τ_a = l_a/U, associated to a convective length scale l_a. Here, U = Δρκ g/μ is the buoyancy-induced convective velocity, where Δρ = ρ̃_1 - ρ̃_2. Evidently, diffusion occurs on a slow time scale as compared to the convection, i.e., τ_d≪τ_a. We are interested in the stability of the diffusive front to infinitesimal perturbations. We render the above system of equations dimensionless with the following scaling <cit.>



    u = U/Δρκ g/μ,  ρ = ρ̃/Δρ,   c = C/C_0, 
       x = X/ϕ D/U  = X/ϕ D μ/Δρκ g , 
        t = τ/ϕ^2 D/U^2  = τ/ϕ^2 D μ^2/( Δρκ g )^2 , 
        p =  P /ϕμ D/κ - ρ̃_1/Δρx,


from which we obtain 


    ∇·u = 0, 
    
       ∇ p = - [ u - ρ(c) e_x ], 
       ∂_t c + u·∇ c = ∇^2 c, 
       u = (0, 0),  ∂_x c = 0,     [x = - Ra/2], 
       u = (0, 0),  ∂_x c = 0,     [x =  Ra/2], 
       u = (0, 0),   c(x, t = 0) = {[    1, x ≤ 0;    0, x > 0;       ]., ∀ y,


wherein e_x is the unit vector in the x-direction, ∇ = ( ∂/∂ x , ∂/∂ y ), and Ra = Δρκ g H /ϕμ D is the Rayleigh number. Note that the boundary conditions in the y-direction are periodic. 



 §.§ Linear stability analysis 
 

As mentioned above, we are interested in the stability analysis of a time-dependent base flow. Following earlier studies, we assume no convection (pure diffusion of c) as the base flow <cit.>. Subject to this condition, as described by equation (<ref>), the mean squared displacement of the diffusive front of the species concentration grows as 4 t in time, t. For Ra ≫√(4 t), the base flow can be expressed as 



    u _b = 0, 
        c_b(x, t) =  1 / 2 erfc( x / 2 √(t)),


wherein `erfc' represents the complementary error function. Diffusive characteristic length and time scales (introduced in equations (<ref>) and (<ref>)) allow to invoke QSSA based <cit.> linear stability analyses (LSA) of the miscible RT fingers in porous media. 

QSSA method is classical in the stability analysis of time-dependent base flow, but in order to make our treatment reasonably self-contained we outline the main steps in the development of the key equations of the normal mode analysis. We introduce an infinitesimal perturbation to the governing equations (<ref>) through (<ref>) and linearized and sought solution of the following form for the perturbations in the concentration and velocity: 


    [ c̃(x, y, t), ũ(x, y, t) ] = [ c̅(x), u̅(x) ] e^i k y + ω (t_0) t,

where ω(t_0) corresponds the growth rate of the perturbations to the base state frozen at time t_0, and k is the wavenumber of the perturbations. Eliminating the pressure (p) and the transverse velocity component (v), we obtain an eigenvalue problem <cit.> 



    (  d^2/ d x^2 - ω(t_0) - k^2 ) c̅ =  d c_b(t_0)/ d xu̅, 
       (  d^2/ d x^2 - k^2 ) u̅ = -k^2 (  dρ/ d c)_c_bc̅.


For a finite t_0 and k, the eigenvalue problem, described by Eqs. (<ref>) and (<ref>), is solved numerically using a finite difference method on a uniform mesh <cit.>. In section <ref>, we present numerically computed dispersion curves ω(t_0) = ω(k; t_0) and the eigenfunction of the most unstable modes. For a detail discussion of QSSA method, interested readers are redirected to Refs. <cit.>. 



 §.§ Nonlinear simulations 
 

We use Fourier pseudo-spectral method to numerically solve the nonlinear partial differential equations <cit.>



    ∇^2 ψ = ρ^'(c) ∂ c /∂ y , 
       ∂ c /∂ t  + ∂ψ/∂ y ∂ c /∂ x  - ∂ψ/∂ x ∂ c /∂ y  = ∇^2 c,


in a doubly periodic domain of size [Ra, Ra/A], where A is the aspect ratio of the computational domain. Here ψ is the stream function, which is defined as u = ( ∂ψ/∂ y , -∂ψ/∂ x ). We consider the Fourier expansion of the unknown variables and the nonlinear terms 


    M(x, y, t) = ρ^'(c) ∂ c/∂ y, 
        N(x, y, t) = ∂ψ/∂ y ∂ c /∂ x  - ∂ψ/∂ x ∂ c /∂ y ,


as follows: 



    c(x, y, t) = ∑_j, lĉ_j, l (t) e^i(k_j x + k_l y), 
       ψ(x, y, t) = ∑_j, lψ̂_j, l (t) e^i(k_j x + k_l y), 
        M(x, y, t) = ∑_j, lM̂_j, l (t) e^i(k_j x + k_l y), 
        N(x, y, t) = ∑_j, lN̂_j, l (t) e^i(k_j x + k_l y),

 

In the Fourier space, the resulting differential-algebraic equations 



    dĉ_j, l(t)/ d t + N̂_j, l(t) = -(k_j^2 + k_l^2) ĉ_j, l(t), 
       ψ̂_j, l(t) = -M̂_j, l(t)/(k_j^2 + k_l^2),


are solved using an operator splitting predictor-corrector method. Adams-Bashforth method is used as a predictor followed by a trapezoidal method as the corrector <cit.>. 

Note that since Ra appears in the boundary conditions, the (in)stability of the system (i.e., the onset of instability, the unstable wave numbers, etc.) is independent of the choice of the value of Ra. It is verified that the numerical eigenvalues of the linear stability problem is independent for Ra ⩾ 200. Therefore, unless mentioned otherwise, we choose Ra = 200 in our linear stability analysis. On the other hand, for nonlinear simulations, the value of Ra imposes a cutoff time of our simulations–owing to the periodicity, simulations are terminated as soon as a finger reaches either of the two horizontal boundaries. The nonlinear simulation results shown in this paper are obtained using Ra = 4096 and A = 4. Grid independence has been verified and we use 1024 × 256 spectral points and Δ t = 0.2 that assure an 𝒪(10^-4) accuracy of our numerical results. 





 §.§ Density profile 


To close the mathematical description of the problem, we need to define the density profile. In various reaction-diffusion and double diffusion models, a non-monotonicity in the concentration-dependent density profile emerges due to chemical reaction and differential diffusion <cit.>. Chemical reaction induces asymmetry with regard to the initial contact line to an otherwise symmetric Rayleigh-Taylor or double-diffusive patterns <cit.>. We are interested to understand the basic fluid mechanical processes to understand the asymmetric finger patterns free from the ravages of such complicated situations of chemical reactions and differential diffusion. In this direction, we consider a phenomenological model for the density profile as a non-monotonic function of c (= C/C_0): 


    ρ̃(c) = ρ̃_1 [ 1 +  1 - ρ_m / c_m^2  c (c - 2 c_m) ],

such that, ρ̃_2 = ρ̃_1 [ 1 +  1 - ρ_m / c_m^2  (1 - 2 c_m) ]. Wherein ρ_m = [ max_C {ρ̃(C) : 0 ≤ C ≤ C_0 }] / ρ̃_1, and c_m is the corresponding dimensionless value of the concentration, i.e., ρ_m = ρ̃(C/C_0 = c_m) /ρ̃_1. In figure <ref>(a), we plot the spatial variation of the non-monotonic density profile defined in equation (<ref>) for ρ_m = 1.8,   t_0 = 25, and Ra = 200 along with the base-state concentration, c_b(x) in figure <ref>(b). As seen from figure <ref>(a), our model non-monotonic density profiles are qualitatively similar to the true density profiles obtained in various reaction-diffusion and double diffusive systems <cit.>. 



§ RESULTS AND DISCUSSION 






 §.§ Dispersion curves 
 

Linear stability analysis is performed for the density profiles shown in figure <ref>(a) and the corresponding dispersion curves are shown in figure <ref>. 
For the non-monotonic density profiles ρ̃(c), ρ_m (alternatively, Δρ = ρ_m - 1) remains unchanged, but the corresponding value of the concentration c_m changes that in turn changes the density gradient in both the stable and unstable zones. As c_m decreases, the location of the density maximum appears more towards the initial interface (x = 0)–the unstable density stratification reduces, while the stable density stratification increases. Note that the density profiles obtained corresponding to c_m = 1 (c_m = 0.5) is monotonic (symmetric about x = 0). Whereas, reducing c_m from 1^- through 0.5^+ we obtain a family of asymmetric non-monotonic density profiles. Here, we choose c_m = 0.9, 0.7 and 0.5. The corresponding dispersion curves shown in figure <ref>(c) depict that the least unstable system is obtained for c_m = 0.5–the symmetric density profile. The most unstable wavenumber and the corresponding growth rate increases as c_m increases. In response, in the nonlinear regime we anticipate long wave, rounded convective fingers to form at a later time for c_m = 0.5, whereas short wave, sharp convective fingers for c_m = 0.9, and it is evident from our numerical simulation results shown in figure <ref>. The qualitative response of the dispersion curves to the variation of the density gradients are in accordance to that discussed above. By parity of reasoning we can argue that the instability is further delayed as c_m is further decreased within the interval c_m ∈ (0, 0.5), most probably a linearly stable state is obtained as c_m → 0^+. It is noteworthy that c_m < 1/2 imposes the following restriction on the choice of ρ_m: 


    ρ_m < 1 +  c_m^2 / 1 - 2 c_m .

This is essential to insure that the dimensional density is always positive and physically realistic. 





 §.§ Eigenfunctions 
 

In figure <ref>, we plot the contours of the eigenfunction corresponding to the most unstable wavenumber estimated from the dispersion curves shown in figure <ref>. As the unstable density gradient increases (i.e., c_m increases from 0.5 through 0.9 in figures <ref>(a) through <ref>(c)), the eigenfunctions are more localized near the initial position of the interface x = 0. Moreover, as the stable density gradient increases and the unstable density gradient decreases, the eigenfunction becomes more asymmetric about the initial interface. 





 §.§ Nonlinear simulations 
 

Our linear stability results are supported with nonlinear simulations. When the density maximum occurs at the interface between two fluids, the density profile is symmetric about the interface–the density gradient at the same distance on either side of the interface is of the same magnitude, but of different sign. Although the stable density gradient is of equal strength to that of the unstable density gradient, it does not completely suppress the instability. Indeed, one anticipates long wave finger formation at a later time. Physically speaking, as the unstable density gradient moves towards a smaller value–the density maximum occurs in the defending fluid–convective motion sets in within the diffusive boundary below the density maximum and convection pushes the boundary of the density maximum slightly upwards against the resistance from the stable density gradient. This in turn reduces the strength of the stable density gradient and hence the fingers penetrate into the displacing fluid. Figure <ref> shows the concentration maps for c_m = 0.5, 0.7 and 0.9, which depict the behavior explained above. 

Therefore, the Rayleigh-Taylor instability in a miscible porous media flow is better described when Ra is combined with the density gradient, ∂ρ/∂ x = ρ^'(c) ∂ c /∂ x, alternatively in terms of Ra_g as defined in equation (<ref>). Since ρ^'(c) is fixed for a given system, from the linear theory we conclude that instability depends on ∂ c /∂ x, which varies continuously with t_0. It is observed that an increasing magnitude of the stable density gradient reduces the growth rate. As a resultant, the fluid mixing is mitigated, which we quantify in terms of the variance <cit.> 


    σ^2(t) = ⟨ c^2 ⟩ - ⟨ c ⟩^2,

where ⟨·⟩ represents an average over the spatial domain. Figure <ref> depicts that the variance decreases faster when the density maximum appears in the penetrating fluid, representing an enhanced mixing between the defending and penetrating fluids as the convective instability increases. The temporal variation of ⟨σ^2 ⟩ corresponding to a buoyantly unstable case deviates from that of a stable flow when the convection starts. This demarcates the onset of instability. It is observed from figure <ref> that for c_m = 0.5, the convective instability is yet to set in at t ≈ 2.5 × 10^4, wherein convection starts as early as t ≈ 5 × 10^3 for c_m = 0.9. We verified that when density varies linearly in c, instability sets in at an earlier time than the non-monotonic profiles defined in equation (<ref>). 





§ DISCUSSION AND CONCLUSIONS 
  

Motivated by the basic questions in porous media convection, we developed a phenomenological theory to treat density gradients in the context of the dynamics of convective instability in porous media flows. The principal phenomenon capturing our interest is the influence of the local density gradients on the convective instability. By parity of reasoning as the location of the density maximum appears above (below) the fluid interface the instability enhances (diminishes) as compared to the situation when the density maximum appears at the fluid interface. In common practice, the stability characteristics of the miscible RT fingers in porous media are described in terms of the dimensionless group Rayleigh-Darcy number (Ra) defined in equation (<ref>). For a steady-state base flow, convection starts for Ra above a critical value, Ra > Ra_c. For a time-dependent base flow, the concept of critical Rayleigh-Darcy number is replaced by the critical time of instability (onset of instability) t_c <cit.>. Surprisingly, for the same Ra, one obtains different t_c depending on the density-concentration relation, which may lead to different interpretations of the miscible RT fingers. 

Here, we present a unified theory in terms of the density gradients. Through a systematic analysis of a series of density profiles, we are able to identify the effects of density gradients on miscible Rayleigh-Taylor fingering. Our density profiles can be uniquely defined as 


    ρ̃(c) = ρ̃_1 + Δρ f(c),

wherein 


    f(c) = -c(c-2c_m),    Δρ = ρ̃_1 (ρ_m - 1)/c_m^2,

for the non-monotonic density profiles (<ref>). The competition between a locally stable and a locally unstable density gradients are discussed through various non-monotonic profiles. Our model construction allows to single out the effects of variation in stable as well as unstable density gradients. Interestingly, the explanation of miscible Rayleigh-Taylor fingers through the density gradients presented here has important implications for geological CO_2 sequestration and other porous media convection driven out of chemical reactions <cit.>, where the fluid density varies non-monotonically on the solute concentration. One of the most important results of our theory is the following. The role of density gradients can give a possible explanation of the asymmetric RT finger formation in double diffusive reactive system <cit.> without solving the multi-component reaction-diffusion-convection equations.  

As shown in figure <ref>(a), the unstable density gradient, [ ∂ρ/∂ x ]_x > x_m = [ f^'(c) ∂ c /∂ x ]_x > x_m, decreases as the density maximum appears higher above the initial fluid interface, which we mark at c = 0.5 (i.e., x = 0, which is obtained by parity of symmetry of the error function profile). Here, x_m is obtained inverting the relation c_m = c(x_m, t_0), i.e., x_m depends on t_0. Additionally, the magnitude of the stable density gradient, [ ∂ρ/∂ x ]_x < x_m = [ f^'(c) ∂ c /∂ x ]_x < x_m, increases as c_m decreases. Note that although Ra = 200 remains unchanged for the three density profiles considered (figure <ref>), the corresponding dispersion curves are different (figure <ref>). This signifies that if the density profile of a physical system evolves non-monotonically having different local density gradients without affecting Ra, one may over/under predict the stability. Therefore, an a priori estimate about the stability of a system for a given Ra may be misleading. Contrary to that, a quantification of the stability of a system in terms of the local density gradient as discussed above improves our understanding. This confirms that when the solute concentration field possesses transient nature, the Rayleigh-Taylor instability in porous media is better explained in terms of the local density gradients, but not the difference between the maximum and minimum density. 

Recently, Toppaladoddi and Wettlaufer <cit.> studied high Rayleigh numbers penetrative convection of a fluid confined between two plates. One of the major results of their work was to identify the essence of a dimensionless parameter (Λ)–the ratio of the top and bottom plate temperatures relative to the temperature that maximizes the water density–apart from Rayleigh and Prandtl numbers. In particular, it was reported that the the penetration of the plumes generated from the hot bottom plate is hindered by a strong stable layer stratification for a large Λ. Similarly, our results corresponding to the density profile ρ̃_ nm, 3 can be explained in terms of the dimensionless parameter 

    Γ = 1-c_m/c_m > 0.

It is clear from figures <ref>(a-c), which depict time evolution of the concentration field for Γ = 1, 3/7, 1/9, respectively, that due to strong stable layer stratifications the fingers are impeded to penetrate into the stable layer as Γ increases, and the flow remains stable for a longer period. As Γ decreases, the fingers penetrate the stable upper layer and the patterns approach toward an up-down symmetric situation. Therefore, it can be expected that the standard Rayleigh-Taylor fingers with a linear density-concentration relation is neared as Γ→ 0. 

However, we must note that although the qualitative dependences were understood in the course of this study, the quantitative relation between the growth rate and the density gradients remained unclear. Note that most of the recent studies on stability analysis of miscible fingering instabilities have emphasized on various other stability methods–QSSA with respect to a suitably defined self-similar variable, ξ (SS-QSSA) <cit.>, initial value calculation (IVC) <cit.>, non-modal stability analysis (NMA) <cit.>, etc.–than the classical QSSA. Recall, the instantaneous growth rates obtained using SS-QSSA and QSSA are related via <cit.>


    ω_x = ω_ξ +  1 / 4 t_0 ,

where t_0 is the frozen diffusive time <cit.>. The present study focuses upon the role of local density gradients on Rayleigh-Taylor fingers in a homogeneous porous medium. A qualitative understanding without any direct comparison of with experiments allows to choose QSSA method for the stability analysis. We believe the qualitative results presented here remains unaffected by the choice of the method. Interested researchers are encouraged to choose the most appropriate stability method for a more quantitative measurements, in particular while comparing with the experiments. 

We show that the onset of instability, and hence the fluid mixing can be significantly altered by varying the stable and unstable density gradients parameterized in terms of Γ. This indicates that chemical reactions, that result a continuous evolution of the density gradients, control convective instability and fluid mixing. Mathematically, convective instability occur for any density profile as long as density decreases along the direction of gravity–the onset of instability depends on the gradients of the density in the stable as well as unstable zones. Diffusion causes to change the density gradient in both the stable and unstable zones. Convection sets in early for a steep (shallow) unstable (stable) density gradient and the onset time  increases as the unstable (stable) density gradient becomes shallow (steep). Present analysis can be extended to understand miscible viscous fingering with non-monotonic viscosity profiles <cit.>. A more complicated situations of non-monotonic density effects on the convective instability in porous media awaits understanding when the solute concentration is localized within a finite region <cit.>. 


SP thanks Prof. John Wettlaufer, Cristobal Arratia, and Francesco Macarella for many helpful discussions. This work started when SP was at the Department of Mathematics, IIT Ropar and was generously funded a Ph.D. fellowship by the National Board for Higher Mathematics (NBHM), Department of Atomic Energy (DAE), Govt. of India. 






99

Rayleigh1882
Rayleigh. 1882  Investigation of the Character of the Equilibrium of an
  Incompressible Heavy Fluid of Variable Density. Proceedings of the
  London Mathematical Society 1-14, 170–177.

Gopal2017
Gopalakrishnan SS, Carballido-Landeira J, De Wit A, Knaepen B. 2017  Relative
  role of convective and diffusive mixing in the miscible Rayleigh-Taylor
  instability in porous media. Phys. Rev. Fluids 2, 012501.

Huppert2014
Huppert HE, Neufeld JA. 2014  The Fluid Mechanics of Carbon Dioxide
  Sequestration. Annu. Rev. Fluid Mech. 46, 255–272.

Slim2014
Slim AC. 2014  Solutal-convection regimes in a two-dimensional porous medium.
  J. Fluid Mech. 741, 461–491.

Nakanishi2016
Nakanishi Y, Hyodo A, Wang L, Suekane T. 2016  Experimental study of 3D
  Rayleigh-Taylor convection between miscible fluids in a porous medium. Advances in Water Resources 97, 224 – 232.

Binda2017
Binda L, Hasi CE, Zalts A, D'Onofrio A. 2017  Experimental analysis of density
  fingering instability modified by precipitation. Chaos: An
  Interdisciplinary Journal of Nonlinear Science 27, 053111.

Riaz2006
Riaz A, Hesse M, Tchelepi HA, Orr FM. 2006  Onset of convection in a
  gravitationally unstable diffusive boundary layer in porous media. J.
  Fluid Mech. 548, 87–111.

Horton1945
Horton CW, Jr. FTR. 1945  Convection Currents in a Porous Medium. Journal
  of Applied Physics 16, 367–370.

Sabet2017b
Sabet N, Hassanzadeh H, Abedi J. 2017  A New Insight into the Stability of
  Variable Viscosity Diffusive Boundary Layers in Porous Media under Gravity
  Field. AIChE Journal.

Lemaigre2013
Lemaigre L, Budroni MA, Riolfo LA, Grosfils P, Wit AD. 2013  Asymmetric
  Rayleigh-Taylor and double-diffusive fingers in reactive systems. Physics of Fluids 25, 014103.

Loodts2014
Loodts V, Thomas C, Rongy L, De Wit A. 2014  Control of Convective Dissolution
  by Chemical Reactions: General Classification and Application to
  CO_2 Dissolution in Reactive Aqueous Solutions. Phys.
  Rev. Lett. 113, 114501.

Loodts2016
Loodts V, Trevelyan PMJ, Rongy L, De Wit A. 2016  Density profiles around
  A+BC reaction-diffusion fronts in partially
  miscible systems: A general classification. Phys. Rev. E 94,
  043115.

Trevelyan2011
Trevelyan PMJ, Almarcha C, De Wit A. 2011  Buoyancy-driven instabilities of
  miscible two-layer stratifications in porous media and Hele-Shaw cells. J. Fluid Mech. 670, 38–65.

Srikanth2017
Toppaladoddi S, Wettlaufer JS. 2017  Penetrative Convection at High
  Rayleigh Numbers. ArXiv e-prints.

Tan1986
Tan CT, Homsy GM. 1986  Stability of miscible displacements in porous media:
  Rectilinear flow. Phys. Fluids (1958-1988) 29, 3549–3556.

Homsy1987
Homsy GM. 1987  Viscous Fingering in Porous Media. Annu. Rev. Fluid Mech.
  19, 271–311.

Pramanik2015e
Pramanik S, Wit AD, Mishra M. 2015a  Viscous fingering and deformation of a
  miscible circular blob in a rectilinear displacement in porous media. J.
  Fluid Mech. 782, R2.

Pramanik2015d
Pramanik S, Hota TK, Mishra M. 2015b  Influence of viscosity contrast on
  buoyantly unstable miscible fluids in porous media. J. Fluid Mech.
  780, 388–406.

Pramanik2016a
Pramanik S, Mishra M. 2016  Coupled effect of viscosity and density gradients
  on fingering instabilities of a miscible slice in porous media. Phys.
  Fluids 28, 084104.

Manickam1995
Manickam O, Homsy GM. 1995  Fingering instabilities in vertical miscible
  displacement flows in porous media. J. Fluid Mech. 288,
  75–102.

Pramanik2013
Pramanik S, Mishra M. 2013  Linear stability analysis of Korteweg stresses
  effect on miscible viscous fingering in porous media. Phys. Fluids
  (1994-present) 25, 074104.

Pramanik2015f
Pramanik S. 2015 Analysis of hydrodynamic instabilities in miscible
  displacement flows in porous media.
PhD thesis Indian Institute of Technology Ropar.

Tan1988
Tan CT, Homsy GM. 1988  Simulation of nonlinear viscous fingering in miscible
  displacement. Phys. Fluids (1958-1988) 31, 1330–1338.

Loodts2015
Loodts V, Rongy L, De Wit A. 2015  Chemical control of dissolution-driven
  convection in partially miscible systems: theoretical classification. Phys. Chem. Chem. Phys. 17, 29814–29823.

Pramanik2015c
Pramanik S, Mishra M. 2015  Viscosity scaling of fingering instability in
  finite slices with Korteweg stress. EPL (Europhysics Letters)
  109, 64001.

DeWit2016
De Wit A. 2016  Chemo-hydrodynamic patterns in porous media. Philosophical
  Transactions of the Royal Society of London A: Mathematical, Physical and
  Engineering Sciences 374.

Hota2015a
Hota TK, Pramanik S, Mishra M. 2015a  Onset of fingering instability in a
  finite slice of adsorbed solute. Phys. Rev. E 92, 023013.

Hota2015b
Hota TK, Pramanik S, Mishra M. 2015b  Nonmodal linear stability analysis of
  miscible viscous fingering in porous media. Phys. Rev. E 92,
  053007.

Kim2011
Kim MC, Choi CK. 2011  The stability of miscible displacement in porous media:
  Nonmonotonic viscosity profiles. Phys. Fluids (1994-present)
  23, 084105.

Manickam1993
Manickam O, Homsy GM. 1993  Stability of miscible displacements in porous media
  with nonmonotonic viscosity profiles. Phys. Fluids A: Fluid Dynamics
  (1989-1993) 5, 1356–1367.

Haudin2016
Haudin F, Callewaert M, De Malsche W, De Wit A. 2016  Influence of nonideal
  mixing properties on viscous fingering in micropillar array columns. Phys. Rev. Fluids 1, 074001.

DeWit2005
De Wit A, Bertho Y, Martin M. 2005  Viscous fingering of miscible slices. Phys. Fluids (1994-present) 17, 054114.

Mishra2008
Mishra M, Martin M, De Wit A. 2008  Differences in miscible viscous fingering
  of finite width slices with positive or negative log-mobility ratio. Phys. Rev. E 78, 066306.













