Mass conservation — continuity equation

Mass conservation evolves the ice thickness H in time given the depth-averaged horizontal velocity \bar{\mathbf u} = (\bar u, \bar v) from the momentum balance and a set of mass-balance inputs (surface, basal, frontal, calving). The implementation lives in src/physics/mass_conservation.f90, with the advection step delegated to src/physics/solver_advection.f90 and the per-timestep driver in src/yelmo_topography.f90.

Continuity equation

Following Robinson et al. (2020), Eq. (2), the ice thickness obeys the vertically integrated mass-conservation equation

\frac{\partial H}{\partial t} \;=\; -\,\nabla\!\cdot\!\bigl(H\,\bar{\mathbf u}\bigr) \;+\; \dot a \;+\; \dot b_g \;+\; \dot b_f \;-\; \dot c,

with

  • H the ice thickness;
  • \bar{\mathbf u} = (\bar u, \bar v) the depth-averaged horizontal velocity;
  • \dot a the surface mass balance (bnd%smb in the code, units \mathrm{m\,a^{-1}} ice equivalent);
  • \dot b_g the basal mass balance on grounded ice (bnd%bmb_grnd);
  • \dot b_f the basal mass balance on floating ice (bnd%bmb_shlf);
  • \dot c the calving rate (positive for mass loss).

All five forcings on the right-hand side are inputs to the mass-conservation step: Yelmo does not generate them — they are supplied either from a coupled boundary class (bnd), from sub-models (basal hydrology, sub-shelf melt), or from a calving law. The calving laws themselves are documented separately and are not covered here.

The continuity equation is written in flux-divergence form so that mass is conserved by construction up to the boundary fluxes: the divergence operator \nabla\!\cdot(H\,\bar{\mathbf u}) moves mass from cell to cell without creating or destroying it, and all source / sink terms appear explicitly on the right-hand side. This is the property the numerical schemes are required to preserve at the discrete level — see Numerical solution.

Operator-split form used in the code

The Yelmo timestep updates H in an operator-split sequence, applying the advective divergence first and the mass-balance source terms in turn. Schematically, per (sub-)timestep \Delta t:

H^{n+1} \;=\; H^{n} \;+\; \Delta t\,\dot H_\mathrm{dyn} \;+\; \Delta t\,\bigl(\dot a + \dot b_g + \dot b_f - \dot c\bigr),

where \dot H_\mathrm{dyn} = -\nabla\!\cdot(H\,\bar{\mathbf u}) is computed by the advection solver and the mass-balance contributions are added one by one through the helper apply_tendency, which also enforces non-negativity of H and clips mass-loss rates to \dot m \ge -H/\Delta t. (A few additional terms — sub-grid discharge dmb, relaxation mb_relax, and a boundary residual mb_resid — are added in the same loop; they are diagnostic in most configurations and do not change the continuum equation above.) The same advection solver is used for tracers such as the ice area fraction f_ice and the level-set field used by some calving schemes.

The key choice is how the advective tendency \dot H_\mathrm{dyn} is discretised. Yelmo provides several schemes, selected at runtime by the namelist parameter ytopo.solver; the two recommended choices — an explicit donor-cell upwind scheme (expl-upwind) and an implicit upwind scheme solved with LIS (impl-lis) — are described in Numerical solution.

Margin-front mass balance

Following Pollard and DeConto (2012,2016), an ice-margin front melting scheme has been implemented that accounts for the melt rate along the vertical face of ice submerged by seawater.

The frontal mass balance (\dot{f}, m yr^{-1}) is calculated as:

\dot{f} = \dot{b}_{\rm eff} \frac{A_f}{A_{\rm tot}} \theta_f

where \dot{b}_{\rm eff} is the effective basal mass balance (the mean of the basal mass balance calculated for the ice-free neighbors), A_{\rm tot}=\Delta x \Delta x is the horizontal grid area and A_f is the area of the submerged faces (i.e., the sum of the depth of submerged ice for each face of the grid cell adjacent to an ice-free cell – potentially four faces in total). \theta_f=10 is a scaling coefficient that implies the face mass balance should be ~10 times higher than the basal mass balance (Pollard and DeConto, 2016, appendix).