1 Introduction↩︎

In this paper, we are interested in addressing the problem of viscous fingering in two-phase, compositional displacement process, where the invading and displaced phases are not fully miscible, but mass transfer occurs between the phases (partially miscible flow). Our work focuses on the prediction of unstable modes of perturbation using linear stability analysis.

Flow displacement processes in porous media are critical to numerous industrial and environmental applications. In enhanced oil recovery operations, the injection of water, gas, or chemical agents displaces oil from reservoir rocks, directly impacting production efficiency and economic viability. Similarly, in carbon sequestration and storage (CCS) projects, CO\(_2\) is injected into deep saline aquifers or depleted hydrocarbon reservoirs, where the stability of the displacement front determines plume migration patterns. Groundwater remediation efforts also rely on displacement processes to remove or contain contaminants in subsurface formations. In all these applications, the stability of the displacement front is essential for process efficiency, as unstable flows lead to viscous fingering—a phenomenon where the displacing fluid penetrates the displaced fluid in preferential pathways, creating finger-like intrusions. This instability significantly reduces sweep efficiency, causes premature breakthrough of the displacing fluid, and results in substantial bypassing of the target fluid, diminishing the effectiveness of enhanced oil recovery, compromising CO\(_2\) storage security, and limiting remediation success.

Unlike flow channelling, which arises from spatial heterogeneity in the porous medium, viscous fingering stems from variations in fluid properties—such as viscosity and density—across the displacement front. These property differences may result from displacement between two immiscible phases or from mixing of two fully miscible fluids [1]. When a small, localized perturbation at the displacing front (with the displacing fluid advancing across the front) leads to steeper (negative) pressure gradient, local flow rate is increased, amplifying the initial perturbation and leading to instability. If, on the other hand, the perturbation reduces the local pressure gradient, the displaced fluid resists the local advance, and the displacing front is stabilized. The primary drivers of instability are, for most scenarios, adverse mobility and density ratios. Conversely, favourable mobility ratios or density contrasts act as stabilizing forces.

Besides viscous and gravity forces, dissipative forces also influence the growth rate and length scale of finger-like perturbations. In macroscopic, immiscible flow, the primary dissipative mechanism is capillary pressure, while in miscible flow, molecular diffusion and mechanical dispersion play this role. These dissipative forces counteract any gradients in composition (whether saturation or molecular concentration), and act as stabilizing mechanisms by (\(i\)) smoothing the compositional gradient across the displacing front, limiting the growth rate of perturbations, and (\(ii\)) by smoothing transverse gradients along the front, effectively dispersing small-scale perturbations—which, in the context of modal analysis, results in a wavelength cutoff value below which flow is stabilized.

Linear stability analysis provides a powerful mathematical framework for predicting the onset of viscous fingering and characterizing the early stages of instability development. The analysis for displacement processes typically begins with a base state—the one-dimensional displacement profile at the displacement front—and superimposes small-amplitude, wave-like disturbances in the transverse direction (normal mode analysis). By linearizing the governing equations around this base state, one obtains an eigenvalue problem whose solution yields the growth rate as a function of perturbation wavelength. A positive growth rate indicates a mode is unstable, with its amplitude growing exponentially over time, while negative growth rates indicate stable modes that decay over time.

Since the pioneering study by [2], extensive research has been published on viscous fingering. A classical review on the subject was given by [1]. More recently, [3] provided an extensive review on experimental and computational studies on miscible and immiscible viscous fingering. A review focusing on variational methods for numerical simulation of unstable miscible displacement was given by [4], and a shorter, but well organized review on immiscible viscous fingering was provided by [5].

Linear stability analysis alone has generated substantial literature. Following the foundational work of [6] and [7], several works have addressed linear stability analysis for immiscible viscous fingering (see e.g., for a recent review [5]), including notable contributions by [8], [9], [10], [11], [12], and [13]. Similarly, extensive research exists for miscible displacement, with important examples including [14], [15], [16], [17], [18], and [19]. An important subset of miscible displacement is gravity-driven instability—typically the result of concentration gradients producing adverse density gradients in the aqueous phase—with notable work by [20], [21], [22], and [23].

The literature on viscous fingering in multiphase compositional displacement (i.e., partially miscible flow), however, is considerably more limited, particularly regarding linear stability analysis. Many practical applications involve displacement processes at conditions that are neither immiscible nor fully miscible, with limited mass transfer occurring between invading and displaced fluids such as CO\(_2\)-enhanced oil recovery below the miscibility pressure or carbon sequestration in saline aquifers. This underscores the importance of studying the viscous fingering in the context of multiphase compositional model.

Most studies of unstable multiphase compositional flow address fingering through numerical simulation of the nonlinear governing equations, analysing instability either by visual inspection of compositional fields or through indirect metrics such as sweep efficiency. In this category, [24] examined the impact of several factors—including capillarity, molecular diffusion, flow rates, domain geometry and medium properties—on viscous fingering, and concluded that displacement in partially miscible conditions differs profoundly from the fully miscible case. In the context of linear stability analysis, the study by [23], which examines evaporation-driven density instabilities in porous media, considers a model incorporating both immiscible displacement between two phases (air and water) and miscible transport of a solute (NaCl) within the aqueous phase. A key distinction from the present work, however, is that interphase mass transfer does not occur in their analysed case. We are aware of only one work, by [25], that accounts for mass transfer between phases. The authors focused on CO\(_2\) displacement of heavy oil under conditions that are predominantly immiscible, but with limited interphase mass transfer. Their analysis assumed that: (\(i\)) displacement occurs in an idealized step-saturation profile with constant volume fractions on either side of the displacement front; (\(ii\)) capillary forces are neglected; and (\(iii\)) mass transfer follows a prescribed function rather than being determined from phase-equilibrium relationships. Their results show that mass transfer of CO\(_2\) into the oil phase has a stabilizing effect on the displacement by reducing both the front velocity and the viscosity gradient across the front—findings that will be shown to be in agreement with our own analysis.

In the present work, we focus the analysis of linear stability to partially miscible displacement in a two-phase, two-component system with mass transfer between phases. In contrast to [25], our model accounts for: gravity forces; fractional flow and capillary forces in the two-phase region; variable viscosity and mechanical dispersion in single-phase regions. We limit our study to the common scenario where the invading fluid is gaseous and the displaced fluid is liquid. While the two-component system does not fully capture the physics of many compositional flow problems—especially when the liquid phase contains multiple volatile components—it provides a useful framework that considerably simplifies the mathematical description of the problem (see Section 2.1). By accounting for partial miscibility, we attempt to the bridge the gap between the well-studied limiting cases of immiscible and miscible displacement. To the best of our knowledge, this is the first work dealing a linear stability analysis of a genuine two-phase, two-component system taking into account interphase mass transfer.

For the idealized system where pressure variation is small compared to the system pressure and liquid density is weakly dependent on composition, our two-phase, two-component model becomes quite similar to the immiscible case solved by [12]. The key differences introduced by mass-transfer effects in the two-component case are: (\(i\)) shock calculations deviate from the classical Buckley-Leverett approach; (\(ii\)) a single-phase flow regime appears downstream of the shock region; and (\(iii\)) flow functions become discontinuous at the transition from two-phase to single-phase flow. The basic steps of the linear stability analysis for the proposed model are:

1. From the general equations for multiphase compositional flow, assuming thermodynamic equilibrium, derive the governing equations for the simpler case of a two-phase, two-component system with additional simplifying assumptions (see Section 2.1);

2. From the one-dimensional, purely advective case, obtain the shock front configuration following the procedure developed by [26];

3. Determine the steady-state composition profile around the shock (the travelling-wave solution of the one-dimensional problem)—the base state;

4. Derive the perturbation equations by decomposing the primary variables into base-state solutions and normal mode perturbations, then linearizing the flow functions around the base state—the eigenvalue problem;

5. Solve the eigenvalue problem using the matched initial value problem (MIVP) method (adapted from [11]), which consists of:

   1. Calculate the far-field solutions of the linearized flow equations for a given wavenumber to obtain the initial conditions on the far-upstream and far-downstream sides, characterized by arbitrary constants;

   2. Integrate the linearized equations from both the far-upstream and far-downstream sides toward the matching point (the flow regime transition point) using the calculated initial conditions;

   3. Calculate the coupling terms (jump conditions) for the discontinuous derivatives to arrive at the matched initial value system, which yields a non-trivial solution characterized by a real-valued growth rate.

2 Mathematical Formulation↩︎

In this section, we develop the mathematical formulation for two-dimensional flow in an isotropic porous medium with uniform properties, for a two-phase, two-component system in thermodynamic equilibrium. We identify the components in our binary mixture as “\(a\)” (the lighter component) and “\(b\)” (the heavier component).

We are interested in the case of a gas (rich in component \(a\)) displacing an undersaturated liquid (rich in component \(b\)). Flow is considered uniform in the \(x\)-direction with a shock at the interface between the invading and displaced fluids. We assume that if instability occurs, it will be at the shock where the total mobility gradient is highest—which [12] shows is a reasonable assumption. Therefore, after defining our base-state solution as the one-dimensional, travelling-wave solution around the shock, our analysis then consists in subjecting the base-state to two-dimensional, small wave-like perturbations, and checking the stability of the system.

Figure 1 shows a sketch of the system configuration. Although not explicitly shown in our diagram, gravity is accounted for in our formulations and results shown in Section 3 encompass three displacement directions: horizontal, upward and downward.

Figure 2 shows an example of a full Buckley-Leverett profile for the overall composition of component \(a\) (\(z^a\)), and the solution at a localized region around the shock, which at large displacement periods converges to our base-state solution.

In the following subsections, we first discuss the flash equations for a system of only two components, list the simplifying assumptions, write the two main flow equations, briefly discuss the calculation of the shock in the presence of mass transfer, develop the base-state equations given by the moving-frame transformation, and, finally, we derive the perturbation equations and the coupling solution around the phase-change discontinuity.

2.1 Phase equilibrium↩︎

Phase equilibrium calculations are considerably simpler for binary systems than for systems containing three or more components. To illustrate this, we first consider a mixture of \(N_c\) components with its overall composition given by \({z^\gamma=n^\gamma/n_t}\), \({\gamma=1,\ldots,N_c}\), with \(n^\gamma\) as the number of moles of component \(\gamma\) and \(n_t\) as the total number of moles in the system. Suppose that at a certain pressure and temperature (\(P,T\)) this system splits into two phases, gas and liquid, with their respective compositions given by \({y^\gamma=n^\gamma_g/n_g}\) and \({x^\gamma=n^\gamma_\ell/n_\ell}\), where \(n_\alpha\) is the number of moles of phase \(\alpha\) and \(n^\gamma_\alpha\) the number of moles of \(\gamma\) in \(\alpha\). The two phases are in thermodynamic equilibrium if the fugacity of each component is the same in both phases [27]: \[\hat{f}_g^\gamma(P, T, \{y\}) = \hat{f}_\ell^\gamma(P, T, \{x\}),\quad \text{for}\quad \gamma=1,\ldots,N_c, \label{eq:fugacity}\tag{1}\] where \(\hat{f}_\alpha^\gamma\) is the fugacity of component \(\gamma\) in phase \(\alpha\), \({\{y\}=\{y^1,\ldots,y^{N_c-1}\}}\) and \({\{x\}=\{x^1,\ldots,x^{N_c-1}\}}\) are the sets of independent values of the gas and liquid compositions (since, by definition, \({\sum{y^\gamma}=1}\) and \({\sum{x^\gamma}=1}\)). The expressions for \(\hat{f}_\alpha^\gamma\) can be written with an appropriate equation of state (EOS).

Assuming the system is in equilibrium, the composition of the gas and liquid phases can be found by a flash calculation, which consists of solving the thermodynamic equilibrium equations (1 ) together with the material balance, or tie-line, equations that link the overall composition of the mixture with the equilibrium compositions: \[z^\gamma= L_g\,y^\gamma+ (1-L_g)\,x^\gamma,\quad \text{for}\quad \gamma=1,\ldots,N_c-1, \label{eq:tieline}\tag{2}\] where \({L_g=n_g/n_t}\) is the gas mole fraction.

In total, there are \(2N_c-1\) primary unknowns and \(2N_c-1\) independent equations, so the system is perfectly determined. If \(N_c\geq3\), Equations 1 and 2 need to be solved simultaneously. However, for the binary case (\({N_c=2}\)), we have two unknowns and two equations in Equation 1 , and the equilibrium compositions can be solved independently of the overall composition—i.e. \({y^\gamma=y^\gamma(P,T)}\) and \({x^\gamma=x^\gamma(P,T)}\). Equation 2 is then decoupled from Equation 1 , and with \(z^\gamma\), \(y^\gamma\) and \(x^\gamma\) given, we can solve for \(L_g\): \[L_g = \frac{z^\gamma- x^\gamma_e}{y^\gamma_e - x^\gamma_e},\quad \text{for}\quad \gamma=a,b, \label{eq:tieline2}\tag{3}\] where \(y_e^\gamma\) and \(x_e^\gamma\) are the constant equilibrium compositions of the gas and liquid phases inside the two-phase region.

