Enthalpy-Porous Media Phase Change Model

• The model is a framework that represents phase change in porous media by integrating latent heat into the enthalpy variable for both numerical and physical closures.
• It couples thermal transport, phase change, and momentum balance across scales, employing methods such as fixed-grid, REV-scale, and phase-field formulations.
• Numerical implementations range from finite-volume and lattice Boltzmann to implicit Newton–KKT schemes, enabling robust simulation of complex multiphase flows.

Searching arXiv for the cited papers to ground the article in recent literature. The enthalpy–porous media phase-change model denotes a class of formulations in which phase change in a porous medium is represented through an enthalpy variable, together with porous-medium closures for momentum, transport, or upscaled effective behavior. In the cited literature, this class includes the fixed-grid enthalpy porosity method for melting with natural convection, REV-scale solid–liquid phase-change models based on generalized non-Darcy flow, double-distribution-function lattice Boltzmann implementations, a persistent-variable isenthalpic flash formulation for thermal compositional multiphase flow with phase separation, a quasistatic thermomechanical freezing–melting system with hysteresis, and a phase-field model for evaporation with Darcy-scale homogenization (Schüller et al., 2018, Liu et al., 2015, Liu et al., 2016, Gavioli et al., 2021, Ghosh et al., 2021, Lipovac et al., 3 Dec 2025, He et al., 2018).

1. Governing balance laws in porous media

A central feature of these models is the simultaneous treatment of thermal transport, phase change, and porous-medium flow. The specific set of primary unknowns depends on the formulation. In the fixed-grid enthalpy porosity model for water–ice melting, the fields are a single velocity field $u$, pressure $p$, temperature $T$, and a scalar phase field $\phi(x,t)\in[0,1]$, identified with the local liquid volume fraction $f_L$ (Schüller et al., 2018). In the thermal compositional multiphase formulation of Lipovac et al., the primary flow–transport unknowns are

$X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$

where $p$ is pore pressure, $H$ the specific enthalpy of the mixture, and $z_\xi$ the overall mass fraction of component $\xi$ (Lipovac et al., 3 Dec 2025). In the quasistatic thermomechanical model of Gavioli and Krejčí, the unknowns are capillary pressure $p$0, solid displacement $p$1, absolute temperature $p$2, and liquid phase fraction $p$3 (Gavioli et al., 2021).

The mass, momentum, and energy balances likewise vary by scale and physics. Under the Boussinesq approximation, the enthalpy porosity model uses incompressible continuity,

$p$4

Navier–Stokes momentum with a Darcy-type drag in the mushy zone,

$p$5

and an energy equation in enthalpy form,

$p$6

This is a mixture-theory formulation on a fixed mesh (Schüller et al., 2018).

At the REV scale, Liu and He formulate incompressible flow in a rigid, isotropic porous medium with generalized non-Darcy momentum,

$p$7

together with an enthalpy-form energy equation

$p$8

Here the porous-medium closure enters through permeability $p$9, porosity $T$0, and the Forchheimer coefficient $T$1 (Liu et al., 2015).

The thermal compositional model extends these balances to multicomponent, multiphase flow with phase separation. Summing phase-wise component balances yields

$T$2

while the energy balance is written as

$T$3

The term $T$4 arises from summing phase-wise enthalpy balances and converting $T$5 (Lipovac et al., 3 Dec 2025).

The most general model in the set is the quasistatic thermomechanical system of PDEs for a deformable porous solid. There, mass, momentum, and energy are coupled to a relaxed Stefan-type phase dynamics inclusion and a Preisach pressure-saturation hysteresis relation. This places phase change within a Biot-type coupling with visco-elasto-plastic deformation rather than within a rigid-medium framework (Gavioli et al., 2021).

2. Enthalpy as thermal state variable

Across these formulations, enthalpy is the device by which sensible and latent heat are combined. In the enthalpy porosity model,

$T$6

with $T$7 the latent heat and $T$8 the liquid fraction. Rewriting the energy balance in terms of $T$9 gives a latent-heat source term

$\phi(x,t)\in[0,1]$0

showing explicitly how the enthalpy formulation accounts for melting and solidification (Schüller et al., 2018).

In Liu and He’s REV-scale solid–liquid formulation,

$\phi(x,t)\in[0,1]$1

