Well-Balanced Gas-Kinetic Scheme

• Well-balanced GKS is a numerical framework that integrates kinetic flux evolution with source discretization to exactly preserve key steady states like hydrostatic and lake-at-rest equilibria.
• It utilizes BGK-type interface evolution alongside adaptive techniques such as ALE and compact high-order reconstructions tailored for gravitational gas and shallow water systems.
• The scheme demonstrates high-order accuracy and computational efficiency by maintaining equilibrium preservation across multiscale flows, wet/dry fronts, and moving mesh applications.

Well-balanced gas-kinetic scheme (GKS) denotes a family of kinetic-flux numerical formulations in which BGK-type interface evolution, source discretization, and, in recent variants, mesh motion are coupled so that prescribed steady states are preserved exactly or to roundoff. In the current literature, the equilibria of interest include hydrostatic isothermal states for gravitational gas dynamics, lake-at-rest states for shallow water equations (SWEs), and layered rest states for two-layer SWEs. The topic spans symplecticity-preserving BGK methods for gravitational hydrodynamics, unified gas-kinetic schemes (UGKS) for multiscale flow transport under gravity, well-balanced Navier–Stokes GKS with global reconstruction, compact high-order shallow-water GKS on unstructured meshes, and recent arbitrary Lagrangian–Eulerian (ALE) and adaptive-mesh extensions (Luo et al., 2010, Xiao et al., 2016, Chen et al., 2019, Zhao et al., 2023, Liu et al., 19 Aug 2025, Zhao et al., 16 Oct 2025).

1. Definition, equilibrium targets, and problem setting

The phrase well-balanced is not tied to a single PDE or a single equilibrium. In gravitational gas dynamics, Xiao, Cai, and Xu define a well-balanced scheme as one that evolves an isolated gravitational system under any initial condition to an isothermal hydrostatic equilibrium state and keeps that solution thereafter (Xiao et al., 2016). In shallow-water formulations, the target is typically the lake-at-rest solution, characterized by zero velocity and constant free-surface elevation. In the two-layer shallow-water setting, the preserved equilibrium is layered rest, with the upper-layer thickness constant and the lower-layer free surface balanced against bottom topography (Zhao et al., 2023).

System	Equilibrium target	Representative condition
Gravitational gas dynamics	Isothermal hydrostatic equilibrium	$\nabla p = -\rho \nabla \Phi$, $\mathbf{U}=0$
Single-layer SWE	Lake-at-rest	$u=v=0$, $\eta=h+b=\mathrm{const}$
Two-layer SWE	Layered rest	$U_1=V_1=U_2=V_2=0$, $h_2=\mathrm{const}$, $h_1+B=\mathrm{const}$
ALE shallow water on moving mesh	Lake-at-rest on arbitrarily moving meshes	Discrete flux, source, mesh-motion, and bottom-advection terms cancel while preserving GCL

For the single-layer inviscid SWE with bathymetry $b(x,y)$, the governing equations are

$$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$

$$\partial_t (h u) + \partial_x \left(h u^2 + \frac{1}{2} g h^2\right) + \partial_y (h u v) = - g h \partial_x b,$$

$\mathbf{U}=0$0

with lake-at-rest given by $\mathbf{U}=0$1 and $\mathbf{U}=0$2 (Zhao et al., 16 Oct 2025). For gas dynamics under gravity, the hydrostatic balance couples pressure to gravitational potential $\mathbf{U}=0$3 through $\mathbf{U}=0$4, and isothermal equilibria yield exponential density profiles (Xiao et al., 2016, Chen et al., 2019).

The central numerical difficulty is that naive discretizations typically do not preserve these equilibria because flux gradients, source terms, and, on moving meshes, geometry-induced terms do not cancel exactly at the discrete level. Well-balanced GKS addresses this by embedding the equilibrium structure directly into reconstruction variables, kinetic interface solutions, or both (Chen et al., 2019, Liu et al., 19 Aug 2025).

2. Kinetic formulation and interface evolution

A defining feature of GKS is that fluxes are not obtained from an approximate Riemann solver alone, but from moments of a time-dependent distribution function. For gravitational gas dynamics, the BGK-type kinetic equation with external force is

$\mathbf{U}=0$5

