Balancing the books¶

2D diagnostics¶

PISM provides a number of 2D diagnostics to keep track of mass conservation.1

All of them are computed as time-averaged fluxes over requested reporting intervals. Positive values correspond to mass gain.

For ice mass, at every grid point we have

ice_mass_accounting_error = tendency_of_ice_mass -
  (tendency_of_ice_mass_due_to_flow +
   tendency_of_ice_mass_due_to_conservation_error +
   tendency_of_ice_mass_due_to_surface_mass_flux +
   tendency_of_ice_mass_due_to_basal_mass_flux +
   tendency_of_ice_mass_due_to_discharge);

ice_mass_accounting_error@long_name = "ice mass accounting error";

Click here to download this ncap2 script.

All names on the right-hand side above correspond to valid PISM diagnostic quantities. They were generated using run-time option -extra_file etc.; see Spatially-varying diagnostic quantities.

To check that all changes in mass are accounted for, download the script above and run2

ncap2 -v -S ice_mass_accounting_error.txt \
      pism_output.nc mass_accounting_error.nc

The variable ice_mass_accounting_error in mass_accounting_error.nc will contain ice mass accounting errors at each point. All values of this variable should be close to or equal to zero. They are not zero (in general) due to rounding errors, though.

Use a shortcut

pismr -extra_file ex.nc -extra_times N -extra_vars mass_fluxes,...

to save all fluxes needed to “balance the books” in terms of ice mass.

Alternatively, use fluxes in terms of “ice amount” (mass per unit area):

ice_amount_accounting_error = tendency_of_ice_amount -
  (tendency_of_ice_amount_due_to_flow +
   tendency_of_ice_amount_due_to_conservation_error +
   tendency_of_ice_amount_due_to_surface_mass_flux +
   tendency_of_ice_amount_due_to_basal_mass_flux +
   tendency_of_ice_amount_due_to_discharge);

ice_amount_accounting_error@long_name = "ice amount accounting error";

Click here to download this ncap2 script.

To save these, use the shortcut

pismr -extra_file ex.nc -extra_times N -extra_vars amount_fluxes,...

Comments

tendency_of_ice_mass_due_to_flow is the change in ice mass corresponding to flux divergence
tendency_of_ice_mass_due_to_conservation_error is the artificial change in ice mass needed to “balance the books”. It is uniformly zero in most simulations.
tendency_of_ice_mass_due_to_surface_mass_balance is the change due to the surface mass balance; note that this is not the same as the provided SMB: in ablation areas this is the effective mass balance taking into account the amount of ice present.
tendency_of_ice_mass_due_to_basal_mass_balance is the effective change due to basal (grounded and sub-shelf) melting.
tendency_of_ice_mass_due_to_discharge combines changes due to calving and frontal melt.

Scalar diagnostics¶

Diagnostics listed above are also available as scalars, integrated over the whole computational domain. The “integrated” mass accounting error can be computed using the ncap2 script below.

ice_mass_accounting_error = tendency_of_ice_mass -
  (tendency_of_ice_mass_due_to_flow +
   tendency_of_ice_mass_due_to_conservation_error +
   tendency_of_ice_mass_due_to_basal_mass_flux +
   tendency_of_ice_mass_due_to_surface_mass_flux +
   tendency_of_ice_mass_due_to_discharge);

ice_mass_accounting_error@long_name = "ice mass accounting error";

Click here to download this ncap2 script.

All names on the right-hand side above correspond to valid PISM diagnostic quantities. They were generated using run-time option -ts_file etc.; see Scalar diagnostic quantities.

Comments

tendency_of_ice_mass_due_to_flow is the integral of \(-\nabla \cdot Q\) over the computational domain. This should be zero (up to the effect of rounding errors) in simulations that do not use Dirichlet boundary conditions for ice thickness. Prescribing ice thickness creates sources and sinks, and this diagnostic describes their influence.
tendency_of_ice_mass_due_to_conservation_error should be zero (or close to zero) in most simulations

Mass accounting in subglacial hydrology models¶

PISM’s hydrology models provide all the diagnostic fields needed to keep track of changes in subglacial water thickness.

Note

We keep track of \(W_{\text{till}} + W\), i.e. the sum of the effective thickness of subglacial water stored in till and the effective thickness of subglacial water in the transport layer (if applicable).

At every grid point we have

water_mass_accounting_error = tendency_of_subglacial_water_mass -
  (tendency_of_subglacial_water_mass_due_to_input +
   tendency_of_subglacial_water_mass_due_to_flow +
   tendency_of_subglacial_water_mass_due_to_conservation_error +
   tendency_of_subglacial_water_mass_at_grounded_margins +
   tendency_of_subglacial_water_mass_at_grounding_line +
   tendency_of_subglacial_water_mass_at_domain_boundary);

water_mass_accounting_error@long_name = "subglacial water mass accounting error";

Click here to download this ncap2 script.

All names on the right-hand side correspond to valid PISM diagnostic quantities.

Use a shortcut

pismr -extra_file ex.nc -extra_times N -extra_vars hydrology_fluxes,...

to save all diagnostics mentioned above.

See Subglacial hydrology for more information about hydrology models.

Mass accounting in the PDD model¶

PISM’s PDD model provides diagnostics needed to compare computed accumulation, melt, and runoff to the effective mass balance. Use diagnostic quantities surface_accumulation_flux, surface_melt_flux, and surface_runoff_flux (units of mass per area per time) and surface_accumulation_rate, surface_melt_rate, surface_runoff_rate (units of mass per time).

To save all these, use -extra_vars shortcuts pdd_fluxes and pdd_rates.

Mass conservation and “rough” bed topography¶

Jarosch and others [113] show that Mahaffy’s [64] SIA discretization used by PISM suffers from mass conservation errors near sufficiently abrupt changes in bed elevation. (It may overestimate ice fluxes through the boundary of a grid cell and remove more ice than available, producing a negative ice thickness.)

PISM uses a “projection step” to ensure non-negativity of ice thickness \(H\):

(43)¶\[H^{n+1}_{i,j} = \max(\widetilde H^{n+1}_{i,j}, 0),\]

where \(\widetilde H^{n+1}_{i,j}\) is the “tentative” ice thickness at a grid point \((i,j)\) and the time step \(n+1\) computed using an explicit-in-time finite-volume discretization of the mass continuity equation.

This step is performed after computing the change in ice thickness due to flow and before applying top surface and basal mass balance fluxes (i.e. PISM uses operator splitting in its approximation of the mass continuity equation).

Prior to version 2.0.6 PISM fully relied on (43) to maintain non-negativity of ice thickness and tendency_of_ice_mass_due_to_conservation_error reported the rate at which mass is created by the projection step.

The current mass transport scheme includes a flux limiter (see section 3 and the appendix of [114]) that ensures non-negativity of \(\widetilde H^{n+1}_{i,j}\), making the projection step (43) unnecessary.

Note

PISM still performs the projections step to guarantee that \(H \ge 0\) is true even if the flux limiter fails.
See examples/bedrock_step for PISM’s implementation of the “cliff benchmark” described in [113].

1: See Diagnostic quantities for the full list of diagnostics.
2: ncap2 is a part of NCO.