In summary, we need to solve Equation 1 only once to obtain the equilibrium compositions of the gas and liquid phases in the two-phase region, and the flash calculation is reduced to solving Equation 3 for the overall composition \(z^\gamma\). We can also easily identify the two-phase and single-phase regions. Indeed, assuming that component \(a\) is lighter than component \(b\) (i.e. \(y_e^a > x_e^a\) and \(y_e^b < x_e^b\)), [26] makes the following observations regarding Equation 3 :

1. If \(z^a > y_e^a\), then the mixture is an undersaturated (or single-phase) gas, and also \(L_g > 1\);

2. If \(z^a < x_e^a\), the mixture is an undersaturated (or single-phase) liquid, and \(L_g < 0\);

3. \(y_e^\gamma\) and \(x_e^\gamma\) are constant inside the two-phase region (\(x_e^a \leq z^a \leq y_e^a\)).

Item (3) implies that, for a fixed (\(P,T\)), the properties—such as viscosity (\(\mu_\alpha\)) and density (\(\rho_\alpha\))—for the gas and liquid phases are constant in two-phase region. In the single-phase regions, that is usually not true, and \(\mu_\alpha\) and \(\rho_\alpha\) are dependent of the overall composition. However, the density of the liquid phase may only be weakly dependent on \(z^a\), and for simplicity we assume that \(\rho_\ell=constant\).

Normally, to calculate the phase density we calculate its molar volume \(V_\alpha\) using our EOS, and then use the relation: \[\rho_\alpha= \bar{M}_\alpha c_\alpha, \label{eq:density1}\tag{4}\] where \(\bar{M}_\alpha\) is the average molar weight of phase \(\alpha\) and \(c_\alpha=1/V_\alpha\) its molar concentration (or molar density).

For the binary case (\(\gamma=a,b\)), in the pure-liquid region, we can write: \[\bar{M}_\ell = \sum_\gamma M^\gamma z^\gamma= M^b + (M^a - M^b) z^a. \label{eq:Mwmix}\tag{5}\]

With \(\rho_{\ell e}\) the value of the liquid density in the two-phase region, we can ensure that \(\rho_\ell = \rho_{\ell e}\) in the single-phase liquid region by calculating \(c_\ell\) with Equation 6 , instead of using the EOS.

\[c_\ell = \frac{\bar{M}_\ell}{\rho_{\ell e}}. \label{eq:density2}\tag{6}\]

2.2 Flow equations↩︎

For the derivation of the flow equations in the next sections, we make the following simplifying assumptions:

1. The fluid injected is in the two-phase region (\(y^a < z^a < x^a\)), and the single-phase gas region is bypassed;

2. The liquid phase density in the single-phase liquid region is independent of the liquid composition and equal to \(\rho_{\ell e}\);

3. The temperature is constant;

4. The pressure drop along flow is small compared to the initial pressure such that the fluid properties are evaluated at initial pressure;

5. The porous medium is uniform and isotropic;

6. Mechanical dispersion dominates over molecular diffusion; can be described by a Fickian-type equation; and the dispersivities are independent of the molecular concentrations (more details in Section 2.2.2);

7. The capillary pressure derivative with respect to gas saturation is a function of the phases’ relative permeabilities (see Section 2.7.3).

2.2.1 Mass balance↩︎

The material balance equation of component \(\gamma\) in phase \(\alpha\), for a multiphase, compositional flow in porous media, for a bounded domain \(\Omega \subset \mathbb{R}^n\, (n=1,2)\), without the presence of source/sink terms, is given by [28]: \[\phi \frac{\partial S_\alpha\rho_\alpha^\gamma}{\partial t} + \phi \nabla \cdot \left[S_\alpha\left(\rho_\alpha^\gamma\mathbf{V}_\alpha+ \mathbf{J}_\alpha^{\gamma}\right) \right] - \sum_{\beta\neq\alpha} f_{\beta\to\alpha}^\gamma= 0, \label{eq:matbal}\tag{7}\] where \(\phi\) is the medium porosity, \(S_{\alpha}\) the saturation of phase \(\alpha\), \(\rho_\alpha^\gamma\) the mass density of component \(\gamma\) in phase \(\alpha\), \(\mathbf{V}_\alpha\) is the phase velocity, \(\mathbf{J}_\alpha^{\gamma}\) represents the sum of diffusive and dispersive fluxes of component \(\gamma\) in phase \(\alpha\) and \(f_{\beta\to\alpha}^\gamma\) represents the transport of \(\gamma\)-mass from phase \(\beta\) to phase \(\alpha\).

Recognizing that in a two-phase system \(f_{\beta\to\alpha}^\gamma=-f_{\alpha\to\beta}^\gamma\), we arrive at the equation for the total mass balance of \(\gamma\) by summing over phases \(g\) (gas) and \(l\) (liquid) in Equation 7 : \[\phi \frac{\partial\left(S_g \rho_g^\gamma+ S_\ell \rho_\ell^\gamma\right) }{\partial t} + \nabla \cdot \left(\rho_g^\gamma\mathbf{u}_g + \rho_\ell^\gamma\mathbf{u}_\ell \right) + \phi \nabla \cdot \left[ \left( S_g \mathbf{J}_g^{\gamma} + S_\ell \mathbf{J}_\ell^{\gamma} \right) \right] = 0, \label{eq:matbal95total1}\tag{8}\] where \(\mathbf{u}_\alpha=\phi S_\alpha\mathbf{V}_\alpha\) is the specific discharge—or Darcy velocity—of phase \(\alpha\).

2.2.2 Dispersive flux↩︎

For simplicity, we assume that mechanical dispersion dominates over molecular diffusion. For the case of uniform flow in an isotropic medium, we can express the mechanical dispersion flux in a Fickian-type law as [28]: \[\mathbf{J}_{\alpha}^\gamma= -M^\gamma\mathbf{D}_{\alpha}^{\gamma\delta} \nabla c_\alpha^\gamma, \label{eq:fickdisp}\tag{9}\] where \(\mathbf{D}_{\alpha}^{\gamma\delta}\) is the dispersion tensor given by: \[\mathbf{D}_{\alpha}^{\gamma\delta} = \begin{bmatrix} a_L & 0 \\ 0 & a_T \\ \end{bmatrix} V_{x,\alpha}, \label{eq:disptensor}\tag{10}\] where \(a_L\) and \(a_T\) are the longitudinal and transverse dispersivities.

We assume the dispersivities are independent of \(c_\alpha^\gamma\), and so the condition \({\mathbf{J}_\alpha^\gamma+ \mathbf{J}_\alpha^\delta = 0}\) implies that \(\mathbf{D}_{\alpha}^{\gamma\delta}=\mathbf{D}_{\alpha}^{\delta\gamma}=\mathbf{D}_{\alpha}\). The total dispersion flux becomes: \[\mathbf{J}_{\alpha}^\gamma= -M^\gamma\mathbf{D}_{\alpha} \nabla c_\alpha^\gamma,\quad\text{with}\quad \mathbf{D}_{\alpha} = \begin{bmatrix} a_L & 0 \\ 0 & a_T \\ \end{bmatrix} \frac{u_{x,\alpha}}{\phi}. \label{eq:ficktotal}\tag{11}\]

In the two-phase region, \(\nabla c_\alpha^\gamma=\mathbf{0}\), so the dispersion flux term is only non-zero at the single-phase regions. Furthermore, since we assumed that the injected fluid is in the two-phase region, the dispersion in the gas phase is inconsequential. We can then write the material balance equation (8 ) as: \[\phi \frac{\partial\left(S_g c_g^\gamma+ S_\ell c_\ell^\gamma\right) }{\partial t} + \nabla \cdot \left(c_g^\gamma\mathbf{u}_g + c_\ell^\gamma\mathbf{u}_\ell \right) - \phi \nabla \cdot \left(\mathbf{D}_{\ell} \nabla c_\ell^\gamma\right) = 0, \label{eq:matbal95total2}\tag{12}\]

2.2.3 Darcy’s law↩︎

For multiphase flow in an isotropic porous media, Darcy’s law can be written as:

\[\mathbf{u}_\alpha= \lambda_\alpha\left(\nabla P_\alpha+ \rho_\alpha\mathbf{g} \right), \label{eq:darcy}\tag{13}\] where \(P_\alpha\) is the pressure and \(\rho_\alpha\) the mass density of phase \(\alpha\), \(\mathbf{g}=g\mathbf{e}_x\) is the acceleration vector due to gravity, \(\mathbf{e}_x\) being the unit vector in the vertical direction \(x\), and \(\lambda_\alpha\) the mobility function given by: \[\lambda_\alpha= \frac{k_{r\alpha}}{\mu_\alpha} k, \label{eq:lambda}\tag{14}\] where \(k\) is the absolute permeability, \(k_{r\alpha}\) is the relative permeability function and \(\mu_\alpha\) is the viscosity of phase \(\alpha\).

The total Darcy velocity \(\mathbf{u}_t=\mathbf{u}_g + \mathbf{u}_\ell\) can be written as: \[\mathbf{u}_t = -\lambda_t \nabla P_\ell - \lambda_g \nabla P_c - (\lambda_g \rho_g + \lambda_\ell \rho_\ell) \mathbf{g}, \label{eq:ut}\tag{15}\] where \(\lambda_t=\lambda_g+\lambda_\ell\) and \(P_c\) is the capillary pressure given by: \[\begin{align} P_c = P_g-P_\ell = \sqrt{\frac{\phi}{k}}\;\sigma_{g\ell}\,\cos\theta_{g\ell}\, J_c, \label{eq:Pc} \end{align}\tag{16}\] where \(\sigma_{g\ell}\) is the gas-liquid interfacial tension, \(\theta_{g\ell}\) is the contact angle and \(J_c\) is the Leverett \(J\)-function.

We can write the phase velocities as functions of the total velocity: \[\begin{align} \mathbf{u}_g &= f_g \mathbf{u} - \lambda \nabla P_c + \lambda \Delta\rho \mathbf{g}, \\ \mathbf{u}_\ell &= f_\ell \mathbf{u} + \lambda \nabla P_c - \lambda \Delta\rho \mathbf{g}, \end{align} \label{eq:phasevelo}\tag{17}\] where \(f_\alpha=\lambda_\alpha/\lambda_t\) is the fractional flow of phase \(\alpha\), \(\lambda=\lambda_g\lambda_\ell/\lambda_t\) and \(\Delta\rho=\rho_\ell-\rho_g\).

2.2.4 Mass balance in terms of total velocity↩︎

By substituting Equation 17 in Equation 12 , we arrive at the material balance equation in terms of the total Darcy velocity.