where $\mathbf{U}=0$6 and $\mathbf{U}=0$7 is the collision time (Xiao et al., 2016). In the one-dimensional notation used in that work,

$\mathbf{U}=0$8

The corresponding time-evolving interface distribution is obtained from the analytic integral solution along phase-space characteristics:

$\mathbf{U}=0$9

This scale-dependent interface solution is the basis of UGKS multiscale behavior, since the ratio $u=v=0$0 continuously connects free transport and hydrodynamic relaxation (Xiao et al., 2016).

For SWE, the kinetic model is adapted to hydrostatic pressure through a Maxwellian-like equilibrium. In single-layer formulations,

$u=v=0$1

with collision invariants $u=v=0$2 and moments

$u=v=0$3

(Liu et al., 19 Aug 2025). In the ALE compact shallow-water scheme, the same kinetic logic is used with a Maxwellian-like equilibrium

$u=v=0$4

and interface moments yield both fluxes and interface states together with their time derivatives (Zhao et al., 16 Oct 2025).

Several well-balanced GKS variants exploit the fact that kinetic evolution supplies more than a flux. The compact GKS for two-layer SWEs and the ALE compact GKS both use the interface distribution to obtain time-dependent interface macroscopic states, gradients, and temporal derivatives; these are then recycled into source-term discretization, gradient updates, and one-stage or multistage time integration (Zhao et al., 2023, Zhao et al., 16 Oct 2025). In the ALE formulation, face motion enters directly through the relative particle speed, and the face-normal kinetic flux becomes

$u=v=0$5

which is consistent with the hydrodynamic ALE flux $u=v=0$6 (Zhao et al., 16 Oct 2025).

3. Mechanisms used to achieve well-balancedness

The literature contains several distinct, technically precise mechanisms for enforcing discrete balance.

The symplecticity-preserving BGK formulation models the gravitational potential as a piecewise step function with a jump at each cell interface. Particle transmission and reflection across the potential barrier are governed by single-particle energy conservation,

$u=v=0$7

together with Liouville’s theorem and the symplecticity-preserving relation $u=v=0$8. This yields exact high-order moment mappings across the potential jump. That construction is used to prove the necessity of an exact Maxwellian for preserving hydrostatic isothermal states, as well as conservation of total mass and total energy including gravitational potential energy (Luo et al., 2010).

A different route is taken by the well-balanced Navier–Stokes GKS with gravitational potential. There, an auxiliary variable $u=v=0$9 is introduced through

$\eta=h+b=\mathrm{const}$0

so that $\eta=h+b=\mathrm{const}$1 becomes constant at isothermal hydrostatic equilibrium. The momentum equation is reformulated so that, at equilibrium, the numerical momentum flux and the numerical source term each vanish separately rather than merely balancing one another. The resulting well-balanced property therefore relies on global reconstruction in the variables $\eta=h+b=\mathrm{const}$2 and on a flux-source splitting that enforces zero numerical flux and zero numerical source at the target state (Chen et al., 2019).

In shallow-water schemes, the dominant mechanism is hydrostatic reconstruction at the level of free-surface variables and source discretization. On arbitrary quadrilateral meshes with hanging nodes, the well-balanced GKS with space-time adaptive mesh refinement reconstructs the free-surface elevation $\eta=h+b=\mathrm{const}$3 rather than $\eta=h+b=\mathrm{const}$4,

$\eta=h+b=\mathrm{const}$5

and uses a source discretization on two triangular subcells inside each quadrilateral. On hanging edges, exact balancing under lake-at-rest requires two Gauss points rather than a single midpoint, because otherwise coarse-fine pressure-flux contributions become mismatched (Liu et al., 19 Aug 2025).

The two-layer shallow-water compact GKS extends this philosophy to coupled hydrostatic sources. Its target equilibrium satisfies

$\eta=h+b=\mathrm{const}$6

and the source terms contain both bed-slope and interfacial pressure coupling. The method uses water-level reconstruction so that $\eta=h+b=\mathrm{const}$7 and $\eta=h+b=\mathrm{const}$8 at edge quadrature points, causing the discrete flux-gradient and source integrals to cancel within each control volume (Zhao et al., 2023).

