Linear Stability Analysis of Buoyancy-Driven Convection of Two-Phase Fluid in a Porous Medium

Abstract

Understanding the onset and stability conditions of buoyancy-driven convection in two-phase immiscible flows through the porous medium is paramount for a wide range of subsurface engineering applications, such as density-driven \(\textrm{CO}_2\) convective dissolution in saline aquifers and water flooding for enhanced oil recovery. Motivated by these imperative processes, we present a rigorous linear stability analysis of buoyancy-driven multiphase immiscible thermal convection in a porous medium. Starting from the extended Darcy’s law, the local thermal equilibrium assumption and the Oberbeck–Boussinesq approximation, and incorporating appropriate relative permeability and capillary pressure models, we derive the mathematical equations governing two-phase Rayleigh–Darcy convection. Steady-state profiles for saturation and temperature are established, and small perturbations are introduced to formulate the linearized perturbation equations, which are solved as a generalized eigenvalue problem. A dedicated numerical solver is developed to map the neutral stability curves, critical wavenumbers, and critical Rayleigh–Darcy numbers \(Ra_c\) governing convective onset. We systematically evaluate the influence of key dimensionless parameters, including the gravity number G, the non-wetting to wetting phase ratios of volumetric thermal expansion coefficient \(B_n\), viscosity \(M_n\), heat conductivity \(\Lambda _n\), and volumetric heat capacity \(C_n\), the porous medium ratios (\(\Lambda _s\), \(C_s\)), as well as capillary pressure parameters and boundary conditions. The results demonstrate that G exerts a modest stabilizing effect, elevating the threshold for convective onset. Furthermore, the stability threshold increases with decreasing \(B_n\) and \(C_n\), and with increasing \(M_n\), \(\Lambda _n\), and \(\Lambda _s\), whereas \(C_s\) has negligible impact. Additionally, stronger capillary pressure further stabilizes the system, thereby raising the threshold for convective initiation. Finally, by combining the individual power-law dependencies obtained from the single-parameter sensitivity analyses, we establish an empirical scaling relation between the critical Rayleigh–Darcy number and the dominant governing parameters, \(Ra_c \propto B_n^{-1} M_n \Lambda _n C_n^{-1} \Lambda _s\), providing a practical predictive compass for the onset of buoyancy-driven multiphase convection in subsurface storage and geo-energy systems.