and

$\phi(x,t)\in[0,1]$2

is added as a time-dependent source term in the thermal lattice Boltzmann equation. The liquid fraction $\phi(x,t)\in[0,1]$3 is therefore the latent-heat carrier in the REV-scale enthalpy balance (Liu et al., 2015).

The three-dimensional DDF-MRT model writes an effective volumetric enthalpy as

$\phi(x,t)\in[0,1]$4

leading to the REV-scale energy equation

$\phi(x,t)\in[0,1]$5

Here porosity multiplies the thermal storage and latent term directly (Liu et al., 2016).

In the deformable porous-medium model, the total specific enthalpy is introduced as

$\phi(x,t)\in[0,1]$6

or in differential form,

$\phi(x,t)\in[0,1]$7

This enthalpy variable is embedded in a setting with phase relaxation, dissipation, and hysteresis (Gavioli et al., 2021).

The compositional formulation of Lipovac et al. uses mixture enthalpy $\phi(x,t)\in[0,1]$8 as both a PDE unknown and a thermodynamic equilibrium constraint. A specific point made there is that latent heat is built in via the jump in $\phi(x,t)\in[0,1]$9 across phase boundaries, so no explicit extra term is needed once $f_L$0 is modelled correctly (Lipovac et al., 3 Dec 2025). This does not contradict the explicit-source formulations above; it indicates that the placement of the latent-heat contribution depends on how enthalpy is defined and closed.

The evaporation model extends the same idea to diffuse-interface multiphase flow. It defines mixture enthalpy as

$f_L$1

and also notes the split

$f_L$2

At Darcy scale, periodic homogenization yields

$f_L$3

with $f_L$4 the saturation and $f_L$5 the effective heat capacity (Ghosh et al., 2021).

3. Phase representation, mushy zones, and thermodynamic closure

One major distinction among enthalpy–porous media phase-change models is how the phase state is represented. In fixed-grid enthalpy porosity methods, the liquid fraction is a scalar field tied directly to temperature through a piecewise-linear mushy-zone law: $f_L$6 To suppress motion in partially or fully solidified regions, the momentum equation includes the classical Kozeny–Carman drag

$f_L$7

so that $f_L$8 in the liquid and $f_L$9 in the solid (Schüller et al., 2018).

REV-scale porous-medium phase-change models use essentially the same mushy-zone closure for $X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$0. In Liu and He’s two-dimensional model,

$X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$1

and the interface is traced through the liquid fraction determined by the enthalpy method (Liu et al., 2015). The review by Liu et al. presents the same closure in total-enthalpy form and emphasizes that the solid–liquid interface is handled without explicit front-tracking; $X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$2 varies smoothly across a few grid cells, and no geometric front-tracking or re-meshing is needed (He et al., 2018).

A different closure appears in thermal compositional flow with phase separation. For given $X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$3, the local fluid must split into up to $X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$4 phases by solving an isenthalpic flash: $X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$5 subject to overall-mass constraints, unity constraints, and nonnegativity of the phase-mass fractions $X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$6. The local variables are

$X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$7

A Lagrangian/KKT treatment produces fugacity equalities, an enthalpy match $X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$8, and complementarity conditions $X=[p,\;H,\;z_2,\ldots,z_{n_C}]^\top,$9. The persistent-variable idea fixes $p$0 upfront, retains all unknowns whether that phase is present or not, and enforces $p$1 with complementarity instead of performing extra stability tests (Lipovac et al., 3 Dec 2025).

The quasistatic thermomechanical model adopts yet another phase representation. Its phase-fraction variable $p$2 obeys a relaxed Stefan-type inclusion,

$p$3

with the driving terms coupled to pressure, hysteresis, and deformation (Gavioli et al., 2021).

The phase-field evaporation model replaces a sharp interface by a diffuse layer of width $p$4 and evolves the phase field $p$5 with an Allen–Cahn equation. In the limit $p$6, matched asymptotics recover the classical Stefan problem with jump conditions and Laplace pressure (Ghosh et al., 2021). This suggests that “enthalpy–porous media phase-change model” is not a single closure law but a framework in which the phase variable may be a liquid fraction, a persistent set of phase fractions and compositions, a hysteretic phase fraction, or a diffuse phase field.

4. Numerical formulations and solver architectures