On moving meshes, well-balancedness requires one further compatibility relation: bathymetry is stationary in the physical frame but advects relative to the moving mesh according to

$\eta=h+b=\mathrm{const}$9

The ALE compact GKS therefore advances bottom topography with the same space-time discretization as the hydrodynamic variables. For lake-at-rest with piecewise linear $U_1=V_1=U_2=V_2=0$0, the paper proves

$U_1=V_1=U_2=V_2=0$1

for hydrodynamic and source terms, and

$U_1=V_1=U_2=V_2=0$2

so that $U_1=V_1=U_2=V_2=0$3 and the lake-at-rest state is preserved to roundoff on arbitrarily moving meshes (Zhao et al., 16 Oct 2025).

4. High-order formulations, compact reconstruction, ALE geometry, and adaptivity

Well-balancedness in GKS is now tightly coupled to high-order space-time discretization. The two-layer compact GKS on triangular meshes uses a compact Hermite-type reconstruction based on cell averages and cell-averaged gradients. For fourth-order spatial accuracy it reconstructs a complete degree-3 polynomial

$U_1=V_1=U_2=V_2=0$4

and combines it with lower-order polynomials through a compact WENO-like nonlinear blending. The same kinetic evolution supplies the time-dependent interface states needed to update cell-averaged gradients by Gauss’s theorem (Zhao et al., 2023).

The recent ALE shallow-water scheme adopts a nonlinear fourth-order compact reconstruction on unstructured moving meshes. Its cubic polynomial

$U_1=V_1=U_2=V_2=0$5

is obtained from constrained least squares on a compact stencil consisting of face- and vertex-neighbor cells. Because the interface kinetic model provides $U_1=V_1=U_2=V_2=0$6, $U_1=V_1=U_2=V_2=0$7, $U_1=V_1=U_2=V_2=0$8, $U_1=V_1=U_2=V_2=0$9, $h_2=\mathrm{const}$0, and $h_2=\mathrm{const}$1 in one stage, the method attains second-order temporal accuracy with the one-stage update

$h_2=\mathrm{const}$2

This is paired with a direct physical-mesh update, so the method avoids remapping (Zhao et al., 16 Oct 2025).

ALE consistency is not only a matter of replacing $h_2=\mathrm{const}$3 by $h_2=\mathrm{const}$4. The moving-mesh formulation also requires exact accounting of the swept face area during each time step. In the cited ALE compact GKS, straight faces with linearly moving endpoints yield an explicit geometric contribution involving

$h_2=\mathrm{const}$5

in the time derivative of the mesh-motion operator. For a constant state, the total mesh-motion-induced update on a moving edge is exactly equal to the swept area $h_2=\mathrm{const}$6, which is the discrete GCL identity (Zhao et al., 16 Oct 2025).

Adaptive mesh refinement introduces a different geometric complication. The well-balanced shallow-water GKS with space-time adaptive mesh refinement is built on arbitrary quadrilateral meshes with hanging nodes, managed by a quadtree structure via p4est with Morton indexing and strict 2:1 balance. Temporal adaptivity is achieved by level-dependent local time stepping:

$h_2=\mathrm{const}$7

Coarse-fine interface fluxes are stored during subcycling and used for conservative synchronization when the coarse cell reaches the corresponding intermediate time (Liu et al., 19 Aug 2025).

Across these formulations, compactness is repeatedly emphasized. In the SWE literature, compact GKS means that only immediate neighbors are needed for reconstruction, while discontinuous left/right interface states and their spatiotemporal derivatives are produced by the same kinetic evolution. This reduces stencil width on unstructured meshes and directly couples reconstruction, flux evaluation, and source treatment (Zhao et al., 2023, Zhao et al., 16 Oct 2025).

5. Numerical behavior and validated application domains

The empirical record reported in the cited papers covers hydrostatic preservation, small-perturbation propagation, strong discontinuities, wet/dry fronts, moving meshes, and multiscale non-equilibrium transport.

