Global Buckley–Leverett Framework

• The Global Buckley–Leverett Framework is an advanced model that extends classical multiphase flow theory by incorporating hyperbolic conservation laws, kinetic methods, and data-driven surrogates.
• It employs nonlocal regularizations, phase-field techniques, and multiscale upscaling to ensure rigorous solvability and accurate simulation in heterogeneous and fractured media.
• The framework has practical applications in petroleum engineering, carbon sequestration, and hydrology by merging theoretical insights with cutting-edge numerical and stochastic approaches.

The Global Buckley–Leverett Framework encompasses the evolution, extension, and coupling of the classical Buckley–Leverett (BL) theory for multiphase flow in porous media, supporting modern mathematical, computational, and physical approaches. It integrates hyperbolic conservation-law theory, kinetic/entropy methods, compressible-to-incompressible limits, phase-field regularizations, stochastic representations, upscaling in heterogeneous domains, deep-learning surrogates, and conservative multicomponent balances with rigorous solvability and well-posedness results. These developments directly impact fields including petroleum engineering, carbon sequestration, contaminant hydrology, and geomechanics.

1. Classical Buckley–Leverett Formulation and Analytical Structure

The BL model describes immiscible two-phase displacement (e.g., water/oil, gas/oil) by combining Darcy flow and mass conservation. For scalar saturation $u(t, x) \in [0,1]$ in domain $\Omega \subset \mathbb{R}^d$, the canonical equation is: $$\phi \frac{\partial u}{\partial t} + \nabla \cdot (q f(u)) = 0$$ where $q$ is total flux, $\phi$ is porosity, and $f(u)$ is the fractional-flow function derived from phase mobilities: $$f(u) = \frac{\lambda_w(u)}{\lambda_w(u) + \lambda_o(u)},\quad \lambda_i(u) = k_{r,i}(u)/\mu_i$$ Shock and rarefaction development are predicted via Rankine–Hugoniot conditions: $$V_s = \frac{f(S_R) - f(S_L)}{S_R - S_L}, \qquad \text{with entropy admissibility}\ \lambda(S_L) > V_s > \lambda(S_R)$$ where $\lambda(u) = df/du$ serves as the characteristic velocity. Classical solutions proceed via the method of characteristics and Welge construction, capturing the nonlinear wave structure of saturation fronts.

2. Generalizations: Nonlocal Regularizations, Solvability, and Pseudo-Parabolic Couplings

Global solvability and regularity of BL-type systems require regularization. The generalized Stokes–Buckley–Leverett (S–BL) system (Chemetov et al., 2010) augments classical BL with Brinkman (diffusive) and Cattaneo–Forchheimer (time-delay) terms: $$\begin{cases} \partial_t u + \nabla_x \cdot (v g(u)) = 0 \ T\partial_t v - \nu \Delta_x v + h(u) v = - \nabla_x p \ \nabla_x \cdot v = 0 \end{cases}$$ Here, $g(u), h(u)$ derive from phase-relative permeabilities, $T > 0$ is the time-delay, and $\nu > 0$ the viscosity parameter. The global existence of weak entropy solutions is established for arbitrary time horizon $T > 0$ with kinetic formulation and compactness arguments. The kinetic method involves introducing $$f^\varepsilon(t, x, v) = \chi_{u^\varepsilon(t, x) > v}$$ and controlling entropy-dissipation via nonnegative Radon measures.

Nonlocal and fractional regularizers (diffusive and conservative) further generalize BL (Burczak et al., 2015): $$u_t + \partial_x f(u) + \nu\,\Lambda^\alpha u + \mu\,\Lambda^\beta u = 0$$ where $\Lambda^s = (-\Delta)^{s/2}$, and global existence is proved for regularizing order $\max\{\alpha, \beta\} > 1$, or under explicit smallness conditions if $\max\{\alpha, \beta\} = 1$. Blow-up for insufficient regularization is supported by numerical evidence. Entropy (Lyapunov) functionals underpin $L^\infty$–in–time bounds.

3. Upscaling and Physical Heterogeneity: Macroscale Formulations

Heterogeneities and vertical layering in porous media are incorporated via upscaling of Darcy flow (Benham et al., 2020). The layer-averaged mobilities $\lambda_w^{\mathrm{eff}}(S)$ and $\lambda_n^{\mathrm{eff}}(S)$ interpolate between capillary and viscous limits, parameterized by $N_c$, the capillary number:

• Viscous limit: $\lambda_{w}^{\mathrm{visc}}(S) = k_{rw}(S)$; $\lambda_{n}^{\mathrm{visc}}(S) = k_{rn}(S)/M$.
• Capillary limit: $\lambda_{w}^{\mathrm{cap}}(S)$, $\lambda_{n}^{\mathrm{cap}}(S)$ are spatial integrals over permeability-weighted relative permeabilities.

