Earth deformation models

The option -bed_def [iso, lc, given] (flag bed_deformation­.model) turns one of the three available bed deformation models.

Point-wise isostasy

The first model -bed_def iso, is instantaneous pointwise isostasy. This model assumes that the bed at the starting time is in equilibrium with the load. Then, as the ice geometry evolves, the bed elevation is equal to the starting bed elevation minus a multiple of the increase in ice thickness from the starting time:

\[b(t,x,y) = b(0,x,y) - f \left[H(t,x,y) - H(0,x,y)\right].\]

Here \(f\) is the density of ice divided by the density of the mantle, so its value is determined by the values of bed_deformation­.mantle_density and constants­.ice­.density in the configuration file. For an example and verification, see Test H in Verification.


The second model -bed_def lc is much more physical. It is based on papers by Lingle and Clark [91] and Bueler and others [32]. It generalizes and improves the most widely-used earth deformation model in ice sheet modeling, the flat earth Elastic Lithosphere Relaxing Asthenosphere (ELRA) model [92]. It imposes essentially no computational burden because the Fast Fourier Transform is used to solve the linear differential equation [32]. When using this model in PISM, the rate of bed movement (uplift) and the viscous plate displacement are stored in the PISM output file and then used to initialize the next part of the run. In fact, if gridded “observed” uplift data is available, for instance from a combination of actual point observations and/or paleo ice load modeling, and if that uplift field is put in a NetCDF variable with standard name tendency_of_bedrock_altitude in the input file, then this model will initialize so that it starts with the given uplift rate.

All parameters (except for constants­.ice­.density) controlling the Lingle-Clark model are listed below (they all have the prefix bed.deformation.).


Prefix: bed_deformation.

  1. lc­.elastic_model (yes) Use the elastic part of the Lingle-Clark bed deformation model.

  2. lc­.grid_size_factor (4) The spectral grid size is (Z*(grid.Mx - 1) + 1, Z*(grid.My - 1) + 1) where Z is given by this parameter. See [91], [32].

  3. lc­.update_interval (10 365days) Interval between updates of the Lingle-Clark model

  4. lithosphere_flexural_rigidity (5e+24 Newton meter) lithosphere flexural rigidity used by the bed deformation model. See [91], [32]

  5. mantle_density (3300 kg meter-3) half-space (mantle) density used by the bed deformation model. See [91], [32]

  6. mantle_viscosity (1e+21 Pascal second) half-space (mantle) viscosity used by the bed deformation model. See [91], [32]

Here are minimal example runs to compare these models:

mpiexec -n 4 pismr -eisII A -y 8000 -o
mpiexec -n 4 pismr -eisII A -bed_def iso -y 8000 -o
mpiexec -n 4 pismr -eisII A -bed_def lc -y 8000 -o

Compare the topg, usurf, and dbdt variables in the resulting output files. See also the comparison done in [32].

To include “measured” uplift rates during initialization, use the option -uplift_file (parameter bed_deformation­.bed_uplift_file) to specify the name of the file containing the field dbdt (CF standard name: tendency_of_bedrock_altitude).

Use the option -topg_delta_file (parameter bed_deformation­.bed_topography_delta_file) to apply a correction to the bed topography field read in from an input file. This sets the bed topography \(b\) at the beginning of a run as follows:

(31)\[b = b_{0} + \Delta b.\]

Here \(b_{0}\) is the bed topography (topg) read in from an input file and \(\Delta b\) is the topg_delta field read in from the file specified using this option.

A correction like this can be used to get a bed topography field at the end of a paleo-climate run that is closer to observed present day topography. The correction is computed by performing a “preliminary” run and subtracting modeled bed topography from present day observations. A subsequent run with this correction should produce bed elevations that are closer to observed values.


The variable viscous_bed_displacement does not correspond to any measured physical quantity. Do not even attempt to analyze it without a careful reading of [32].

Trying to provide a “hand-crafted” viscous_bed_displacement field to PISM is not a good idea.

Keep in mind that zero viscous_bed_displacement does not mean that the bed deformation model is in equilibrium.

Given bed deformation history

The last option -bed_def given can be used if a bed deformation history (i.e. bed elevation changes relative to a reference topography) is known from an external solid-Earth model1. This can be useful when running simulations using offline coupling to such a model.

The bed topography \(b\) is set to

\[b(t,x,y) = b_{\text{ref}}(x,y) + \Delta b(t,x,y),\]

which is a time-dependent version of (31).

This class uses two input files:

  1. Reference topography \(b_{\text{ref}}(x,y)\) (variable topg, in meters).

  2. Time-dependent history of bed elevation changes \(\Delta b(t,x,y)\) relative to the reference topography (variable topg_delta, in meters).

Use the following configuration parameters (prefix: bed_deformation.given.) to set them.

  1. file Name of the file containing time-dependent topg_delta.

  2. reference_file Name of the file containing the reference bed topography topg.


It is possible to combine high-resolution reference bed topography with low-spatial-frequency bed elevation changes: both files have to use the same grid projection and cover the modeling domain but they do not have to use the same grid.



E.g. a relative sea level model.

Previous Up Next