For gravitational gas dynamics, UGKS preserves hydrostatic backgrounds while resolving small perturbations and regime transitions. In the one-dimensional pressure-perturbation test with $h_2=\mathrm{const}$8, $h_2=\mathrm{const}$9, and $h_1+B=\mathrm{const}$0, the method maintains the hydrostatic background and resolves the perturbation. In the gravity-modified Sod problem, simulations at $h_1+B=\mathrm{const}$1, $h_1+B=\mathrm{const}$2, and $h_1+B=\mathrm{const}$3 display the continuous transition from Euler to collisionless Boltzmann behavior. In the lid-driven cavity under gravity, the method reports a non-equilibrium phenomenon in the transition regime: heat flux correlated with gravity, flowing from the colder upper region to the hotter lower region, contrary to Fourier’s law (Xiao et al., 2016).

For Navier–Stokes under gravity, the auxiliary-variable well-balanced GKS preserves hydrostatic solutions exactly in one-dimensional tests with the potentials $h_1+B=\mathrm{const}$4, $h_1+B=\mathrm{const}$5, and $h_1+B=\mathrm{const}$6 on 100 uniform cells up to $h_1+B=\mathrm{const}$7. In long-time viscous runs from a highly non-balanced initial condition, the system converges to zero velocity and constant temperature, with final velocities below $h_1+B=\mathrm{const}$8 at $h_1+B=\mathrm{const}$9. The same formulation accurately resolves small perturbations riding on equilibrium; normalized solutions for perturbation amplitudes $b(x,y)$0 and $b(x,y)$1 collapse, with deviations at $b(x,y)$2 occurring at approximately $b(x,y)$3, consistent with double-precision round-off (Chen et al., 2019).

In shallow-water applications, machine-precision lake-at-rest preservation is reported on both adaptive and moving meshes. With dynamic AMR on arbitrary quadrilateral meshes with hanging nodes, the well-balanced GKS produces $b(x,y)$4–$b(x,y)$5 and $b(x,y)$6, $b(x,y)$7–$b(x,y)$8 at $b(x,y)$9, $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$0, and $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$1. In a two-dimensional dam-break efficiency study, the uniform mesh uses 97,792 cells and 240 s CPU time, whereas AMR uses 18,682 cells and 100 s, and space-time AMR uses 18,682 cells and 60 s, corresponding to speedups of approximately $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$2 and $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$3 (Liu et al., 19 Aug 2025).

On moving unstructured meshes, the ALE compact GKS verifies both discrete GCL and well-balancedness. For lake-at-rest over the linear bottom

$$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$4

errors for $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$5 and momentum stabilize below $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$6 on a moving mesh. For a nonlinear Gaussian hump, the $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$7 and momentum errors stabilize below $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$8, with the slight discrepancy localized near the crest and attributed to the reconstruction of nonlinear bathymetry on moving cells. In a small-perturbation test over non-flat bottom, comparisons with a fixed-mesh reference show discrepancies $$\partial_t h + \partial_x (h u) + \partial_y (h v) = 0,$$9 in $$\partial_t (h u) + \partial_x \left(h u^2 + \frac{1}{2} g h^2\right) + \partial_y (h u v) = - g h \partial_x b,$$0 away from the bottom crest. In an irregular-domain dam-break, the adaptive moving mesh attains accuracy comparable to a refined fixed mesh near the dam region while using about one quarter of the degrees of freedom (Zhao et al., 16 Oct 2025).

The two-layer compact GKS further demonstrates that well-balancedness is compatible with strongly coupled hydrostatic sources, internal-wave dynamics, and wet/dry fronts on triangular meshes. The paper reports well-balanced lake-at-rest errors below $$\partial_t (h u) + \partial_x \left(h u^2 + \frac{1}{2} g h^2\right) + \partial_y (h u v) = - g h \partial_x b,$$1 up to $$\partial_t (h u) + \partial_x \left(h u^2 + \frac{1}{2} g h^2\right) + \partial_y (h u v) = - g h \partial_x b,$$2, accurate resolution of weak and strong discontinuity Riemann problems, good agreement with references in equal-density and light-over-dense dam-breaks, and robust bore capture in irregular wet and dry bed configurations (Zhao et al., 2023).

6. Limitations, common misconceptions, and extensions