Effective fractional-flow $F(S)$ is constructed as: $$F(S) = \frac{\lambda_w^{\mathrm{eff}}(S)}{\lambda_w^{\mathrm{eff}}(S) + \lambda_n^{\mathrm{eff}}(S)}$$ The global BL equation is preserved: $$\frac{\partial S}{\partial t} + \frac{\partial}{\partial x} F(S) = 0$$ but wave structure (shock speed, rarefaction) and displacement dynamics are significantly altered. In 1D and axisymmetric domains, front speeds shift by 13%–44%, with heterogeneous predictions matching field CO₂ breakthrough for sequestration.

4. Extensions: Compressibility, Free Boundaries, Multicomponent and Fractured Media

Compressible approximations with stiff pressure $p_\gamma(u_\gamma) = u_\gamma^\gamma$ yield rigorous incompressible Hele–Shaw-type free-boundary limits for BL systems (Gomes et al., 15 Apr 2024). In the $\gamma \to \infty$ regime: $$\partial_t u - \nabla \cdot (\Theta(u) \nabla p) = 0,\qquad \operatorname{div}\left( \frac{1}{h(u)} \nabla p \right) = 0$$ subject to $p(u - 1) = 0$ almost everywhere, resulting in a moving front that obeys Rankine–Hugoniot across the interface and admits global (BV) weak solutions.

Multicomponent BL formulations (Tantardini et al., 9 Nov 2025) introduce Maxwell–Stefan diffusion and dynamic capillarity, yielding pseudo-parabolic transport: $$-\nabla \cdot [M_\alpha(S)\, \tau_\alpha(S)\, \nabla(\partial_t S_\alpha)]$$ A single global-pressure equation coordinates total flux, while strict component-wise mass conservation is maintained using equation-of-state flashes and conservative balances. Stress-sensitive permeability, non-Darcy fracture flow (Forchheimer), buoyancy, and capillary drifts are rigorously incorporated, providing a backbone for storage and transport simulations in fractured media.

5. Stochastic and Data-Driven Approaches

Physical uncertainty in $k(\mathbf x)$ is represented as a log-Gaussian random field, expanded via Karhunen–Loève and projected onto generalized polynomial chaos (gPC) bases (Pettersson et al., 2016). The stochastic Galerkin approach expands both $S(\mathbf x, t, \boldsymbol{\xi})$ and $\mathbf u(\mathbf x, \boldsymbol{\xi})$, yielding a hyperbolic system for mode coefficients. Locally reduced basis methods (thresholding) allow adaptive complexity without modifying system matrices, resulting in robust, accurate front mean/variance capture versus Monte Carlo, with speedups $6\times$–$10\times$.

Deep learning implementations employ physics-informed neural networks (PINNs) and physics-informed GANs (PI-GANs) (Fraces et al., 2020). Neural network surrogate $S(x, t; \theta)$ employs automatic differentiation for PDE-residual regularization: $$r(x, t; \theta) = \partial_t \hat{S} + \partial_x f(\hat{S})$$ Loss functions combine data misfit and residuals at collocation points. GANs introduce latent variables for ensemble statistics and improve shock capture. The frameworks accommodate both forward (surrogate) and inverse (parameter-learning) BL problems, requiring minimal data for extrapolation and uncertainty quantification. Hyperparameters and network depth/width must be carefully tuned for strong discontinuities.

6. Numerical and Algorithmic Advances

Phase-field (Cahn–Hilliard) regularizations recast BL equations into energy-gradient flows (Cogswell et al., 2016), supporting large-timestep stable solvers: $$\phi\,\frac{\partial S_w}{\partial t} + \nabla\cdot (\mathcal{F}_w \mathbf{v}_t) = \nabla \cdot (M(S_w) \nabla \mu_w) + q_w$$ with convex-splitting and homotopy $\lambda$ scaling ($\lambda \to 0$ recovers BL limit), solved by semi-implicit multigrid methods. Timesteps $\sim10^8$ are possible, with near-linear computational speedup, enabling efficient simulation of highly heterogeneous reservoirs.

7. Physical Regimes, Applicability, and Practical Guidelines

The incompressibility assumption for BL models holds for water/oil with density changes $\Delta \rho/\rho < 10\%$ under typical pressure ranges (Zhuravljov et al., 2017). For gases (CO₂, CH₄, N₂), this assumption fails when $\Delta \rho/\rho$ exceeds $20\%$ (low initial pressure or high injection pressure), requiring fully compressible simulators. In higher dimensions, mismatch between incompressible and compressible solutions is magnified—errors can double in 2D. A simple criterion: $$\frac{\rho_g(P_{\text{ini}} + \Delta p) - \rho_g(P_{\text{ini}})}{\rho_g(P_{\text{ini}})} < 0.10 \Longrightarrow \text{BL applies}$$ otherwise, the compressible extension must be used.

The Global Buckley–Leverett Framework thus encompasses rigorous mathematical generalizations, effective upscaling in heterogeneous and fractured environments, advanced numerical and data-driven methods, and detailed regimes of physical applicability. It serves both as a theoretical foundation and a practical toolbox for multiphase and multicomponent flow modeling in complex geologic media.