\[\phi \frac{\partial C^\gamma}{\partial t} + \nabla \cdot F^\gamma- \nabla \cdot \left( \mathbf{D}_t^\gamma\nabla C^\gamma\right) = 0,\quad \left\{ \begin{align} &C^\gamma= S_g c_g^\gamma+ S_\ell c_\ell^\gamma, \\ &F^\gamma= \left( f_g c_g^\gamma+ f_\ell c_\ell^\gamma\right) \mathbf{u}_t + \lambda\Delta\rho \Delta c^\gamma\mathbf{g}, \\ &\mathbf{D}_t^\gamma= \lambda\Delta c^\gamma\frac{dP_c}{dC^\gamma}\mathbf{I} + \phi \delta_\ell \mathbf{D}_{\ell}, \\ &\Delta c^\gamma= c_g^\gamma-c_\ell^\gamma, \end{align}\right. \label{eq:matbal95ut}\tag{18}\] where \(C^\gamma\) is the overall concentration, \(F^\gamma\) is the total flow function and \(D_t^\gamma\) the total dissipative tensor for component \(\gamma\). The term \(\delta_\ell\) is the indicator function for the single-phase liquid (or pure-liquid) region defined by: \[\delta_\ell = \left\{ \begin{align} & 0,\quad z^b < x^b_e,\\ & 1,\quad z^b \geq x^b_e. \end{align} \right.\]

In the one-dimensional, incompressible case, the total velocity becomes a known constant, and so the one-dimensional version of Equation 18 becomes quite useful as the base-state equation.

2.2.5 Mass balance in terms of phase pressure↩︎

For the two-dimensional case, \(\mathbf{u}_t\) is not known a priori, and it is preferable to write the mass balance equation in terms of the overall concentration and the pressure. By substituting Equation 13 in Equation 12 , we arrive at the material balance equation in terms of the total liquid phase pressure.

\[\phi \frac{\partial C^\gamma}{\partial t} - \nabla \cdot \left( \bar{\lambda}^\gamma\nabla P_\ell + \bar{\lambda}_\rho^\gamma\mathbf{g} \right) - \nabla \cdot \left(\boldsymbol{\Lambda}_t^\gamma\nabla C^\gamma\right) = 0,\quad \left\{ \begin{align} &\bar{\lambda}^\gamma= \lambda_g c_g^\gamma+ \lambda_\ell c_\ell^\gamma, \\ &\bar{\lambda}_\rho^\gamma= \lambda_g c_g^\gamma\rho_g^\gamma+ \lambda_\ell c_\ell^\gamma\rho_\ell^\gamma, \\ &\boldsymbol{\Lambda}_t^\gamma= \lambda_g c_g^\gamma\frac{dP_c}{dC^\gamma} \mathbf{I} + \phi \delta_\ell \mathbf{D}_\ell. \end{align}\right. \label{eq:matbal95pl}\tag{19}\]

2.2.6 Total velocity↩︎

The assumptions made in Section 2.1 mean that fluid properties are constant inside the two-phase region, and that liquid density is also constant in the pure-liquid region. That leads to the useful result (\(\nabla \cdot \mathbf{u}_t =0\)). In the two-phase region, this can be proven by rewriting Equation 18 with the knowledge that \(c_\alpha^\gamma\) is constant in the two-phase region, to obtain:

\[\phi\frac{\partial S_g}{\partial t} + \nabla \cdot \left[ \mathbf{u}_t \left( f_g + \frac{c_\ell^\gamma}{\Delta c^\gamma} \right) + \lambda \Delta\rho \mathbf{g} + \lambda \nabla P_c\right] = 0. \label{eq:matbal95sg}\tag{20}\]

By subtracting Equation 20 with \(\gamma=b\) from Equation 20 with \(\gamma=a\), we arrive at:

\[\left( \frac{c_\ell^a}{\Delta c^a} - \frac{c_\ell^b}{\Delta c^b} \right) \nabla \cdot \mathbf{u}_t = 0\quad \implies \quad \nabla \cdot \mathbf{u}_t=0.\]

For the pure-liquid region, where \({\mathbf{u}_t=\mathbf{u}_\ell}\), and assuming \({\rho_\ell=constant}\), Equation 12 becomes: \[\phi \frac{\partial\rho_\ell^\gamma}{\partial t} + \nabla \cdot \left( \rho_\ell^\gamma\mathbf{u}_t \right) + \phi \nabla \cdot \mathbf{J}_\ell^\gamma= 0, \label{eq:const95u95proof2}\tag{21}\] and by summing over \(\gamma\), and recognizing that \(\rho_\ell^a + \rho_\ell^b =\rho_\ell\) and \(\mathbf{J}_\ell^a + \mathbf{J}_\ell^b = 0\), we arrive at: \[\rho_\ell \nabla \cdot \mathbf{u}_t = 0.\]

With the result \(\nabla\cdot\mathbf{u}_t=0\) holding for both the two-phase and pure-liquid regions, we can then write: \[-\nabla \cdot \mathbf{u}_t = \nabla \cdot \left( \lambda_t \nabla P_\ell + \lambda_g \nabla P_c + \lambda_\rho \mathbf{g} \right) = 0,\quad \text{with}\quad \lambda_\rho = \lambda_g \rho_g + \lambda_\ell \rho_\ell. \label{eq:gradut}\tag{22}\]

Equations 19 and 22 form the system of equations for the two-dimensional, perturbed flow case, and can be solved for \(C^\gamma\) and \(P_\ell\).

Finally, because the term \(c_\ell^\gamma/\Delta c^\gamma\) is a constant, for the two-phase region, Equation 20 can be rewritten as: \[\phi\frac{\partial S_g}{\partial t} + \nabla \cdot \left( f_g \mathbf{u}_t + \lambda \Delta\rho \mathbf{g} + \lambda \nabla P_c \right) = 0, \label{eq:matbal95sg952}\tag{23}\] which is the well-known balance equation for incompressible, immiscible two-phase flow. This means that, in the case where the initial fluid to be displaced is a saturated liquid, and therefore flow is restricted entirely in the two-phase region, our two-phase, two-component model reduces to the incompressible, immiscible two-phase model.

2.3 Non-dimensionalization↩︎

We can write our equations in dimensionless by the following normalizations: \[\begin{align} \mathbf{x}^*&=\frac{\mathbf{x}}{L}, &\mathbf{u}^*=&\frac{\mathbf{u}}{u_{inj}}, &t^*=&\frac{u_{inj}}{\phi L}t,\\ c_\alpha^{\gamma*}&=\frac{c_\alpha^\gamma}{c_{\ell e}}, &\mu_{\alpha*}=&\frac{\mu_\alpha}{\mu_{g e}}, &\lambda_\alpha^*=&\frac{\mu_{g e}}{k}\lambda_\alpha= \frac{k_{r\alpha}}{\mu_\alpha^*},\\ P_\alpha^*&=\frac{k}{u_{inj}\,\mu_{ge}\, L}P_\alpha, &\rho_\alpha^*=&\frac{kg}{u_{inj}\,\mu_{g e}}\rho_\alpha, &\mathbf{D}_\alpha^*=&\frac{\phi}{u_{inj}\,L} \mathbf{D}_\alpha, \end{align} \label{eq:normalization}\tag{24}\] where \(u_{inj}\) is the injection specific discharge, and \(c_{\ell e}\) and \(\mu_{ge}\) are the molar concentration of the liquid phase and the viscosity of the gas phase, respectively, inside the two-phase region. The length scale \(L\) is chosen as: \[L = \frac{\sqrt{\phi\, k}}{N_{ca}},\quad N_{ca} = \frac{u_{inj}\, \mu_{ge}}{A_{pc}\,\sigma_{g\ell}\,\cos\theta_{g\ell}}, \label{eq:lenghtscale}\tag{25}\] where \(N_{ca}\) is the capillary number and \(A_{pc}\) is a multiplying constant in the \(J_c\) function.

With the normalizations in 24 , and dropping the superscript “*”, Equations 18 , 19 and 15 can be rewritten as: \[\begin{align} &\frac{\partial C^\gamma}{\partial t} + \nabla \cdot F^\gamma- \nabla \cdot \left( \mathbf{D}_t^\gamma\nabla C^\gamma\right) = 0,\quad \left\{ \begin{aligned} &C^\gamma= S_g c_g^\gamma+ S_\ell c_\ell^\gamma, \\ &F^\gamma= \left( f_g c_g^\gamma+ f_\ell c_\ell^\gamma\right) \mathbf{u}_t + \lambda\Delta\rho \Delta c^\gamma\mathbf{e}_g, \\ &\mathbf{D}_t^\gamma= \lambda\Delta c^\gamma\frac{dP_c}{dC^\gamma} \mathbf{I} + \delta_\ell\mathbf{D}_{\ell}, \\ &\Delta c^\gamma= c_g^\gamma-c_\ell^\gamma, \end{aligned}\right. \tag{26} \\ &\frac{\partial C^\gamma}{\partial t} - \nabla \cdot \left( \bar{\lambda}^\gamma\nabla P_\ell + \bar{\lambda}_\rho^\gamma\mathbf{e}_g \right) - \nabla \cdot \left(\boldsymbol{\Lambda}_t^\gamma\nabla C^\gamma\right) = 0,\quad \left\{ \begin{align} &\bar{\lambda}^\gamma= \lambda_g c_g^\gamma+ \lambda_\ell c_\ell^\gamma, \\ &\bar{\lambda}_\rho^\gamma= \lambda_g c_g^\gamma\rho_g^\gamma+ \lambda_\ell c_\ell^\gamma\rho_\ell^\gamma, \\ &\boldsymbol{\Lambda}_t^\gamma= \lambda_g c_g^\gamma\frac{dP_c}{dC^\gamma} + \delta_\ell\mathbf{D}_\ell, \end{align}\right. \tag{27} \\ &\mathbf{u}_t = \lambda_t \nabla P_\ell + \lambda_\rho \mathbf{e}_g + \lambda_c \nabla C^\gamma,\quad \left\{ \begin{align} & \lambda_c = \lambda_g \frac{dP_c}{dC^\gamma}, \\ & \lambda_\rho = \lambda_g \rho_g + \lambda_\ell \rho_\ell, \end{align}\right. \tag{28} \end{align}\] where \(\mathbf{e}_g\) is the unit vector of gravity.

Henceforth, all equations are assumed to be in dimensionless form, unless indicated otherwise.

2.4 Shock in a two-component system↩︎

In Section 2, we briefly introduced the base-state solution as the one-dimensional, travelling-wave solution (to Equation 18 ) that follows the shock between the injected and displaced fluids. This solution, naturally, requires knowledge of the shock front properties, particularly the shock velocity \(v_s\). However, when mass-transfer effects are present, calculating the shock conditions becomes more involved than in the classic immiscible case. [26] has demonstrated how to calculate the flow profile for a two-phase, two-component problem using the method of characteristics. In this section, we provide a brief overview of this procedure. For the two-component case, and inside the two-phase region where the equilibrium concentrations are constants, one can write: \[S_g = \frac{C^\gamma- c_\ell^\gamma}{c_g^\gamma- c_\ell^\gamma}. \label{eq:tieline95sat}\tag{29}\]

Equation 29 shows that phase saturation is a function only of the overall concentration \(C^\gamma\), and therefore we can write that \({F^\gamma=F^\gamma(C^\gamma, S_g(C^\gamma))=F^\gamma(C^\gamma)}\). By ignoring the dissipative terms, we write the purely advective, one-dimensional version of Equation 26 as: \[\frac{\partial C^\gamma}{\partial t} + \frac{d F^\gamma}{d C^\gamma} \frac{dC^\gamma}{d x} = 0.\]

Along a characteristic curve \(\left(x(\eta),t(\eta)\right)\) emanating from a starting position \(x_0\), the characteristic equations are given by: \[\begin{align} \frac{dt}{d\eta} = 1,\quad \frac{dx}{d\eta}=\frac{dF^\gamma}{dC^\gamma},\quad \frac{dC^\gamma}{d\eta}=0. \end{align}\]

So the overall concentration \(C^\gamma\) is constant along the characteristic curves, which in this case are straight lines given by: \[x = \frac{dF^\gamma}{dC^\gamma}\,t + x_{t0},\] where \(x_{t0}\) is the initial position at which \(dF^\gamma/dC^\gamma\) is evaluated.

We assume that the initial system is homogeneous so that, at \(x_{t0}>0\), \(C^\gamma(x)=C^\gamma_i\), \(C^\gamma_i\) being the initial overall concentration of \(\gamma\). At the origin \(x_{t0}=0\), \(C^\gamma=\left[C^\gamma_i, C^\gamma_{inj}\right]\), where \(C^\gamma_{inj}\) is the overall concentration of the injected fluid, and an infinite number of characteristic lines emanate from it. Not all characteristic lines are physically possible, however, and shocks will occur when characteristic lines cross-over each other.

[26] shows that for the two-phase, two-component system, where the initial in-situ fluid is a pure liquid, and the injected fluid is a pure gas, two shocks will appear: a leading shock at the transition from the pure-liquid to the two-phase region, and a trailing shock at the transition from the two-phase to the pure-gas region. However, we are only interested in the leading shock, as the range for \(C^\gamma\) in the solution of the base-state equation (see Section 2.5, Equation 38 ) is \(\left[C^\gamma_{s}, C^\gamma_0\right]\), with \(C^\gamma_s\) the leading shock overall concentration in the two-phase region. To simplify matters, we assume that our injected fluid is a two-phase fluid, and therefore only the leading shock needs to be considered.

The solution for the shock (\(C^\gamma=C^\gamma_s\)) must satisfy three conditions [26]:

1. Jump condition: the Ranking-Hugoniot relation, given by \[v_s = \left. \frac{dx}{dt} \right|_s = \frac{F^\gamma_0 - F^\gamma_1}{C^\gamma_0 - C^\gamma_1}, \label{eq:jump1}\tag{30}\] where the subscripts 0 and 1 indicates the upstream and downstream sides of the shock, respectively;

2. Entropy condition: \[\left. \frac{dx}{dt} \right|_1 \leq \left. \frac{dx}{dt} \right|_s \leq \left. \frac{dx}{dt} \right|_0;\]

3. Velocity constraint: plane wave velocities must decrease monotonically for the solution from the shock composition to the injection composition—i.e. \(v_s\) must be maximal.

Conditions (2) and (3) imply that: \[\left. \frac{dx}{dt} \right|_0 = v_s,\quad \left. \frac{dx}{dt} \right|_1 < v_s,\] and conditions (1), (2) and (3) imply that: \[\left. \frac{dF^\gamma}{dC^\gamma} \right|_s = \left. \frac{dF^\gamma}{dC^\gamma} \right|_0 = \frac{F^\gamma_0 - F^\gamma_1}{C^\gamma_0 - C^\gamma_1}. \label{eq:jump2}\tag{31}\]

We assume that the shock occurs at the two-phase region, where we can write: \[\left. \frac{dF^\gamma}{dC^\gamma} \right|_0 = \left. \frac{dF^\gamma}{dS_g} \frac{dS_g}{dC^\gamma} \right|_0 = u_t \left. \frac{df_g}{dS_g} \right|_0 + \Delta\rho \left. \frac{d\lambda}{dS_g} \right|_0, \label{eq:2p95dFdC}\tag{32}\] and the jump condition becomes: \[u_t \left. \frac{df_g}{dS_g} \right|_0 + \Delta\rho \left. \frac{d\lambda}{dS_g} \right|_0 = \frac{F^\gamma_0 - F^\gamma_1}{C^\gamma_0 - C^\gamma_1}. \label{eq:jump3}\tag{33}\]

We need to evaluate the function \(u_t(C^\gamma)\) at the upstream and downstream compositions to evaluate \(F^\gamma_0\) and \(F^\gamma_1\). As was established earlier, \(u_t\) is constant in both the two-phase and pure-liquid regions. However, if the components undergo a change in volume when they are transferred from one phase to another, the total velocity is not conserved along the shock, meaning that \(u_t(C^\gamma_0)\neq u_t(C^\gamma_1)\). Since we assumed the injected fluid to be in the two-phase region, we simply have \(u_{t0}=u_t(C^\gamma_0)=u_{inj}\). The value for \(u_{t1}=u_t(C^\gamma_1)\) remains undetermined, and therefore an extra equation is needed, and is provided by the fact that Equation 30 is independently true for both \(\gamma=a\) and \(\gamma=b\), which allows us to write: \[\frac{F^a_0 - F^a_1}{C^a_0 - C^a_1} = \frac{F^b_0 - F^b_1}{C^b_0 - C^b_1}. \label{eq:shock952ndcond}\tag{34}\]

Finally, we can solve for (\(C^\gamma_0,u_{t1}\)) using Equations 33 and 34 .

2.5 The base-state equation↩︎

From Reynold’s theorem, the material balance of an extensive quantity \(E\), with density \(e\), in a moving-frame of reference \(\mathbf{\xi} = \mathbf{x} - \mathbf{v} t\), with \(\mathbf{V}^E=\partial\mathbf{x}/\partial t\) as the material element velocity, can be written as: \[\begin{align} \int_{\Omega_{E(t)}} e\, \text{d}\mathbb{V}&= \int_{\Omega_{E(t)}} \frac{\partial e}{\partial t} + \frac{\partial}{\partial x_i} \left(e \frac{\partial x_i}{\partial t}\right) \, \text{d}\mathbb{V}\\ &= \int_{\Omega_{E(t)}} \frac{\partial e}{\partial t} - v_i\frac{\partial e}{\partial \xi_i} + \frac{\partial}{\partial \xi_i} (e V^E_i) \, \text{d}\mathbb{V}, \end{align}\] and the conservation differential equation for an arbitrary elementary volume \(\Omega_{E(t)}\) becomes: \[\frac{\partial e}{\partial t} - \mathbf{v} \cdot \nabla_{\xi}\, e + \nabla_\xi \cdot \left( e \mathbf{V}^E \right) = 0 \label{eq:matbal95emf}\tag{35}\]

From Equation 35 , we observe that the moving-frame equation—in terms of \((\xi,t)\)—is obtained simply by subtracting the term \({\mathbf{v}\cdot\nabla_{\xi}e}\) from the static-frame equation. So for the one-dimensional case of uniform flow, the moving-frame formulation to Equation 26 (where \({e=C^\gamma}\)) that follows the shock front is: \[\frac{\partial C^\gamma}{\partial t} - v_s \frac{\partial C^\gamma}{\partial\xi} + \frac{\partial F^\gamma}{\partial\xi} - \frac{\partial}{\partial\xi} \left( D_{t,xx}^\gamma\frac{\partial z^\gamma}{\partial\xi}\right) = 0. \label{eq:matbal95ut95mf}\tag{36}\]

We could numerically solve Equation 36 to obtain a time-dependent, Buckley-Leverett-type solution for the overall concentration \(C^\gamma\). By assuming that perturbations grow or decay faster than such solution changes with time (the quasi-steady-state assumption), the base state would then be set for a given time \(t\). We opt, however, to use the steady-state version of Equation 36 (where for \({t\to\infty}\), \({\partial C^\gamma/\partial t \to 0}\)) as our base-state equation.

We can rearrange Equation 36 and integrate to obtain: \[\int_{\xi_0}^\xi \frac{d}{d \xi} \left( D_{t,xx}^\gamma\frac{d C^\gamma}{d \xi}\right) d\xi = \int_{\xi_0}^\xi \frac{d}{d \xi} \left(v_s C^\gamma- F^\gamma\right) d\xi, \label{eq:matbal95ut95mf95ss}\tag{37}\] where \(D_{t,xx}^\gamma\) is the diagonal term in \(x\)-direction of tensor \(\mathbf{D}_t^\gamma\).

In the steady-state case, the transition zone around the shock is infinitely large, with \(C^\gamma=C^\gamma_i\) (\(C^\gamma_i\) being the initial overall concentration of \(\gamma\)) and \({\partial_\xi C^\gamma=0}\) when \(\xi_0\to \infty\), and the base-state equation is simplified to: \[D_{t,xx}^\gamma\frac{d C^\gamma}{d \xi} = F_i^\gamma- F^\gamma+ v_s \left( C^\gamma- C_i^\gamma\right), \label{eq:basestate95C}\tag{38}\] where \(F_i^\gamma= F^\gamma(z_i^\gamma)\) and \(C_i^\gamma= C^\gamma(z_i^\gamma)\).

We observe that the solution to Equation 38 varies over the interval \({\left[C^\gamma_{s}, C^\gamma_0\right]}\), and therefore a linear stability analysis around the steady-state solution accounts only for the flow properties across the shock region, ignoring the rest of the Buckley-Leverett solution. This assumption was implicit in many earlier studies [9]–[11]. [12] tested this assumption by perturbing the entire range of saturations in the complete Buckley-Leverett profile and comparing the dispersion relations (growth rate vs wavenumber) for the perturbed saturations. Since any perturbed saturation above or below the shock saturation resulted in smaller growth rates, the authors concluded that instability is indeed governed by the shock region—making the choice of the steady-state solution as the base state a reasonable one. Finally, the expression for the pressure derivative can be derived from Equation 28 : \[\begin{align} \frac{dP_\ell}{d\xi} = -\frac{1}{\lambda_t}\left( u_t - \lambda_\rho \right) - f_g \frac{dP_c}{dC^\gamma} \frac{dC^\gamma}{d\xi}, \label{eq:basestate95pl} \end{align}\tag{39}\] with \(C^\gamma\) the solution for Equation 38 .

2.6 The eigenvalue problem↩︎

Based on Equations 27 and 28 , the two differential equations that define the two-dimensional, moving-frame problem in \((\xi,y)\) is: \[\begin{align} \nabla \cdot \left( \lambda_t \nabla P_\ell + \lambda_\rho \mathbf{e}_g + \lambda_c \nabla C^\gamma\right) &= 0, \tag{40} \\ \frac{\partial C^\gamma}{\partial t} -v_s\frac{\partial C^\gamma}{\partial\xi} - \nabla \cdot \left( \bar{\lambda}^\gamma\nabla P_\ell + \bar{\lambda}_\rho^\gamma\mathbf{e}_g \right) - \nabla \cdot \left(\boldsymbol{\Lambda}_t^\gamma\nabla C^\gamma\right) &= 0 .\tag{41} \end{align}\]

For the linear stability analysis of the system above, we (\(i\)) decompose the primary unknowns (\(C^\gamma,P_\ell\)) in terms of their base-state values and small wave-like disturbances, and (\(ii\)) linearize the differential equations around the same base values. The primary variables (\(C^\gamma,P_\ell\)) can be written as: \[\begin{align} C^\gamma(\xi,y,t) &= \overline{C}(\xi) + \delta C(\xi, y, t) = \overline{C}(\xi) + \hat{c}(\xi)\, e^{iny+\sigma t}, \\ P_\ell(\xi,y,t) &= \overline{P}(\xi) + \delta P(\xi, y, t) = \overline{P}(\xi) + \hat{p}(\xi)\, e^{iny+\sigma t}. \end{align} \label{eq:pertubed95var}\tag{42}\] where \(\overline{C}\) and \(\overline{P}\) are the base-state values, \(\hat{c}\) and \(\hat{p}\) are the disturbed values—with subscript \(\ell\) and superscript \(\gamma\) suppressed for convenience—, \(n\) is the perturbation wave-number and \(\sigma\) is the growth rate of the perturbation. We assume the disturbance eigenfunctions, \(\hat{c}(\xi)\) and \(\hat{p}(\xi)\), decay at infinity:

\[\lim_{\xi\to\pm\infty} \hat{c} = 0\quad \text{and}\quad \lim_{\xi\to\pm\infty} \hat{p}= 0.\]

The linearization of the flow functions around \(\overline{C}\) is done by expanding these functions in a Taylor series and ignoring the non-linear terms: \[f(C^\gamma) \approx f(\overline{C}) + f'(\overline{C}) \delta C,\] where \(f'\) indicates the derivative of \(f\) with respect to \(C^\gamma\).

By substituting Definition 42 in Equations 40 and 41 , linearizing the flow functions and ignoring the quadratic terms in \(\delta C\) and \(\delta P\), and observing that the base-state solutions for Equations 40 and 41 imply \[\begin{align} \frac{d}{d\xi} \left( \lambda_t\frac{d\overline{P}}{d\xi} + \lambda_\rho + \lambda_c \frac{d\overline{C}}{d\xi} \right) &= 0, \\ \frac{d\overline{C}}{dt} - v_s \frac{d\overline{C}}{d\xi} - \frac{d}{d\xi} \left( \overline{\lambda}\frac{d\overline{P}}{d\xi} + \overline{\lambda}_\rho + \Lambda_{t,xx} \frac{d\overline{C}}{d\xi} \right) &= 0, \end{align}\] we arrive at the perturbation equations: \[\begin{align} \frac{d}{d\xi} \left( \lambda_t\frac{d\hat{p}}{d\xi} + \lambda'_t\frac{d\overline{P}}{d\xi}\hat{c} + \lambda'_\rho\hat{c} + \lambda'_c\frac{d\overline{C}}{d\xi}\hat{c} + \lambda_c\frac{d\hat{c}}{d\xi} \right) - n^2\lambda_t\,\hat{p} - n^2\lambda_c\,\hat{c}& = 0, \tag{43} \\ v_s\frac{d\hat{c}}{d\xi} + \frac{d}{d\xi} \left( \overline{\lambda}\frac{d\hat{p}}{d\xi} + \overline{\lambda}'\frac{d\overline{P}}{d\xi}\hat{c} + \overline{\lambda}'_\rho\hat{c} + \Lambda_{t,xx}'\frac{d\overline{C}}{d\xi}\hat{c} + \Lambda_{t,xx}\frac{d\hat{c}}{d\xi} \right)&\nonumber \\ -n^2\overline{\lambda}\,\hat{p} - n^2\Lambda_{t,yy}\,\hat{c}& = \sigma\hat{c}, \tag{44} \end{align}\] where \(\Lambda_{t,xx}\) and \(\Lambda_{t,yy}\) are the diagonal elements of the dispersion tensor \(\boldsymbol{\Lambda}_t\). The derivatives of \(\overline{C}\) and \(\overline{P}\) are computed directly from Equations 38 and 39 . Once again, the superscript \(\gamma\) was dropped for convenience.

Equations 43 and 44 define the eigenvalue problem in terms of eigenfunctions \(\hat{c}\) and \(\hat{p}\), and eigenvalue \(\sigma\). Expanding on the derivatives, we can rewrite the eigenvalue problem as: \[\begin{align} A_2\frac{d^2\hat{p}}{d\xi^2} + A_1\frac{d\hat{p}}{d\xi} + A_0 \hat{p} + B_2\frac{d^2\hat{c}}{d\xi^2} + B_1\frac{d\hat{c}}{d\xi} + B_0 \hat{c} &= 0, \tag{45}\\ E_2\frac{d^2\hat{p}}{d\xi^2} + E_1\frac{d\hat{p}}{d\xi} + E_0 \hat{p} + F_2\frac{d^2\hat{c}}{d\xi^2} + F_1\frac{d\hat{c}}{d\xi} + F_0 \hat{c} &= 0, \tag{46} \end{align}\] where \[\begin{align} A_2 &= \lambda_t,\quad A_1 = \lambda'_t \frac{d\overline{C}}{d\xi},\quad A_0 = -n^2 \lambda_t,\quad B_2 = \lambda_c, \\ B_1 &= \lambda'_t \frac{d\overline{P}}{d\xi} + \lambda'_\rho + 2\lambda'_c\frac{d\overline{C}}{d\xi}, \\ B_0 &= \lambda''_t\frac{d\overline{C}}{d\xi}\frac{d\overline{P}}{d\xi} + \lambda'_t\frac{d^2\overline{P}}{d\xi^2} + \lambda''_\rho\frac{d\overline{C}}{d\xi} + \lambda''_c{\left(\frac{d\overline{C}}{d\xi}\right)}^2 + \lambda'_c\frac{d^2\overline{C}}{d\xi^2} - n^2\lambda_c, \\ E_2 &= \overline{\lambda},\quad E_1 = \overline{\lambda}' \frac{d\overline{C}}{d\xi},\quad E_0 = -n^2 \overline{\lambda},\quad F_2 = \Lambda_{t,xx}, \\ F_1 &= \overline{\lambda}' \frac{d\overline{P}}{d\xi} + \overline{\lambda}'_\rho + 2\Lambda'_{t,xx}\frac{d\overline{C}}{d\xi} + v_s,\\ F_0 &= \overline{\lambda}''_t\frac{d\overline{C}}{d\xi}\frac{d\overline{P}}{d\xi} + \overline{\lambda}'_t\frac{d^2\overline{P}}{d\xi^2} + \overline{\lambda}''_\rho\frac{d\overline{C}}{d\xi} + \Lambda''_{t,xx}{\left(\frac{d\overline{C}}{d\xi}\right)}^2 + \Lambda'_{t,xx}\frac{d^2\overline{C}}{d\xi^2} - n^2\Lambda_{t,yy} - \sigma. \end{align}\]

The linear stability analysis of the moving-frame problem, defined by Equations 40 and 41 , consists in solving numerically the eigenvalue problem (45 46 ) for \(\sigma\), for a given wavenumber \(n\).

2.7 Numerical solution of the eigenvalue problem↩︎

Numerically solving the eigenvalue problem requires careful handling of Equations 45 and 46 around the transition from the two-phase to the pure-liquid region (at \(z^\gamma=x^\gamma_e\)), where the flow functions and base-state derivatives (\(d\overline{C}/d\xi\) and \(d\overline{P}/d\xi\)) are discontinuous.

Solving the base-state equation is straightforward: we set the transition point, where \(C^\gamma_e=C^\gamma(x^\gamma_e)\), at \(\xi=0\), and we integrate Equation 38 in its upstream (\(\xi\leq 0\)) and downstream (\(\xi\geq 0\)) domains separately—where flow is in the two-phase and pure-liquid regions, respectively.

The eigenvalue problem is a more complicated case. The derivatives of \(\hat{c}\) and \(\hat{p}\) are discontinuous at the transition point, and therefore using a conventional method such as the finite-difference method (as employed by [12]) fails to solve our eigenvalue problem. To overcome this difficulty, we adapt the approach of [11], which used the matched initial value problems method. The workflow consists of:

1. Find the far-field solutions of Equations 43 and 44 for their upstream and downstream domains (two-phase and pure-liquid regions, respectively), as functions of arbitrary constants to compute arbitrary initial values, \[\mathbf{m}_{up,iv} = \left\{\hat{p}(\xi_{up}), \hat{p}'(\xi_{up}), \hat{c}(\xi_{up}), \hat{c}'(\xi_{up}) \right\},\quad \xi_{up} < 0,\] and \[\mathbf{m}_{dw,iv} = \left\{\hat{p}(\xi_{dw}), \hat{p}'(\xi_{dw}), \hat{c}(\xi_{dw}), \hat{c}'(\xi_{dw}) \right\},\quad \xi_{dw} > 0,\] where the prime here denotes the derivative with respect to coordinate \(\xi\).

2. Integrate the ordinary differential equations (ODE) for two sets of arbitrary initial values on the upstream domain (\(\xi=[\xi_{up}, 0]\)), and two sets on the downstream domain (\(\xi=[0, \xi_{dw}]\)), to obtain four linearly independent solutions at the transition point \({\xi=0}\): \[\begin{align} \mathbf{m}_{1} &=\left\{\hat{p}_1(0^-), \hat{p}'_1(0^-), \hat{c}_1(0^-), \hat{c}'_1(0^-) \right\},\\ \mathbf{m}_{2} &=\left\{\hat{p}_2(0^-), \hat{p}'_2(0^-), \hat{c}_2(0^-), \hat{c}'_2(0^-) \right\},\\ \mathbf{m}_{3} &=\left\{\hat{p}_3(0^+), \hat{p}'_3(0^+), \hat{c}_3(0^+), \hat{c}'_3(0^+) \right\},\\ \mathbf{m}_{4} &=\left\{\hat{p}_4(0^+), \hat{p}'_4(0^+), \hat{c}_4(0^+), \hat{c}'_4(0^+) \right\}. \end{align}\]

3. The linear combination of each upstream/downstream pair forms the general solution for the upstream/downstream domains, \[\begin{align} \mathbf{m}_{up} &= \alpha_1 \mathbf{m}_1 + \alpha_2 \mathbf{m}_2 = \left\{\hat{p}_{up}^-, \hat{p}_{up}^{\prime-}, \hat{c}_{up}^-, \hat{c}_{up}^{\prime-} \right\},\\ \mathbf{m}_{dw} &= \alpha_3 \mathbf{m}_3 + \alpha_4 \mathbf{m}_4 = \left\{\hat{p}_{dw}^+, \hat{p}_{dw}^{\prime+}, \hat{c}_{dw}^+, \hat{c}_{dw}^{\prime+} \right\}, \end{align}\] where superscripts “-” and “+” indicate that the solution is computed at \(0^-\) or \(0^+\), and \(\alpha_i\) are arbitrary constants.

   Continuity requires that \(\hat{p}_{up}^-=\hat{p}_{dw}^+\) and \(\hat{c}_{up}^-=\hat{c}_{dw}^+\), but such constraint does not apply for the derivatives of \(\hat{p}\) and \(\hat{c}\), since \(\hat{p}_i\) and \(\hat{c}_i\) are discontinuous at \(\xi=0\). To solve the eigenvalue problem we need to derive the jump conditions for the derivatives at \(\xi=0\) (e.g. \(\Delta\left(\hat{p}_i^{\prime}\right)=\lim_{\epsilon\to0}{\left[\hat{p}_i^\prime\right]}_{-\epsilon}^{+\epsilon},\;\epsilon>0\)). More to the point, we seek a pair of expressions that couples the derivatives of \(\hat{p}\) and \(\hat{c}\) to the left and right of \(\xi=0\): \[\hat{p}_i^{\prime+} = f\left( \hat{p}_i^-, \hat{p}_i^{\prime-}, \hat{c}_i^-, \hat{c}_i^{\prime-},\overline{C}(0) \right),\] and \[\hat{c}_i^{\prime+} = g\left( \hat{p}_i^-, \hat{p}_i^{\prime-}, \hat{c}_i^-, \hat{c}_i^{\prime-},\overline{C}(0) \right).\]

   With that, we can match the general solutions at the transition point \(\xi=0\) by defining: \[\mathbf{m}_1^* = \left\{\hat{p}_1^-, \hat{p}_1^{\prime+}, \hat{c}_1^-, \hat{c}_1^{\prime+} \right\},\] and \[\mathbf{m}_2^* = \left\{\hat{p}_2^-, \hat{p}_2^{\prime+}, \hat{c}_2^-, \hat{c}_2^{\prime+} \right\},\] and then write: \[\alpha_1 \mathbf{m}_1^* + \alpha_2 \mathbf{m}_2^* = \alpha_3 \mathbf{m}_3 + \alpha_4 \mathbf{m}_4.\]

   Finally, the non-trivial solution is: \[\begin{vmatrix} \hat{p}_1^- & \hat{p}_2^- & \hat{p}_3^+ & \hat{p}_4^+ \\ \hat{p}_1^{\prime+} & \hat{p}_2^{\prime+} & \hat{p}_3^{\prime+} & \hat{p}_4^{\prime+} \\ \hat{c}_1^- & \hat{c}_2^- & \hat{c}_3^+ & \hat{c}_4^+ \\ \hat{c}_1^{\prime+} & \hat{c}_2^{\prime+} & \hat{c}_3^{\prime+} & \hat{c}_4^{\prime+} \end{vmatrix} = 0. \label{eq:nontrivial}\tag{47}\]

2.7.1 Far-field solutions↩︎

We are interested in the far-field asymptotic solutions of the perturbation equations (43 and 44 ) for the far-upstream (\(\xi<\xi_{up}\)) and far-downstream (\(\xi>\xi_{dw}\)) domains, where \(d\overline{C}/d\xi\) approaches zero, and the base-state variable \(\overline{C}\) and the flow functions have approximately constant values. The perturbation equations can then be expressed for the far-field domains as:

\[\begin{align} \lambda_t \left(\frac{d^2\hat{p}}{d\xi^2} - n^2 \hat{p} \right) + \lambda_c\frac{d^2\hat{c}}{d\xi^2} + \left( \lambda'_t\frac{d\overline{P}}{d\xi} + \lambda'_\rho + \lambda'_c\frac{d\overline{C}}{d\xi}\right) \frac{d\hat{c}}{d\xi} - n^2 \lambda_c \hat{c}& = 0, \tag{48}\\ \overline{\lambda} \left(\frac{d^2\hat{p}}{d\xi^2} - n^2 \hat{p} \right) + \Lambda_{t,xx}\frac{d^2\hat{c}}{d\xi^2} + \left( v_s + \overline{\lambda}'\frac{d\overline{P}}{d\xi} + \overline{\lambda}'_\rho + \Lambda'_{t,xx}\frac{d\overline{C}}{d\xi} \right) \frac{d\hat{c}}{d\xi}& \nonumber\\ -\left(\sigma+ n^2 \Lambda_{t,yy} \right) \hat{c}& = 0. \tag{49} \end{align}\]

By substituting Equation 48 in 49 , we obtain the homogeneous ODE for overall concentration \(\hat{c}\):

\[\quad\theta_2\frac{d^2\hat{c}}{d\xi^2} + \theta_1\frac{d\hat{c}}{d\xi} + \theta_0\hat{c} = 0,\] where \[\def\m{\frac{\overline{\lambda}}{\lambda_t}} \begin{align} &\theta_2 = \Lambda_{t,xx} - \m \lambda_c, \\ &\theta_1 = v_s + \overline{\lambda}'\frac{d\overline{P}}{d\xi} + \overline{\lambda}'_\rho + \Lambda'_{t,xx}\frac{d\overline{C}}{d\xi} - \m \left( \lambda'_t\frac{d\overline{P}}{d\xi} + \lambda'_\rho + \lambda'_c\frac{d\overline{C}}{d\xi}\right), \\ &\theta_0 = -\sigma+ n^2 \left( \m \lambda_c - \Lambda_{t,yy} \right). \end{align}\]

The boundary condition for the far-field domains are: \[\lim_{\xi \to -\infty} \hat{c}_{up}\left(\xi\right) = 0,\quad \text{and}\quad \lim_{\xi \to +\infty} \hat{c}_{dw}\left(\xi\right) = 0.\]

The general solutions—up to an arbitrary constant—for \(\hat{c}\) in the far-upstream and far-downstream domains are: \[\begin{align} \hat{c}_{up} &= A_1 e^{r_{pos} \xi}, \\ \hat{c}_{dw} &= A_2 e^{r_{neg} \xi}, \end{align}\] where \(r_{pos}\) and \(r_{neg}\) are the positive and negative roots of the characteristic equation, \[\theta_2 r^2 + \theta_1 r + \theta_0 = 0.\]

By substituting the exponential solution for \(\hat{c}\) in Equation 48 , we obtain the non-homogeneous ODE for the pressure \(\hat{p}\): \[\frac{d^2\hat{p}}{d\xi^2} - n^2\frac{d\hat{p}}{d\xi} = -\frac{A}{\lambda_t} \left[ \lambda_c r^2 + \left( \lambda'_t\frac{d\overline{P}}{d\xi} + \lambda'_\rho + \lambda'_c\frac{d\overline{C}}{d\xi}\right) r - n^2 \lambda_c \right] e^{r\xi}.\]

The boundary conditions for the far-field domains are: \[\lim_{\xi \to -\infty} \hat{p}_{up}\left(\xi\right) = 0,\quad \text{and}\quad \lim_{\xi \to +\infty} \hat{p}_{dw}\left(\xi\right) = 0,\] and the general solutions for \(\hat{p}\) in the far-upstream and far-downstream domains are: \[\begin{align} \hat{p}_{up} &= B_1 e^{n \xi} + \frac{A_1}{\lambda_t}\, \frac{\lambda_c r^2_{pos} + \left( \lambda'_t\frac{d\overline{P}}{d\xi} + \lambda'_\rho + \lambda'_c\frac{d\overline{C}}{d\xi}\right) r_{pos} - n^2 \lambda_c }{r^2_{pos} - n^2}\, e^{r_{pos}\xi}, \\ \hat{p}_{dw} &= B_2 e^{-n \xi} + \frac{A_2}{\lambda_t}\, \frac{\lambda_c r^2_{neg} + \left( \lambda'_t\frac{d\overline{P}}{d\xi} + \lambda'_\rho + \lambda'_c\frac{d\overline{C}}{d\xi}\right) r_{neg} - n^2 \lambda_c }{r^2_{neg} - n^2}\, e^{r_{neg}\xi}. \end{align}\]

2.7.2 Coupling the left and right-side derivatives↩︎

The coupling expressions for the left and right-side derivatives (jump conditions) are obtained by integrating the perturbation equations across the transition point (\(\xi=0\)): \[\begin{align} {\left[ \lambda_t\frac{d\hat{p}}{d\xi} + \lambda'_t\frac{d\overline{P}}{d\xi}\hat{c} + \lambda'_\rho\hat{c} + \lambda'_c\frac{d\overline{C}}{d\xi}\hat{c} + \lambda_c\frac{d\hat{c}}{d\xi} \right]}_{-\Delta\xi}^{+\Delta\xi}\nonumber\\ -\int_{-\Delta\xi}^{+\Delta\xi} \left( n^2\lambda_t\,\hat{p} + n^2\lambda_c\,\hat{c} \right) d\xi &= 0, \label{eq:coupling1a} \end{align}\tag{50}\] \[\begin{align} {\left[ v_s\hat{c} + \overline{\lambda}\frac{d\hat{p}}{d\xi} + \overline{\lambda}'\frac{d\overline{P}}{d\xi}\hat{c} + \overline{\lambda}'_\rho\hat{c} + \Lambda_{t,xx}'\frac{d\overline{C}}{d\xi}\hat{c} + \Lambda_{t,xx}\frac{d\hat{c}}{d\xi} \right]}_{-\Delta\xi}^{+\Delta\xi}\nonumber\\ -\int_{-\Delta\xi}^{+\Delta\xi} \left[ n^2\overline{\lambda}\,\hat{p} + \left(\sigma+ n^2\Lambda_{t,yy}\right)\hat{c} \right] d\xi &= 0. \label{eq:coupling1b} \end{align}\tag{51}\]

Assuming that the flow functions and the base-state derivatives (\(d\overline{C}/d\xi\) and \(d\overline{P}/d\xi\)) can be discontinuous but are bounded in the neighbourhood of \(\xi=0\), and that \(\hat{c}\) is continuous, by taking the limit \(\Delta\xi\to0\), we obtain: \[\begin{align} {\left[ \lambda_t\frac{d\hat{p}}{d\xi} + \lambda'_t\frac{d\overline{P}}{d\xi}\hat{c} + \lambda'_\rho\hat{c} + \lambda'_c\frac{d\overline{C}}{d\xi}\hat{c} + \lambda_c\frac{d\hat{c}}{d\xi} \right]}_{0^-}^{0^+} &= 0, \tag{52}\\ {\left[ \overline{\lambda}\frac{d\hat{p}}{d\xi} + \overline{\lambda}'\frac{d\overline{P}}{d\xi}\hat{c} + \overline{\lambda}'_\rho\hat{c} + \Lambda_{t,xx}'\frac{d\overline{C}}{d\xi}\hat{c} + \Lambda_{t,xx}\frac{d\hat{c}}{d\xi} \right]}_{0^-}^{0^+} &= 0. \tag{53} \end{align}\]

Using the fact that flow functions \(\lambda_t\) and \(\overline{\lambda}\) are continuous at \(\xi=0\), we can write: \[\begin{align} \lambda_{t0} \left(\left.\frac{d\hat{p}}{d\xi}\right|_+ - \left.\frac{d\hat{p}}{d\xi}\right|_-\right) + \hat{c}_0 \beta_1 + \lambda_{c+}\left.\frac{d\hat{c}}{d\xi}\right|_+ - \lambda_{c-}\left.\frac{d\hat{c}}{d\xi}\right|_- &= 0, \tag{54} \\ \overline{\lambda}_0 \left(\left.\frac{d\hat{p}}{d\xi}\right|_+ - \left.\frac{d\hat{p}}{d\xi}\right|_-\right) + \hat{c}_0 \beta_2 + \Lambda_{t,xx+}\left.\frac{d\hat{c}}{d\xi}\right|_+ - \Lambda_{t,xx+}\left.\frac{d\hat{c}}{d\xi}\right|_- &= 0, \tag{55} \end{align}\] where \[\begin{align} \beta_1 &= \lambda'_{t+} \left.\frac{d\overline{P}}{d\xi}\right|_+ - \lambda'_{t-}\left.\frac{d\overline{P}}{d\xi}\right|_- + \lambda'_{\rho+} - \lambda'_{\rho-} + \lambda'_{c+}\left.\frac{d\overline{C}}{d\xi}\right|_+ - \lambda'_{c-}\left.\frac{d\overline{C}}{d\xi}\right|_-, \\ \beta_2 &= \overline{\lambda}'_+ \left.\frac{d\overline{P}}{d\xi}\right|_+ - \overline{\lambda}'_-\left.\frac{d\overline{P}}{d\xi}\right|_- + \overline{\lambda}'_{\rho+} - \overline{\lambda}'_{\rho-} + \Lambda_{t,xx+}'\left.\frac{d\overline{C}}{d\xi}\right|_+ - \Lambda_{t,xx-}'\left.\frac{d\overline{C}}{d\xi}\right|_-, \end{align}\] where the minus and plus sign subscripts indicate that functions are evaluated at \(\xi=0^-\) and \(\xi=0^+\), respectively. For the flow functions and base state derivatives, that translates to being evaluated at \(\overline{C}(z^\gamma=x^{\gamma}_{e-})\) and \(\overline{C}(z^\gamma=x^{\gamma}_{e+})\), respectively.

By substituting Equation 54 into 55 , we can write the relationship between the left-side and right-side derivatives of \(\hat{c}\): \[\left.\frac{d\hat{c}}{d\xi}\right|_+ = \frac{\Lambda_{t,xx-} - \alpha\lambda_{c-}}{\Lambda_{t,xx+} - \alpha\lambda_{c+}} \left.\frac{d\hat{c}}{d\xi}\right|_- - \frac{\hat{c}_0}{\Lambda_{t,xx+} - \alpha\lambda_{c+}} \left( \beta_2 - \alpha\beta_1 \right), \label{eq:coupling4a}\tag{56}\] where \(\alpha= \lambda_{t0}/\overline{\lambda}_0\).

The relationship between the left-side and right-side derivatives of \(\hat{p}\) are simply: \[\left.\frac{d\hat{p}}{d\xi}\right|_+ = \left.\frac{d\hat{p}}{d\xi}\right|_- - \frac{1}{\lambda_{t0}} \left( \hat{c}_0 \beta_1 + \lambda_{c+}\left.\frac{d\hat{c}}{d\xi}\right|_+ - \lambda_{c-}\left.\frac{d\hat{c}}{d\xi}\right|_- \right). \label{eq:coupling4b}\tag{57}\]

Equations 56 and 57 allows us to solve Equation 47 for eigenvalue \(\sigma\).

2.7.3 Capillary function choice↩︎

Since the mechanical dispersion on the two-phase region is zero, we can write the dissipative function on the upstream domain as: \[D_{t,xx}^{up} = \lambda \Delta c^\gamma\frac{dP_c}{dS_g}\frac{dS_g}{dC^\gamma}=\frac{\lambda_g \lambda_\ell}{\lambda_t}\frac{dP_c}{dS_g},\] where we used the result that \(dC^\gamma/dS_g = \Delta c^\gamma\) on the two-phase region.

Now, for a given \(P'_c\) function, it is possible that \[\lim_{C^\gamma\to C^\gamma_e} D_{t,xx} = \lim_{S_g\to 0} \frac{\lambda_g \lambda_\ell}{\lambda_t} \frac{dP_c}{dS_g} = 0,\] and from Equation 38 , \[\lim_{D_{t,xx}\to 0} \frac{d\overline{C}}{d\xi},\] becomes unbounded.

To avoid this result, we take inspiration in previous work—e.g.[11] and [12]—and choose the Leverett J-function in Equation 16 to have its derivative in the form: \[\frac{dJ_c}{dS_g} = \frac{A_{pc}}{k_{rg}\, k_{r\ell}},\] where \(A_{pc}\) is a multiplying constant.

The dimensionless capillary force derivative can then be written as: \[\frac{dP_c}{dS_g} = \frac{1}{k_{rg}\, k_{r\ell}} = \frac{1}{M} \frac{1}{\lambda_g\, \lambda_\ell}, \label{eq:dPc}\tag{58}\] where is the viscosity ratio \(M=\mu_{\ell e}/\mu_{ge}\) on the two-phase region.

This choice for \(P'_c\) means that the dissipative function \(d\overline{C}/d\xi\) is bounded in the neighbourhood of \(C^\gamma_e\). As pointed out by [12], it also has a physically meaningful relationship with the fluid saturations in the two-phase region, and has the added benefit of allowing comparison with past investigations on linear stability of immiscible, two-phase systems.

3 Results↩︎

We refer to the two components in the binary fluid as component \(a\) (the lighter component) and component \(b\) (the heavier component). Phase equilibrium properties were calculated using the Peng-Robinson (PR78) equation of state, and viscosities were computed using the corresponding states (CS) model described in [27].

Two simplifying assumptions were made: (\(i\)) the liquid phase density remains constant in the pure-liquid region (as discussed in Section 2.1), and (\(ii\)) that the liquid viscosity in the pure-liquid phase follows and exponential function given by: \[\mu_\ell = \mu^b e^{A_\mu \left(z^b - 1\right)},\quad \text{for}\quad z^b > x^b_e, \label{eq:mu95exp}\tag{59}\] with \[A_\mu = \frac{\ln\mu^b - \ln\mu_{\ell e}}{1 - x^b_e},\] where \(\mu^b\) and \(\mu_{\ell e}\) are the viscosities of the pure component \(b\) and the liquid phase at the two-phase region, respectively, calculated using the CS model.

When both gravity and mass transfer are present, the instability of our two-phase, two-component model is controlled by the dimensionless functions \(k_{rg}\), \(k_{r\ell}\) and \(\mu_\ell\), and the dimensionless quantities \(M=\mu_{\ell e}/\mu_{ge}\), \(\Delta_\rho^*\), \(D_{\ell,xx}^*\), \(R_a=a_L/a_T\) and \(z^\gamma_i\). The high number of dimensionless variables obscure the physical interpretation of the system. Therefore, to facilitate our analysis, we define a dimensional base case, and work with dimensional and non-dimensional variables as appropriate for each case. To avoid confusion, dimensionless variables are explicitly denoted with a superscript “*” in the following sections.

We note that, unlike the immiscible case studied by earlier authors [9]–[12], for partially miscible flow the maximum growth rate (\(\sigma_{max}\)) and cutoff wavenumber (\(n_{cut}\)) values do not scale linearly with the capillary number (\(N_{ca}\)). This is readily apparent when we observe—from Equations 24 and 25 —that the \(N_{ca}\) appears explicitly in the expression for the dimensionless mechanical dispersion tensor \(\mathbf{D}_\ell^*\). In fact, as we show in Section 3.3, \(\mathbf{D}_\ell^*\) has a very non-linear impact on \(\sigma_{max}\) and \(n_{cut}\).

3.1 Base case↩︎

For our base case fluid we chose a CO\(_2\)-decane mixture. Table 1 lists the EOS parameters for the components. Flow conditions and the porous medium properties for the base case are listed in Table 2. For the purpose of evaluating the fluid properties, pressure and temperature were assumed to be constant at \({50\,bar}\) and \({50^\circ C}\), respectively (see Section 2.1).

Figure 3 shows the gas and liquid properties for varying overall compositions (in terms of \(z^b=1-z^a\)), and a comparison between the viscosities, molar density and mass density with and without the simplifying assumptions above. It can be seen that, for our binary mixture, the assumption of constant liquid density in the pure-liquid region is reasonable, as the liquid density varies only slightly with composition, and that the exponential function for the liquid viscosity is a good approximation of the CS model. The viscosity ratio in the two-phase region is \({M_e=\mu_{\ell e}/\mu_{ge}\approx23}\), while the viscosity ratio between the initial liquid and injected gas is \({M_i=\mu_{\ell i}/\mu_{gi}\approx40}\).

The relative permeability (RP) parameters for the Corey-type base case RP are listed in Table 3. The capillary pressure is a function of the RPs, as described in Section 2.7.3, and the multiplying constant \(A_{pc}\) is also listed in Table 3. Figure 4 shows the RPs and the corresponding \(P_c\) function.

In the following sections, most analyses consist of examining the sensitivity of stability to a set of chosen parameters. Unless otherwise stated, parameters not explicitly varied are held at their base case values.

3.2 Initial composition↩︎

The initial fluid composition (\(z^b_i\)) determines the amount of gas that can be dissolved into the liquid phase. If \(z^b_i \leq x^b_e\), then the initial fluid is a saturated liquid and the displacement process occurs solely in the two-phase region—which, as shown in Section 2.2.6, means that the flow equations reduce to those of a two-phase, immiscible case. If \(z^b_i > x^b_e\), then limited mass transfer is present in the displacement process, with flow being characterized as partially miscible, and transition from the two-phase to the pure-liquid region occurs at the shock. The lower the value of \(x^b_i\), the greater the amount of mass transfer.

Figures 5 and 6 show the base-state profiles for two initial compositions: (\(i\)) \(z^b_i=0.582\), just slightly below \(x^b_e\), with \(S_{gi}=0.001\), and (\(ii\)) \(z^b_i=1\). For the former case, where flow is immiscible, the solution of the base-state equation (38 ) is smooth, while for the latter, where flow is partially miscible (PM), the base-state solution is not smooth and the base-state derivatives are discontinuous.

The results of the linear stability analysis for the two initial compositions are shown in Figures 7 and 8, respectively, which show the plot of the growth rate vs wavenumber curves (dispersion relation) and the eigenfunctions associated with the maximum growth rate (\(\sigma_{max}\)). We point out that while the eigenfunctions for the saturated case are continuous and smooth, those for the PM case are not (see Section 2.7).

An interesting observation is that the maximum growth rate (\(\sigma_{max}\)) for the PM case is lower than \(\sigma_{max}\) for the saturated case, while the opposite is true for the cutoff wavenumber (\(n_{cut}\)) and most unstable mode (\(n_{max}\)).

Many factors influence stability during the transition from saturated to undersaturated conditions (achieved by increasing \(z^b_i\)), making it difficult to isolate their individual effects. The most obvious impact is the increase in initial fluid viscosity and, consequently, in \(M_i\). Another important effect is that increasing \(z^b_i\) enhances mass transfer between phases, which in turn alters the shock configuration.

To understand the impact of \(z^b_i\) on the shock configuration, we first turn to the simplest case where the volume occupied by component \(a\) remains unchanged as it dissolves in the liquid phase (i.e.\(a\) has the same density as \(b\)), and the total velocities upstream (\(u_{t0}\)) and downstream (\(u_{t1}\)) of the shock are equal. From Figure 5a, we see that an increase in the initial concentration of component \(b\) (\(C^b_i\)) results in a decrease of \(dF^b/dC^b\) at the shock, which in turn leads to a decrease of the overall concentration of \(b\) at the shock (\(C^b_s\)). Because \(C^b\) and the gas saturation are inversely related, an increase in \(C^b_i\) means a higher saturation at the shock (\(S_{gs}\)), and finally, with the assumption that \(u_t\) is constant, mass conservation requires that a higher \(S_{gs}\) results in a lower shock front velocity \(v_s\). In the more realistic case where component \(a\) decreases in volume when transferring to the liquid phase, the total velocity downstream (\(u_{t1}\)) of the shock is lower than \(u_{t0}\) (see Figure 6a), and \(S_{gs}\) is increased (and \(v_s\) decreased) even further. These trends can be clearly observed in Figure 9.

The impact of these changes on stability can be inferred from the first-order approximation of the growth rate given by [9] for immiscible displacements: \[\frac{\sigma}{n} = v_s \frac{\lambda_{t0} - \lambda_{t1}}{\lambda_{t0} + \lambda_{t1}} \label{eq:1o95stab}\tag{60}\] where \(\lambda_{t0}\) and \(\lambda_{t1}\) are the total mobilities upstream and downstream of the shock, respectively. From Equation 60 , we can expect a decrease in \(v_s\) to reduce instability. Conversely, increasing \(S_g\) at the shock should increase the mobility contrast \({\lambda_{t0}-\lambda_{t1}}\), which should increase instability. Additionally, the dispersion tensor \(\mathbf{D}_\ell\) is linearly dependent on \(u_{t1}\) and therefore \(D_{\ell x}\) and \(D_{\ell y}\) decrease with \(z^b_i\). So we have competing effects, and the variation of \(\mathbf{D}_\ell\) further complicates the picture.

3.3 Dispersivities↩︎

Figure 10 shows the impact on stability of increasing mass transfer (by increasing \(z^b_i\)) for different values of for the longitudinal dispersivity \(a_L\), with the ratio \({a_L/a_T=5}\). The dispersion relations are, in general, approximately symmetrical, and \(n_{max}\) is approximately half of \(n_{cut}\) for most cases. We see that \(\sigma_{{max}}\) and \(n_{cut}\) exhibit a strongly non-linear behaviour. For the case where \(a_L=0.001\,m\), both \(\sigma_{max}\) and \(n_{cut}\) initially increase, and then decrease, with \(z^b_i\). For the base case (\(a_L=0.01\,m\)), \(\sigma_{max}\) and \(n_{cut}\) follow different trends, while the case \(a_L=0.1\,m\) has the largest values for \(\sigma_{max}\) and \(n_{cut}\) of all three cases, with \(\sigma_{max}\) varying only slightly with \(z^b_i\).

Intuitively, we would expect that higher values of \(a_L\) would lead to a more stable system with lower values of \(\sigma_{max}\) and \(n_{cut}\), but Figure 10 shows the opposite for the cases selected. To understand the physics behind these results, we first need to understand the effect of transverse and longitudinal dispersivities on stability. To that end, we run a simple sensitivity on \(a_L\) and \(a_T\). First, transverse dissipative forces—be it from mechanical dispersion or capillary pressure—stabilize the displacement process by laterally dispersing the perturbations, and have a straightforward impact on stability. This is shown in Figure 11, where \(a_T\) is varied while \(a_L\) is held constant at \(0.01\,m\). Three curves are shown for different values for the capillary pressure multipliers \(A_{pc}\). As expected, increasing \(a_T\) decreases both \(\sigma_{max}\) and \(n_{cut}\), although the effect is quite limited for the cases where \(A_{pc}=0.01\;\text{and}\;0.001\), when mechanical dispersion is small compared to the capillary dispersion.

Longitudinal dissipative forces influence stability by smoothing the shock profile in the flow direction, reducing the overall concentration gradient of the base-state. Because it is the variation of the fluid properties along the displacement front that induces instability, one would expect that smoothing the \(\overline{C}\)-profile would reduce instability. However, Figure 12, where the sensitivity for \(a_L\) is plotted, shows a more complex picture. Two main observations emerge: (\(i\)) there appears to exist an intermediate value of \(a_L\) for which instability is maximized, and (\(ii\)) \(\sigma_{max}\) and \(n_{cut}\) exhibit asymptotic behaviour for low values of \(a_L\) and large values of \(P_c'\).

Observations (\(i\)) and (\(ii\)) are more easily seen in Figure 13, where we plot the same simulations from Figure 12 in dimensionless form, with \(\sigma_{max}^*\) and \(n_{cut}^*\) as functions of the dimensionless dispersion coefficient \(D^*_{\ell,xx}\), recast in Equation 61 for convenience. Note that the expression for \(D^*_{\ell,xx}\) incorporates the ratio between mechanical and capillary dispersion.

\[D^*_{\ell,xx} = \frac{\mu_g}{\sqrt{\phi k}} \frac{a_L\, u_{t1}}{A_{pc}\, \sigma_{g\ell}\,\cos\theta_{g\ell}}. \label{eq:Dstar}\tag{61}\]

Figure 13 shows that there exists an intermediate value of \(D^*_{\ell,xx}\) for which instability is maximized (\(D^{*\,max}_{\ell,xx}\)). For lower values of \(D^*_{\ell,xx}\), where mechanical dispersion is small compared to capillary dispersion, \(\sigma_{max}^*\) and \(n_{cut}^*\) behave asymptotically. For higher values of \(D^*_{\ell,xx}\), where mechanical dispersion dominates, \(\sigma_{max}^*\) and \(n_{cut}^*\) fall steadily towards zero. This explains the behaviour of \(\sigma_{max}\) and \(n_{cut}\) seen in Figure 10.

First, by calculating the values of \(D^*_{\ell,xx}\) for the three Figure 10 cases (\(a_L=0.001\), \(0.01\) and \(0.1\,m\), with \(D^*_{\ell,xx}\approx 10 a_T\) for the chosen parameters), we see from Figure 13 that, in agreement with Figure 10, the case \(a_L=0.1\,m\) (\(D^*_{\ell,xx}=1\)) is indeed the most unstable case, and case \(a_L=0.001\,m\) (\(D^*_{\ell,xx}=0.01\)) the most stable of the three.

Second, as discussed in Section 3.2, increasing \(z^b_i\) has three effects on displacement: the shock velocity decreases, the mobility contrast increases (because \(S_{gs}\) increases) and \(D^*_{\ell,xx}\) decreases (because \(u_{t1}\) decreases). As seen in Figure 13, the impact on \(\sigma_{max}^*\) and \(n_{cut}^*\) of decreasing \(D^*_{\ell,xx}\) is strongly non-linear: if \(D^*_{\ell,xx}<D^{*\,max}_{\ell,xx}\), then decreasing \(D^*_{\ell,xx}\) should reduce instability, while if \(D^*_{\ell,xx}>D^{*\,max}_{\ell,xx}\), decreasing \(D^*_{\ell,xx}\) increases instability. This explains why, for the \(a_L=0.001\,m\) (\(D^*_{\ell,xx}=0.01\)) case, \(\sigma_{max}^*\) and \(n_{cut}^*\) mostly decrease with \(z^b_i\), while the opposite happens for \(a_L=0.1\,m\) (\(D^*_{\ell,xx}=1\)). Further discussion on the impact of \(D^*_{\ell,xx}\) on stability is given in Section 4.

3.4 Miscibility↩︎

Another variable that influences the extent of mass transfer between the two phases is the level of gas/liquid miscibility. A useful device to alter miscibility without changing the pressure, temperature or the components of the system is to alter the binary interaction coefficients (BIC) of the EOS model. Higher BIC values correspond to a lower content of component \(a\) in the liquid phase at equilibrium, \(x^a_e\)—which we use as a proxy for system miscibility in this section.

Figure 14 shows the impact of increasing miscibility (or increasing \(x^a_e\)) on the shock configuration and the viscosity ratio of the two-phase region, \(M_e\). It shows a similar picture to that of Figure 9, where we saw the changes in \(v_s\), \(u_{t1}\) and \(S_{gs}\) due to increasing mass transfer. We note that the decrease in \(M_e\) with \(x^a_e\) further induce higher values of \(S_{gs}\).

As mentioned earlier, changes in the shock configuration have competing effects on stability, making it difficult to evaluate their impact. On the other hand, the decrease in \(M_e\)—which reduces the mobility contrast along the displacing front—has a clear stabilizing effect. Figure 15 shows that, in all cases, the effect of increasing miscibility is to reduce \(\sigma_{max}\).

The behaviour of \(n_{cut}\) is more complex. For the cases with \(a_L=0.1\), 0.3, and 0.003, \(n_{cut}\) also decreases with increasing miscibility. However, for the base case (\(a_L=0.01\)), \(n_{cut}\) follows a very different trend at high values of \(x^a_e\), increasing at an accelerating rate with \(x^a_e\). We also observe that the approximate relation \(n_{cut}\approx2n_{max}\) breaks down at high \(x^a_e\), reflecting increased asymmetry in the dispersion relation curves. This divergent behaviour suggests a complex interplay between the various mechanisms governing system stability.

3.5 Viscosity ratio↩︎

To analyse the impact of the viscosity ratio on stability, we run a series of cases where the viscosity of the liquid phase (in both the two-phase and pure-liquid regions) is altered by a multiplying constant, while keeping all other fluid properties of the base case unchanged. We compare two scenarios, one where immiscible displacement (IM) occurs, and one where flow is partially miscible (PM). The metric compared is the viscosity ratio between the initial and invading fluids, \(M_i\).

Figure 16 compares the shock properties of the IM and PM cases, for various values of \(M_i\), and in Figure 17 we compare the instability of the IM and PM, for different values of \(a_L\). The main observation in Figure 16 is that, for the IM scenario, the shock velocity \(v_s\) grows considerably with \(M\), while in the PM scenario, \(v_s\) fall slightly with \(M\). This behaviour of \(v_s\) explains why in the PM cases the slope of \(\sigma_{max}\) decreases with \(M\), while for the IM case \(\sigma_{max}\) increases approximately linearly at high values of \(M_i\) (as observed by [12]). This dampening effect is stronger for higher values of \(a_L\), and we see that although the case with high \(a_L\) (\(0.001\,m\)) is the most unstable at low \(M_i\), at high values of \(M_i\) it is the most stable. Finally, this dampening effect is less pronounced for the \(n_{cut}\), and for the low \(a_L\) case (\(0.001\,m\)) \(n_{cut}\) increases with \(M_i\) in a similar rate for both the IM and PM cases.

3.6 Gravity↩︎

Before we proceed with the analysis of gravity effects on stability, we look at the impact of the density variation alone (absent gravity, horizontal flow). The key here is that the density ratio between the gas and liquid phases in equilibrium (\(\rho_{ge}/\rho_{\ell e}\)) dictates the amount of change in volume that component \(a\) undergoes when transferring from the gas to the liquid phase. The higher the density (mass or molar) of component \(a\), the higher the density of the gas phase. When \(c^a_e=c^b_e\), the phase densities match and component \(a\) experiences no volume change upon transfer between the phases.

As discussed earlier in Section 3.2, the lower the density of component \(a\), the higher the shock saturation \(S_{gs}\), and the lower the shock and downstream velocities, \(v_s\) and \(u_{t1}\). In Figure 18 we see these trends in all scenarios of \(M_i\) (10, 100 and 1000), and can observe that \(u_{t1}=1\) when \(\rho_{ge}=\rho_{\ell e}\).

In Figure 19, we see that the impact of volume change is higher for larger values of \(M_i\), and that for cases \(M_i=100\) and 1000, the effect of increasing \(\rho_{ge}/\rho_{\ell e}\) is increasing \(\sigma_{max}\) and \(n_{cut}\). For the case \(M_i=10\), we find that the case with no volume change (\(\rho_{ge}/\rho_{\ell e}=1\)) is more stable than the base case (\(\rho_{ge}/\rho_{\ell e}=0.14\)). As discussed earlier, increasing \(v_s\) increases instability, and decreasing \(S_{gs}\) reduces the mobility contrast, decreasing instability. For the \(M_i=10\) case, the latter appears to have a larger impact on overall instability.

Finally, we test the impact of gravity forces on stability. The two flow configurations studied are the vertically upward, where \(\Delta\rho^*>0\), and the vertically downward (\(\Delta\rho^*<0\)) cases. In the case where \(\Delta\rho^*=0\), both the upward and downward flow cases match the horizontal flow scenario. The variable \(\Delta\rho^*\) was introduced in Equation 17 and is recast in its dimensionless form in Equation 62 for clarity: \[\Delta\rho^* = \frac{k\, g}{u_{inj}\, \mu_{ge}} \left(\rho_{\ell e} - \rho_{ge}\right). \label{eq:Delta95rho}\tag{62}\]

We analyse the impact of \(\Delta\rho^*\) on stability for three different values of \(M_i\) (40, 200 and 500), and for two types of displacement, IM and PM. The results are shown in Figure 20.

For the IM scenarios, viscous and gravity forces interact in a more additively fashion compared to the PM scenarios. For the former, we see that the curves \(\sigma_{max}\) vs \(\Delta\rho^*\) for different viscosity ratios have similar monotonic shapes and are approximately parallel, whereas the PM curves are strongly non-linear functions of \(\Delta\rho^*\). Interestingly, for the PM cases with high \(M\) (100 and 1000), where viscous forces become more relevant than gravity, higher values of \(|\Delta\rho^*|\) resulted in lower values of \(\sigma_{max}\) even for upward scenario (\(\Delta\rho^*>0\)). As discussed earlier, such behaviour is explained by mass exchange effects, with decreasing gas density decreasing instability (see Figure 19). For the upward PM case with \(M_i=40\), gravity does initially contribute to instability, with \(\sigma_{max}\) and \(n_{cut}\) reaching maximum values at around \(\Delta\rho^*\approx20\) and \(\Delta\rho^*\approx25\), respectively. For the downward scenario, increasing \({|\Delta\rho^*|}\) always has a stabilizing influence, which is larger for the PM case (compared to the IM case), since \(\Delta\rho^*<0\) not only introduces stabilizing gravity forces, but also mass-transfer effects.

4 Discussion↩︎

Comparison between immiscible and partially miscible displacements shows that mass transfer has a stabilizing effect on instability. In most cases presented in Section 3, adding or enhancing mass transfer reduces the maximum growth rate (\(\sigma_{max}\)), particularly for displacements with large viscosity ratios (see Figure 17). Mass transfer also dampens gravity-induced instability and, in the analysed cases, exerts a larger stabilizing effect than gravity at large density and viscosity contrasts (see Figure 20).

The impact on the cutoff wavenumber—\(n_{cut}\), at which the displacement reaches marginal stability—is neither as pronounced nor as straightforward as that on \(\sigma_{max}\). For example, Figure 17 shows that for the base case (\(a_L=0.01\)), \(n_{cut}\) values for the IM and PM cases are quite similar. The most apparent effect of mass transfer is the alteration of the shock configuration—reducing the front velocity and increasing the invading fluid saturation—along with a reduction in the equilibrium viscosity ratio \(M_e\). While the latter unequivocally promotes stability, the former introduces competing effects: reducing the shock velocity decreases the growth rate, but increasing the gas saturation at the shock increases mobility contrast, thereby promoting instability. One hypothesis is that while \(v_s\) has less impact on \(n_{cut}\) than on \(\sigma_{max}\), \(S_{gs}\) exhibits the opposite behaviour, and if the decrease in \(M_e\) caused by mass transfer is offset by the increase in \(S_{gs}\), then both IM and PM scenarios would exhibit similar values for \(n_{cut}\).

An interesting comparison can be made with results from [25]. Their study—which focused on marginal stability (i.e., \(n_{cut}\)) for a case with highly adverse displacement (\(M_i\approx280\))—concluded that mass transfer promotes stability by reducing total flow velocity and viscosity contrast at the displacing front. While these mechanisms are consistent with those identified in the present study, we found that, although they significantly impact \(\sigma_{max}\), their effects on \(n_{cut}\) are much less pronounced. In some cases, the PM scenario even exhibits larger \(n_{cut}\) values than the IM scenario. A key difference in their model may explain this discrepancy: [25] assumed an idealized step-saturation displacement profile, thereby neglecting the effects of mass transfer on the shock profile.

Arguably, the most interesting result of our linear stability analysis for two-component compositional displacement is shown in Figure 13, which illustrates the impact of \(D^*_{\ell,xx}\) on instability. The behaviour at extreme values of \(D^*_{\ell,xx}\) is physically plausible and intuitive: at very low values, capillary dispersion dominates and further reductions in \(D^*_{\ell,xx}\) produce negligible changes in the displacement process, resulting in asymptotically small variations in \(\sigma_{max}\) and \(n_{cut}\). At large values of \(D^*_{\ell,xx}\), increasing dispersion stabilizes the displacement, as expected. However, the existence of a most dangerous value \(D^{*,max}_{\ell,xx}\)—where the interaction between capillary and mechanical dispersion conspires to increase instability—is difficult to explain. To our knowledge, this complex behaviour has not been reported in the literature on viscous fingering and linear stability analysis.

The closest observation appears in [19], who compared models with constant and variable diffusion coefficients for nanoparticles. In the variable-diffusion case, the diffusion coefficient decreases rapidly with increasing concentration of other solutes, resulting in a sharper base-state concentration profile compared to the constant-diffusion case (Figure 4 in [19]). Interestingly, the dispersion relation showed that the variable-coefficient model—where concentration gradients are steeper—predicted lower values for both \(\sigma_{max}\) and \(n_{cut}\). The authors attributed this finding to reduced maximum viscosity gradients in the transition zone under variable dispersion. Critically, the variable-coefficient model in their study affected only the concentration profile of stabilizing nanoparticles added to increase the invading fluid viscosity. In contrast, in our study the mobility contrast is either unaffected by \(D^*_{\ell,xx}\) in the two-phase region or decreases monotonically with \(D^*_{\ell,xx}\) in the single-phase region.

Four main caveats apply to the analysis presented in this paper. First, one main assumption in this analysis—and in several others in the literature—is that the relative permeability curves are identical for all values of \(M\). This assumption is unlikely to hold in practice, as demonstrated by [5]: although the viscous fingering experiments presented in their study showed a trend of decreasing shock saturations with increasing \(M\), the authors had to select different sets of relative permeability curves for each experiment to capture these changes. Additionally, the linear scaling of \(\sigma_{max}\) with \(M\) for large \(M\), shown by [12] (and in Figure 17), is actually valid only for a particular choice of relative permeability curves.

Second, our analysis in Section 3.6 of gravitational effects on stability accounted for how the dimensionless density contrast—\({\Delta\rho^*}\), given by Equation 62 —impacts the relative strength of viscous versus gravity forces and mass-transfer effects, which counteract gravity in upward displacement (see Figure 20). However, a key assumption was that miscibility between the phases (or, equivalently, \(x^b_e\)) remains independent of density contrast, whereas in many practical cases these properties are inversely correlated—higher density contrasts typically correspond to lower miscibility. This assumption does not invalidate our analysis of how gravity and mass transfer affect instability, but it means Figure 20 does not fully represent real systems where mass transfer extent varies with density contrast.

Third, the assumption of constant liquid density in the pure-liquid region implies that gravity effects on instability arose from stabilizing/destabilizing forces in the two-phase region only. This simplifying assumption was reasonable for our base case given that \(\rho_\ell\) varied little with \(z^a_i\). However, this is not universally true, and it would be valuable to evaluate the impact of variable density in the liquid region. In fact, such a formulation may not add significant complexity to the two-component model. Although the continuity equation \({\nabla\cdot\mathbf{u}_t=0}\) would no longer be valid in the pure-liquid domain, the total mass balance equation—obtained by summing Equation 19 over \(\gamma\)—could be used instead, since it is independent of Equation 19 and valid for non-constant densities.

Finally, mechanical dispersion is strongly scale-dependent, with laboratory experiments showing values of \(a_L\) around \(0.1\) to \(10\,\text{mm}\), whereas field studies using numerical simulations show values up to \(100\,\text{m}\) [29]. Indeed, mechanical dispersion is fundamentally an upscaling technique [24], and larger scales generally correspond to larger dispersivities. Since the dispersion relations may span vastly different scales depending on flow parameters, the appropriate dispersivity values may vary considerably. At larger scales and in weakly wet porous media, it may even be reasonable for mechanical dispersion to dominate over capillary effects.

5 Conclusions↩︎

A mathematical model for two-phase, two-component displacement in porous media was successfully developed and implemented for linear stability analysis. The model accounts for partial miscibility, gravity effects, capillary forces, and mechanical dispersion in a system where a gaseous fluid displaces a heavier, more viscous liquid. By deriving the jump conditions for the eigenfunctions’ derivatives—which are discontinuous at the transition from two-phase to pure-liquid flow at \(\xi=0\)—we were able to solve the eigenvalue problem numerically using the matched initial value problem (MIVP) method.

Results show that mass transfer in the partially miscible (PM) scenario reduces the equilibrium viscosity ratio (\(M_e\)), the shock front velocity (\(v_s\)), and the downstream total velocity (\(u_{t1}\)), while increasing the invading fluid saturation at the front (\(S_{gs}\)). The overall effect is stabilization of the displacement front, and increasing miscibility reduces both the maximum growth rate (\(\sigma_{max}\)) and the cutoff wavenumber (\(n_{cut}\)) in most cases.

Comparison between immiscible (IM) and partially miscible (PM) displacement scenarios shows that mass transfer in PM flow dampens the effect of increasing viscosity ratio (\(M\)) on \(\sigma_{max}\), but has a much less pronounced impact on \(n_{cut}\), particularly when dispersivities (\(a_L\) and \(a_T\)) are small.

For upward displacement, the added instability from gravity forces is mitigated by mass-transfer effects, or even suppressed—the key insight being that the density contrast driving gravity-induced instability also enhances the extent of phase mass transfer. For downward displacement, both gravity and mass transfer act as stabilizing mechanisms.

An intriguing and unexpected finding is the existence of a most dangerous value for the dimensionless longitudinal dispersion coefficient, \(D^{*,max}_{\ell,xx}\), where both \(\sigma_{max}\) and \(n_{cut}\) are maximized. This suggests a complex interaction between capillary forces, mechanical dispersion, and instability, and should be further investigated.

No supplementary information or material.

The authors thank Arne Skauge for his helpful remarks. The authors acknowledge the PhD scholarship for Paulo L. K. Caetano Chang provided by Petrobras. Kundan Kumar acknowledges funding from the Centre of Sustainable Subsurface Resources (CSSR), grant nr., supported by the Research Council of Norway, research partners NORCE Norwegian Research Centre and the University of Bergen, and user partners Equinor ASA, Harbour Energy Norge AS, Sumitomo Corporation, Earth Science Analytics, GCE Ocean Technology, and SLB Scandinavia.

Declarations↩︎

Paulo Lee Kung Caetano Chang received funding from Petrobras. Kundan Kumar received funding from Centre of Sustainable Subsurface Resources, Grant ID 331841.

The authors have no competing interests to declare that are relevant to the content of this article.

Not applicable.

Not applicable.

Not applicable.

Not applicable.

The python code developed to generate and solve the eigenvalue problem will be made available on request.

P. L. K. Caetano Chang: conceptualization; methodology; software; writing—original draft; writing—review & editing. K. Kumar: conceptualization; methodology; writing—review & editing.

References↩︎