Stability of miscible Rayleigh-Taylor fingers in porous media with non-monotonic density profiles

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.


I. 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 [1]. 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 [2, and refs. therein]. The possibility to store anthropogenic CO 2 in underground reservoir has renewed the interest to understand miscible RT fingers in porous media [3][4][5][6].
Notwithstanding the intensive theoretical, numerical and experimental researches of this classical hydrodynamic instability in porous media, several questioned remained unanswered. Recently, Gopalakrishnan et al. [2] 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. * satyajit.math16@gmail.com † manoranjan.mishra@gmail.com 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 [7], 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 [8], which is violated while considering a variation in the dynamic viscosity of the fluid [9]. The scenario is completely different when the basestate 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 [10][11][12]. Similar non-monotonic profile is also evident in double diffusive systems [13] or as an influence of nonideal mixing. Water density maximizes at temperature ≈ 4 • C. This results penetrative convection in thermal Rayleigh-Bénard convection [14, and refs. therein]. It is evident from equation (1) that Ra considers the density difference, not its gradient. Therefore, we define gradient Rayleigh-Darcy number 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ρ dx = ∆ρ l c for l c equal to the domain depth, and hence, as described in equation (2) 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) [15] 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 II, 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 III. Finally, discussion and conclusions are presented in §IV.

II. 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.
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.

B. 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 At the outlet boundary, we have 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, 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 [18,19] from which we obtain wherein e x is the unit vector in the x-direction, ∇ = ∂ ∂x , ∂ ∂y , and Ra = ∆ρκgH φµD is the Rayleigh number.
Note that the boundary conditions in the y-direction are periodic.

D. 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 [7,18]. Subject to this condition, as described by equation (8c), the mean squared displacement of the diffusive front of the species concentration grows as 4t in time, t. For Ra √ 4t, the base flow can be expressed as wherein 'erfc' represents the complementary error function. Diffusive characteristic length and time scales (introduced in equations (7b) and (7c)) allow to invoke QSSA based [15] 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 (8a) through (8c) and linearized and sought solution of the following form for the perturbations in the concentration and velocity: 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 [20] For a finite t 0 and k, the eigenvalue problem, described by Eqs. (11a) and (11b), is solved numerically using a finite difference method on a uniform mesh [21]. In section III, 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. [15,20,22].

E. Nonlinear simulations
We use Fourier pseudo-spectral method to numerically solve the nonlinear partial differential equations [18] as follows: In the Fourier space, the resulting differential-algebraic equations 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 [see 18,19,22,23, and refs. therein for the details of the numerical method]. 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 O(10 −4 ) accuracy of our numerical results.

F. 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 [10, 12, 13, 24, and refs. therein]. Chemical reaction induces asymmetry with regard to the initial contact line to an otherwise symmetric Rayleigh-Taylor or double-diffusive patterns [10]. 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 ): such that, i.e., ρ m =ρ (C/C 0 = c m ) ρ 1 . In figure 1(a), we plot the spatial variation of the non-monotonic density profile defined in equation (16) for ρ m = 1.8, t 0 = 25, and Ra = 200 along with the base-state concentration, c b (x) in figure 1(b). As seen from figure 1(a), our model nonmonotonic density profiles are qualitatively similar to the true density profiles obtained in various reaction-diffusion and double diffusive systems [10,12,13,24].

A. Dispersion curves
Linear stability analysis is performed for the density profiles shown in figure 1(a) and the corresponding dispersion curves are shown in figure 2. 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 2(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 4. 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 : This is essential to insure that the dimensional density is always positive and physically realistic.

B. Eigenfunctions
In figure 3, we plot the contours of the eigenfunction corresponding to the most unstable wavenumber estimated from the dispersion curves shown in figure 2. As the unstable density gradient increases (i.e., c m increases from 0.5 through 0.9 in figures 3(a) through 3(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.

C. 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 4 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 , alterna- tively in terms of Ra g as defined in equation (2). 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 [25] where · represents an average over the spatial domain. Figure 5 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 5 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 (16).

IV. 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 (1). 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 [7]. 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 wherein for the non-monotonic density profiles (16). The competition between a locally stable and a locally unstable density gradients are discussed through various nonmonotonic 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 [12, 26, and refs. therein], where the fluid density varies nonmonotonically 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 [10] without solving the multicomponent reaction-diffusion-convection equations. As shown in figure 1(a), the unstable density gradient, ∂ρ ∂x x>xm = f (c) ∂c ∂x x>xm , 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<xm = f (c) ∂c ∂x x<xm , increases as c m decreases. Note that although Ra = 200 remains unchanged for the three density profiles considered (figure 1), the corresponding dispersion curves are different ( figure 2). 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 [14] 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 It is clear from figures 4(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 updown 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) [21], initial value calculation (IVC) [27], non-modal stability analysis (NMA) [28], etc.-than the classical QSSA. Recall, the instantaneous growth rates obtained using SS-QSSA and QSSA are related via [29] where t 0 is the frozen diffusive time [21]. 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 [30,31]. 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 [17,32,33].