The numerical realization of enthalpy-based phase change in porous media spans finite-volume, Newton–KKT, and lattice Boltzmann strategies. In the enthalpy porosity model for melting with natural convection, the PDEs are solved by finite-volume discretization on a fixed mesh, with implicit backward-Euler treatment of all convective and diffusive terms and PIMPLE pressure–velocity coupling as implemented in OpenFOAM’s buoyantBoussinesqPimpleFoam. The nonlinear latent-heat source is treated by an outer fixed-point, or Gauss–Seidel, iteration over $p$7: given $p$8, solve the linearized energy equation for $p$9, compute the consistent temperature, and update $H$0 with

$H$1

using a relaxation $H$2. The outer $H$3-update loop is iterated to a tight tolerance $H$4 (Schüller et al., 2018).

Lipovac et al. formulate a fully implicit thermal compositional problem in which the discrete residual is

$H$5

with $H$6 collecting pressure, component, and energy residuals, and $H$7 collecting flash residuals. They use cell-centered finite volumes for $H$8, multi-point flux approximation for $H$9 and $z_\xi$0, upwind for convective $z_\xi$1-flux, backward Euler in time, and a semi-smooth Newton step on the block system for $z_\xi$2. Because $z_\xi$3 is block-diagonal, a Schur-complement reduction yields a reduced system for $z_\xi$4, while the local flash subproblem remains cellwise and parallelizable (Lipovac et al., 3 Dec 2025).

The lattice Boltzmann literature realizes the same enthalpy closures with double-distribution schemes. Liu and He use D2Q9-MRT for the flow field and D2Q5-MRT for the temperature field with latent-heat source; the source $z_\xi$5 enters the thermal equation, and an inner iteration ensures consistent $z_\xi$6 (Liu et al., 2015). Their three-dimensional extension uses a D3Q19 or D3Q15 lattice for flow and D3Q7 for temperature/enthalpy, with forcing terms for porous drag and buoyancy in the flow equation and the latent-heat source $z_\xi$7 in the thermal equation (Liu et al., 2016). The review by Liu et al. presents an enthalpy-based BGK or MRT “recipe” in which one LB subsolver advances the flow field and another advances the enthalpy field; after streaming and collision, $z_\xi$8 is used to update $z_\xi$9 and $\xi$0 pointwise, and no outer iteration is required because the latent heat is already absorbed into $\xi$1 and $\xi$2 (He et al., 2018).

The phase-field evaporation model is numerically different again: it solves mass, momentum, vapor transport, phase-field evolution, and energy at pore scale, then derives Darcy-scale effective parameters by periodic homogenization through cell problems for $\xi$3, $\xi$4, $\xi$5, and $\xi$6 (Ghosh et al., 2021).

5. Validation, benchmark behavior, and solver performance

Validation studies show that enthalpy-based porous-media phase-change models can reproduce interface motion, convection-driven asymmetry, and multiscale transport, but the observed accuracy depends strongly on the physical closure. In the water–ice validation study, melting occurred against an isothermal vertical wall in a rectangular plexiglas cavity of inner size $\xi$7. A binary Mumford–Shah segmentation of gray-scale images was used to extract the water–ice interface, and nine probe points equally spaced in $\xi$8 were used to compare simulated and measured interface positions at $\xi$9 and $p$00. The simulation reproduces the overall concave “bowing” of the melting front and the tendency for faster melting near the top, but exhibits a systematic offset, with maximum error near the top of order a few millimeters. The same study concludes that natural convection in the liquid layer induces a clockwise roll, and that values of the mushy-zone constant $p$01 were needed to match the observed sharpness of the water–ice front (Schüller et al., 2018).

The two-dimensional double-MRT LB model of Liu and He was validated against one-dimensional conduction melting in a semi-infinite medium, two-dimensional solidification in a semi-infinite corner, and convection-dominated melting in a square porous cavity. The reported results show excellent agreement with analytical Stefan-problem solutions for $p$02 and interface motion $p$03, a very close match to Lin and Chen (1997) and analytic solutions for the corner-solidification interface at $p$04, and agreement of melting fronts, streamlines, and temperature fields with Beckermann and Viskanta (1988) experiments and earlier numerics at $p$05 (Liu et al., 2015). The three-dimensional DDF-MRT extension reports that analytic velocity and temperature profiles are recovered for mixed convection in a porous channel, that the average Nusselt number $p$06 in a cubical porous cavity matches boundary-element results over $p$07 and $p$08, and that a spatial-convergence study shows second-order accuracy with global $p$09 error scaling as $p$10 (Liu et al., 2016).

