Deriving the Conservative Enthalpy Equation in Sigma Coordinates
This note details the derivation of a strong-conservation form of the enthalpy equation in terrain-following sigma coordinates. Our objective is to formulate the thermodynamics of an ice sheet in a manner that strictly conserves energy when discretized using a finite-volume method (FVM) on a dynamically deforming grid.
To guarantee energy conservation on a moving grid, the conserved quantity for the finite-volume method is the thickness-weighted enthalpy \(HE\). The resulting conservative equation is:
\[ \frac{\partial (HE)}{\partial t} + \frac{\partial (HEu)}{\partial x} + \frac{\partial (HEv)}{\partial y} + \frac{\partial}{\partial \sigma} \left( E\omega - \frac{\kappa}{H} \frac{\partial E}{\partial \sigma} \right) = \frac{H\Phi}{\rho} - \frac{H \rho_w}{\rho} L \, r_d \, \varpi(E). \]
1. Geometric Motivation: The Move to Sigma Coordinates
Ice sheets present a unique geometric challenge for numerical modeling. They are bounded by a complex, irregular basal topography (\(b\)) and a freely evolving upper surface (\(h_s\)).
If we attempt to model this system on a rigid Cartesian grid, these boundaries will arbitrarily intersect our computational cells. This introduces severe numerical complications, requiring either stair-stepped boundaries or fractional cut-cell tracking. Furthermore, because the most critical thermodynamic processes, such as extreme strain heating and phase transitions, occur in a thin boundary layer near the bed, a Cartesian grid struggles to provide adequate vertical resolution without wasting computational effort on empty space above the ice.
We resolve this by transforming the physical domain into terrain-following sigma coordinates:
\[ \sigma = \frac{z - b}{H} \]
where \(H = h_s - b\) is the ice thickness. This mapping neatly abstracts the complex geometry, transforming the time-varying, irregular ice column into a static computational unit cube where \(\sigma \in [0, 1]\).
2. The Cartesian Baseline for an Incompressible Tracer
To understand why this coordinate transformation complicates the thermodynamics, it helps to first review how a tracer is conserved in a standard Cartesian framework.
The finite-volume method requires governing PDEs to take a flux-divergence form (\(\partial_t U + \nabla \cdot \mathbf{F} = Q\)). This allows the Divergence Theorem to convert volume integrals into surface fluxes, ensuring that energy leaving one cell enters its neighbor exactly.
The advective transport of specific enthalpy \(E\) following the fluid flow is given by the material derivative:
\[ \frac{\partial E}{\partial t} + \mathbf{v} \cdot \nabla E = S \]
Because this equation relies on the spatial gradient of \(E\) rather than the divergence of a flux, it is not inherently conservative. However, glacier ice is treated as an incompressible fluid, meaning the velocity field is divergence-free:
\[ \nabla \cdot \mathbf{v} = 0 \]
We can multiply the continuity equation by \(E\) and add it to the advective equation without changing the physical balance:
\[ \frac{\partial E}{\partial t} + \mathbf{v} \cdot \nabla E + E(\nabla \cdot \mathbf{v}) = S \]
Applying the product rule groups the spatial terms, cleanly yielding the required flux-divergence form:
\[ \frac{\partial E}{\partial t} + \nabla \cdot (\mathbf{v} E) = S \]
In a fixed Cartesian grid, FVM can be applied directly to the specific enthalpy \(E\) because the incompressibility of the fluid guarantees that the transport terms perfectly balance.
3. The Compressible Analogy on a Deforming Grid
When we map our equations into sigma space, the computational grid stretches and compresses over time as the ice thickness \(H\) evolves.
From the perspective of this deforming computational grid, the effective velocity field is no longer divergence-free. The geometric stretching of the grid acts as a local expansion or compression of the computational control volume. If we were to evaluate the advective equation in sigma coordinates and update \(E\) independently of these volume changes, total energy would drift.
This situation is mathematically identical to modeling a variable-density fluid (like a compressible gas) on a fixed grid. In our framework, the ice thickness \(H\) (mass per unit area) plays the exact role of fluid density \(\rho\) (mass per unit volume).
In compressible fluid dynamics, there is a standard, mechanical procedure to recover a strict conservation law for any specific tracer \(\phi\). We shift the state variable from the specific quantity (\(\phi\)) to the density-weighted quantity (\(\rho\phi\)) using the following recipe:
- Take the advective tracer equation and multiply it by the density \(\rho\).
- Take the mass continuity equation and multiply it by the specific tracer \(\phi\).
- Add the two equations together. The reverse product rule will perfectly collapse the advective gradients and divergences into a single conservative flux divergence: \(\partial_t(\rho\phi) + \nabla \cdot (\rho\mathbf{v}\phi)\).
The remainder of this note simply executes this exact mathematical recipe. We will define our advective enthalpy equation, rigorously derive our sigma-space continuity equation, and then blend them together.
4. Governing Continuous Equations
The advective transport of specific enthalpy in a moving sigma-coordinate framework is:
\[ \frac{\partial E}{\partial t} + u \frac{\partial E}{\partial x} + v \frac{\partial E}{\partial y} + \dot{\sigma} \frac{\partial E}{\partial \sigma} = S_{\sigma} \]
The source term \(S_{\sigma}\) contains three processes: vertical thermal diffusion, internal strain heating (\(\Phi\)), and meltwater drainage. Applying the chain rule for vertical derivatives (\(\partial_z = \frac{1}{H}\partial_\sigma\)), the source term becomes:
\[ S_{\sigma} = \frac{1}{H} \frac{\partial}{\partial \sigma} \left( \frac{\kappa}{H} \frac{\partial E}{\partial \sigma} \right) + \frac{\Phi}{\rho} - \frac{\rho_w}{\rho} L \, r_d \, \varpi(E) \]
The drainage term removes enthalpy from temperate ice to represent the gravitational percolation of liquid water out of the ice matrix. The water content \(\varpi\) is defined as the liquid fraction above the pressure melting point:
\[ \varpi(E) = \frac{\max(E - E_{\text{pmp}}, \; 0)}{L} \]
where \(E_{\text{pmp}} = c_i(T_{\text{pmp}} - T_{\text{ref}})\) is the local enthalpy at the pressure melting point, \(L\) is the latent heat of fusion, and \(r_d\) is the drainage rate.
5. Transforming the Continuity Equation
To proceed, we map the Cartesian statement of incompressibility (\(\partial_x u + \partial_y v + \partial_z w = 0\)) into sigma coordinates. Applying the multivariate chain rule to the horizontal derivatives accounts for the slopes of the constant-\(\sigma\) layers:
\[ \left(\frac{\partial u}{\partial x}\right)_z = \left(\frac{\partial u}{\partial x}\right)_\sigma - \frac{1}{H} \left( \frac{\partial b}{\partial x} + \sigma \frac{\partial H}{\partial x} \right) \frac{\partial u}{\partial \sigma} \] \[ \left(\frac{\partial v}{\partial y}\right)_z = \left(\frac{\partial v}{\partial y}\right)_\sigma - \frac{1}{H} \left( \frac{\partial b}{\partial y} + \sigma \frac{\partial H}{\partial y} \right) \frac{\partial v}{\partial \sigma} \]
Substituting these alongside \(\partial_z w = \frac{1}{H} \partial_\sigma w\) into the Cartesian continuity equation, and multiplying the entire expression by \(H\) to clear the fractions, gives:
\[ H\left(\frac{\partial u}{\partial x}\right)_\sigma + H\left(\frac{\partial v}{\partial y}\right)_\sigma + \frac{\partial w}{\partial \sigma} - \frac{\partial u}{\partial \sigma} \left( \frac{\partial b}{\partial x} + \sigma \frac{\partial H}{\partial x} \right) - \frac{\partial v}{\partial \sigma} \left( \frac{\partial b}{\partial y} + \sigma \frac{\partial H}{\partial y} \right) = 0 \]
Using the reverse product rule (\(H \partial_x u = \partial_x(Hu) - u \partial_x H\)), we force the horizontal terms into a flux-divergence form. Manually adding and subtracting the time derivative \(\partial_t H\) completes the mass conservation statement, allowing us to regroup the equation into two brackets:
\[ \begin{aligned} \underbrace{ \left[ \frac{\partial H}{\partial t} + \frac{\partial (Hu)}{\partial x} + \frac{\partial (Hv)}{\partial y} \right] }_{\text{Divergence of Horizontal Mass}} &+ \\ \underbrace{ \left[ \frac{\partial w}{\partial \sigma} - \frac{\partial u}{\partial \sigma} \left( \frac{\partial b}{\partial x} + \sigma \frac{\partial H}{\partial x} \right) - \frac{\partial v}{\partial \sigma} \left( \frac{\partial b}{\partial y} + \sigma \frac{\partial H}{\partial y} \right) - \left( \frac{\partial H}{\partial t} + u \frac{\partial H}{\partial x} + v \frac{\partial H}{\partial y} \right) \right] }_{\text{Geometric Residuals}} &= 0 \end{aligned} \]
This isolates a clean horizontal mass divergence, alongside a residual containing the vertical and geometric cross-terms. Notice that the second bracket is perfectly integrable with respect to \(\sigma\). We define a new variable, \(\omega\), such that its vertical derivative \(\partial_\sigma \omega\) is exactly equal to the residual bracket:
\[ \omega = w - \left( \frac{\partial b}{\partial t} + u \frac{\partial b}{\partial x} + v \frac{\partial b}{\partial y} \right) - \sigma \left( \frac{\partial H}{\partial t} + u \frac{\partial H}{\partial x} + v \frac{\partial H}{\partial y} \right) \]
This substitution elegantly collapses the transformed continuity equation into its final sigma-space form:
\[ \frac{\partial H}{\partial t} + \frac{\partial (Hu)}{\partial x} + \frac{\partial (Hv)}{\partial y} + \frac{\partial \omega}{\partial \sigma} = 0 \]
Physically, \(\omega\) represents the scaled vertical velocity of the ice relative to the deforming computational grid, or \(H\dot{\sigma}\).
6. Transforming the Kinematic Boundary Conditions
Evaluating \(\omega\) at the domain boundaries clarifies this physical interpretation. The Cartesian kinematic boundary conditions for the surface (\(z = h_s = b + H\)) and bed (\(z = b\)) are respectively:
\[ w_s - \left( \frac{\partial h_s}{\partial t} + u_s \frac{\partial h_s}{\partial x} + v_s \frac{\partial h_s}{\partial y} \right) = -\dot{a} \] \[ w_b - \left( \frac{\partial b}{\partial t} + u_b \frac{\partial b}{\partial x} + v_b \frac{\partial b}{\partial y} \right) = -\dot{m} \]
Substituting \(\sigma = 1\) into our definition of \(\omega\) and expanding \(H = h_s - b\) reveals that all basal terms cancel, exactly matching the Cartesian surface boundary condition to yield \(\omega_s = -\dot{a}\). Similarly, substituting \(\sigma = 0\) removes the stretching terms, matching the basal boundary condition to yield \(\omega_b = -\dot{m}\).
7. Synthesis: Executing the Variable-Density Recipe
We now execute the exact mathematical recipe established in Section 3 to merge our two equations.
Step 1: We multiply the advective enthalpy equation by our computational density, \(H\). Noting that the vertical advective derivative is \(\dot{\sigma}\partial_\sigma E\), multiplying by \(H\) converts this to \(\omega \partial_\sigma E\). (Applying this multiplier to the source term \(S_{\sigma}\) also cancels the \(1/H\) coefficient on the vertical diffusion divergence).
\[ H \frac{\partial E}{\partial t} + Hu \frac{\partial E}{\partial x} + Hv \frac{\partial E}{\partial y} + \omega \frac{\partial E}{\partial \sigma} = \frac{\partial}{\partial \sigma} \left( \frac{\kappa}{H} \frac{\partial E}{\partial \sigma} \right) + \frac{H\Phi}{\rho} - \frac{H\rho_w}{\rho} L \, r_d \, \varpi(E) \]
Step 2: We multiply our new sigma-space continuity equation by our specific tracer, \(E\):
\[ E \frac{\partial H}{\partial t} + E \frac{\partial (Hu)}{\partial x} + E \frac{\partial (Hv)}{\partial y} + E \frac{\partial \omega}{\partial \sigma} = 0 \]
Step 3: We add the two equations together. As anticipated, applying the reverse product rule (\(A\,dB + B\,dA = d(AB)\)) perfectly collapses the advective gradients into strict flux divergences. Moving the vertical diffusion divergence to the left side groups it with the vertical advective transport to form a single, comprehensive vertical flux.
This yields the complete conservative equation:
\[ \frac{\partial (HE)}{\partial t} + \frac{\partial (HEu)}{\partial x} + \frac{\partial (HEv)}{\partial y} + \frac{\partial}{\partial \sigma} \left( E\omega - \frac{\kappa}{H} \frac{\partial E}{\partial \sigma} \right) = \frac{H\Phi}{\rho} - \frac{H \rho_w}{\rho} L \, r_d \, \varpi(E) \]
8. Thermodynamic Boundary Conditions in Conservative Form
Because the state variable is now \(HE\), boundary conditions must be applied directly to the vertical fluxes rather than as pointwise constraints on \(E\).
Upper Surface Face (\(\sigma = 1\))
The total vertical flux \(F_\sigma = E\omega - \frac{\kappa}{H}\partial_\sigma E\) has two components at the surface:
- Advective flux: Using \(\omega_s = -\dot{a}\), the advective mass flux correctly injects enthalpy proportional to accumulation, resulting in a flux of \(-E_{\text{surface}} \dot{a}\).
- Diffusive flux: A Dirichlet condition (\(E_s\) prescribed) implicitly determines the diffusive flux. The solver evaluates the gradient between the surface value and the first interior node, allowing the conductive flux to enter the sub-surface naturally through the diffusion stencil.
Lower Bed Face (\(\sigma = 0\))
In a standard advective formulation, the bed requires a Neumann condition for thermal diffusion: \(-\frac{\kappa}{H}(\partial_\sigma E)_0 = q_{\text{geo}} + q_{\text{fric}}\).
In conservative form, the total basal flux \(F_{\sigma, \text{bed}}\) automatically combines this diffusive contribution with the advective enthalpy loss caused by basal melting (\(E_b \omega_b = -E_b \dot{m}\)):
\[ F_{\sigma, \text{bed}} = -E_b \dot{m} + q_{\text{geo}} + q_{\text{fric}} \]
9. Conclusion
By conceptually treating our deforming ice column as a variable-density fluid, we have successfully cast the thermodynamics into a strict flux-divergence form: \[ \partial_t U + \nabla_\sigma \cdot \mathbf{F} = Q \]
The advantage of this variable-density framing becomes apparent when we move from continuous PDEs to discrete code. Because we have mapped the glaciological enthalpy problem onto the exact mathematical scaffolding of a compressible fluid carrying a specific tracer, we do not need to invent new numerical stabilization schemes from scratch. We can directly adopt proven finite-volume methods for existing problems.