A common misconception is to treat well-balancedness as synonymous with general robustness, positivity preservation, or exact shock resolution. The cited literature does not support that equivalence. Several shallow-water papers emphasize robustness near bores and wet/dry fronts but state that positivity-preserving details are not the focus or that no formal positivity-preserving proof is given (Zhao et al., 2023, Liu et al., 19 Aug 2025, Zhao et al., 16 Oct 2025). Conversely, the Navier–Stokes gravity formulation explicitly notes that its momentum discretization is non-conservative because part of the gravitational source is merged into the flux; small defects appear near rarefactions in shock-tube comparisons, and the recommended application range is low-speed continuous flows (Chen et al., 2019).

Another misconception is that well-balancedness is only a stationary-mesh issue. The ALE shallow-water formulation shows that equilibrium preservation on moving meshes additionally requires compatibility with the geometric conservation law and with bottom-topography advection relative to the mesh. Uniform-state preservation is tied to exact swept-area accounting, and lake-at-rest preservation on moving cells requires the discrete mesh-motion and bathymetry updates to preserve free-surface elevation (Zhao et al., 16 Oct 2025).

The present literature also makes clear that the target equilibrium matters. The UGKS construction under gravity targets isothermal hydrostatic equilibrium; extensions to non-isothermal steady states would require modified reconstruction that embeds the appropriate steady relations, such as simultaneous gravitational balance and conductive equilibrium (Xiao et al., 2016). Similarly, the auxiliary-variable Navier–Stokes GKS is designed around the fact that $$\partial_t (h u) + \partial_x \left(h u^2 + \frac{1}{2} g h^2\right) + \partial_y (h u v) = - g h \partial_x b,$$3 becomes constant at isothermal hydrostatic equilibrium, not at arbitrary stratified steady states (Chen et al., 2019).

Methodological limitations are formulation-specific. In the shallow-water AMR scheme, the implementation is second-order, and higher-order extensions would require more complex reconstructions; the dissipation parameters in

$$\partial_t (h u) + \partial_x \left(h u^2 + \frac{1}{2} g h^2\right) + \partial_y (h u v) = - g h \partial_x b,$$4

may require tuning, and extreme wet/dry topographies may remain challenging (Liu et al., 19 Aug 2025). In the ALE compact GKS, residual terms associated with spatial nonuniformity and face-length or normal changes are neglected in $$\partial_t (h u) + \partial_x \left(h u^2 + \frac{1}{2} g h^2\right) + \partial_y (h u v) = - g h \partial_x b,$$5 for simplicity, and for nonlinear bathymetry exact piecewise-linear reconstruction is not maintained after moving-mesh updates, which leads to bounded but nonzero free-surface errors of at most $$\partial_t (h u) + \partial_x \left(h u^2 + \frac{1}{2} g h^2\right) + \partial_y (h u v) = - g h \partial_x b,$$6 in the reported tests (Zhao et al., 16 Oct 2025). In the two-layer compact GKS, the source discretization is second-order in space for efficiency, which limits global convergence order despite fourth-order reconstruction and a two-stage fourth-order temporal update (Zhao et al., 2023).

At the same time, the extension paths are explicit in the cited works. The ALE shallow-water paper states that the methodology generalizes to three-dimensional SWEs and to compressible Euler/Navier–Stokes equations on moving meshes (Zhao et al., 16 Oct 2025). The AMR paper notes that high-order compact GKS ideas can be adapted to SWE and that additional physics such as friction, wind stress, and Coriolis terms follow the same source-in-flux framework (Liu et al., 19 Aug 2025). The UGKS under gravity points toward non-isothermal equilibria, radiative coupling, and multi-species kinetics, while the two-layer compact GKS remarks that multilayer SWEs, friction, Coriolis forcing, adaptive mesh refinement, and implicit time stepping are natural continuations (Xiao et al., 2016, Zhao et al., 2023).

Taken together, the well-balanced GKS literature shows a consistent pattern: equilibrium preservation is achieved not by a single universal trick, but by matching the kinetic interface model, reconstruction variables, source quadrature, and geometric treatment to the target steady state. The resulting schemes retain the characteristic GKS advantage of time-dependent interface evolution while extending that structure to gravitational forcing, non-flat bathymetry, unstructured high-order reconstruction, moving meshes, and space-time adaptivity (Luo et al., 2010, Xiao et al., 2016, Chen et al., 2019, Zhao et al., 2023, Liu et al., 19 Aug 2025, Zhao et al., 16 Oct 2025).