For high-enthalpy compositional flow with phase separation, Lipovac et al. emphasize narrow-boiling phenomena and nonlinear solver behavior. On a 2D test of cold $p$11 injection into a hot reservoir, the enthalpy-based flash allowed stable capture of sharp boiling fronts. Embedding the local flash solver every Newton iteration reduced total global Newton steps by up to $p$12 compared to calling flash every 5th iteration, specifically $p$13 versus $p$14 global iterations over 6 months. A local-flash tolerance as loose as $p$15 had negligible impact on global convergence. Using temperature-only flashes $p$16 eventually failed on fine meshes $p$17 because it could not resolve the narrow-boiling layer, whereas $p$18-flashes succeeded uniformly. A more complex viscosity model, Lohrenz–Bray–Clark, roughly doubled the number of global Newton iterations $p$19 and quadrupled local flash iterations, reflecting the extra nonlinearity, but the framework remained robust (Lipovac et al., 3 Dec 2025).

6. Scope, assumptions, and recurrent methodological issues

A recurrent issue in this literature is that “enthalpy-based” does not imply a unique constitutive or numerical strategy. Some models place latent heat in an explicit source term, as in $p$20 for REV-scale LB formulations (Liu et al., 2015) or $p$21 for the three-dimensional DDF-MRT model (Liu et al., 2016). Others absorb latent heat directly into the enthalpy variable so that no explicit extra term is needed once the phase enthalpy is modeled correctly (Lipovac et al., 3 Dec 2025). The review literature makes the same point in algorithmic form by updating $p$22 and $p$23 locally after solving for total enthalpy (He et al., 2018).

A second recurrent issue is interface representation. Fixed-grid enthalpy porosity methods and LB enthalpy methods avoid explicit front-tracking by representing the interface through a liquid fraction on the mesh (Schüller et al., 2018, He et al., 2018). Diffuse-interface phase-field models go further by replacing the sharp interface with an Allen–Cahn layer of width $p$24, while still recovering Stefan-type jump conditions in the sharp-interface limit (Ghosh et al., 2021). This suggests that the central distinction is not whether the interface is tracked geometrically, but how the phase state is encoded and coupled to enthalpy.

The assumptions under which these models are posed also differ materially. REV-scale LB models typically assume incompressible flow in a rigid, isotropic porous medium, local thermal equilibrium between fluid and solid, constant physical properties except in the buoyancy term, and negligible density change during phase change (Liu et al., 2015). The three-dimensional DDF-MRT model likewise assumes LTE in an isotropic, rigid porous medium (Liu et al., 2016). By contrast, the quasistatic thermomechanical model includes visco-elasto-plastic deformation, Preisach hysteresis, and a global existence theorem for the coupled system under Hypothesis 3.1 (Gavioli et al., 2021). The evaporation model introduces explicit vapor transport, Korteweg capillary stress, matched-asymptotic sharp-interface analysis, and periodic homogenization to derive Darcy-scale effective parameters (Ghosh et al., 2021).

A common misconception is that temperature is always the natural thermal state variable for porous-media phase change. The high-enthalpy compositional simulations show a contrary case: temperature-only flashes can fail on fine meshes in narrow-boiling regimes, while $p$25-based flashes succeed uniformly (Lipovac et al., 3 Dec 2025). A plausible implication is that enthalpy is not merely a convenient accounting device for latent heat; in high-enthalpy or compositional settings it can also be the more robust thermodynamic coordinate for coupling local equilibrium to global transport.

Taken together, these works show that the enthalpy–porous media phase-change model is best understood as a broad modeling framework. Its unifying element is the use of enthalpy to encode latent heat within porous-media transport, while its variants differ in phase closure, momentum law, numerical discretization, and physical scope—from water–ice melting with natural convection to multicomponent phase separation, deformable freezing soils, and evaporation with upscaling (Schüller et al., 2018, Lipovac et al., 3 Dec 2025, Gavioli et al., 2021, Ghosh et al., 2021).