CISM governing equations and numerical solution

From Interactive System for Ice sheet Simulation
Revision as of 14:13, 11 March 2011 by Sprice (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search


Quick links

The important links to other topics that are discussed in more detail below are given here:

  1. Discussion of the higher-order dynamics scheme used in Glimmer/CISM
  2. Discussion of the boundary conditions for the higher-order dynamics scheme in Glimmer/CISM
  3. Discussion of the solution to the higher-order dynamics scheme in Glimmer/CISM
  4. Instructions for running the ISMIP-HOM test suite on output from Glimmer/CISM's higher-order dynamics
  5. Discussion of ice sheet evolution calculations using higher-order dynamics schemes
  6. Instructions for adding a simplified mass-advection scheme to Glimmer/CISM (for evolution calculations)


The main distinction between so-called "higher-order" models and "0-order" (or "shallow ice") models is that higher-order models attempt a closer approximation to solving the non-linear Stokes equations. In general, this usually means incorporating some approximation of horizontal-stress gradients - along-flow stretching or compression and across-flow shearing - in addition to the vertical stress gradients that are accounted for in shallow ice models. This is important for several reasons:

The gravitational stress available to move the ice is the driving stress, indicated in green. Because the ice is assumed to be in equilibrium, the sum of the other stresses is equal to the gravitational driving stress. In the 0-order model (shallow ice approximation) only drag at the glacier base is included and the driving stress is assumed to be balanced by basal drag. In high order models, this restriction is relaxed, and the balance of stresses now extends to lateral and longitudinal stresses. Note, however, that these stresses must be computed from conditions outside of the basic computational cell, and increase the running time of the model.
  • For parts of the ice sheet that we are the most interested in - e.g. ice streams, ice shelves, and other regions of fast flow - horizontal-stress gradients are as or more important than vertical stress gradients. To model the flow in these regions accurately, higher-order models are required.
  • Shallow ice, applied to situations in which there is basal sliding, gives rise to a singularity in the the vertical velocity. Models compute the vertical velocity by integrating incompressablility
\frac{\partial w}{\partial z} = -\frac{\partial u}{\partial x}-\frac{\partial v}{\partial y}
but the horizontal gradients at the bed will depend entirely on the grid spacing where the sliding starts.
  • Incomplete knowledge of the stresses near the grounding line makes it unlikely that shallow ice models will ever be able to accurately simulate grounding line advance and retreat.
  • In some regions of very slow flow, horizontal-stress gradients are important or dominant. For example, ice cores are often recovered at ice divides. Flow modeling is important for interpreting ice core records and using information (such as layer thickness) to infer something about the past flow history in the region. In order to model that flow correctly, one must include horizontal stresses (at the ice divide the surface slope is ~0, in which case vertical stress gradients are also ~0 and horizontal stretching dominates).

The term "higher-order" comes from scaling analyses of the Stokes equations for which a scaling parameter λ=H/L - the ratio of the thickness to the horizontal length scale of interest - is used to assign importance to the various terms. Shallow ice models retain only terms of order 0 while "higher-order" models also retain terms of order 1 (and possibly greater).

The constitutive equations before and after the terms of order higher than one are dropped. Modified terms are indicated in blue. Note that all derivatives of w are gone, eliminating a degree of freedom from the equations.

Available schemes

The most basic and fundamental higher-order scheme is a solution to the full non-linear Stokes equations. Because of the computational burden, however, this is currently not implemented for any practical, 3d, large-scale simulation of ice dynamics (although this is an area of active research, e.g. the ELMER group and those using COMSOL.).

  • Probably the most long-lived higher-order approximation in glaciology is the "shallow-shelf approxmation" (SSA) describing ice shelf flow. It was made popular by Doug MacAyeal in the 80's and 90's. It's main disadvantage is that it is not fully 3d, as it assumes uniform velocity throughout the ice thickness. It is, however, adequate for describing fast flow in many parts of the ice sheets, such as on ice shelves and along some ice streams. In this case, not resolving vertical gradients is a computational advantage.
  • The SSA equations are actually a depth-averaged form of a more general higher-order model, which is commonly referred to as the Blatter-Pattyn model (also see Higher Order Physics for an extensive description). This model has been around since the mid 90's and has become increasingly popular ever since. Two different implementations of the Blatter-Pattyn model are currently incorporated in Glimmer/CISM. These models are both currently being developed and tested and will be freely available to the public in a future release of the model. For the summer school, we will do some simple model runs and experiments using these higher-order schemes.
  • Several "hybrid" schemes exist that are computationally cheaper than the Blatter-Pattyn model. These combine solutions to the shallow ice approximation (for resolving vertical gradients) and the SSA approximation (for resolving horizontal gradients) in some clever way so that a fully 3d solution is obtained. It isn't yet known how well these model solutions compare to fully 3d models, or if one approach (hybrid vs. fully 3d solution) is superior to the other. David Pollard of Penn State and Ed Bueler of Univ. of Alaska Fairbanks currently run large-scale implementation of this type of model (see references below).

Practical differences between shallow ice models and higher-order models

The relationship between several common varieties of ice sheet modesl. Complexity increases along the vertical axis.

By "practical differences", we mean (1) how do we deal with solving the momentum equations in each case (the dynamics, or stress balance) and (2) how do we use the relevant information we derive from that solution (the kinematics, or velocity fields) to evolve the ice sheet geometry in time? There are large differences in how both of these issues are handled - in shallow ice models versus in higher-order models - for two main reasons:

  • The numerical solution of the dynamical equations is fundamentally different in each case. For the shallow ice case, we need only local information (slope and thickness) to solve for the velocity as a function of depth in a single column of ice. We do this pointwise for every location on our model domain (in map view), which is a relatively easy numerical problem; each column of velocities leads to a banded coefficient matrix that is relatively easy to invert (to solve for the velocities). This problem is also what we might call "embarrassingly parallel". That is, very easy to alter for solving efficiently on multi-processor machines. Each column of unknown velocities results in its own tridiagonal matrix, which could be inverted for on its own processor. For higher-order models, however, we cannot do this since the solution at any point also depends on the solution at neighboring points (in map plane). The velocity at any point depends on non-local information, leading to an elliptic system of equations, and every velocity must be solved for simultaneously with every other velocity. The result is a much larger system of equations to solve, which is a more difficult numerical problem to solve on one processor and a much more difficult problem to solve on multiple processors. Because large-scale applications of higher-order models (e.g. whole-ice sheet models and coupling with climate models) will require efficient solution and parallelization techniques, this is a very active area of current research.
  • The equations governing dynamics AND evolution in a shallow ice model can be recast together as a single, non-linear, diffusion equation for ice thickness. A single system of equations is solved to calculate the velocity field and evolve the ice sheet geometry. For higher-order models, we must first solve the momentum balance equations to obtain the velocity field. Then, we need to use some other calculation scheme to evolve the ice thickness.

Both of these differences mean that a model based on shallow ice physics may be built in a fundamentally different way than one based on higher-order physics. Most of the work on Glimmer/CISM over the last ~2 years has had to do with "upgrading" the model so that it can be used effectively and efficiently with higher-order dynamics schemes.

The equations describing a flow that is vertically integrated. The equations are

\frac{\partial}{\partial x}\left ( 2 \eta H 
\left(2\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right)
+\frac{\partial}{\partial y}\left(\eta H\left(
\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
=\rho_w gH \frac{\partial s}{\partial x}

\frac{\partial}{\partial y}\left ( 2 \eta H 
\left(2\frac{\partial v}{\partial y}+\frac{\partial u}{\partial x}\right)\right)
+\frac{\partial}{\partial x}\left(\eta H\left(
\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
=\rho_w gH \frac{\partial s}{\partial y}

The equations describing a flow that is first order are

\frac{\partial}{\partial x}\left ( 2 \eta  
\left(2\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right)
+\frac{\partial}{\partial y}\left(\eta \left(
\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
+\frac{\partial}{\partial z}\left(\eta \frac{\partial u}{\partial z}\right)
=\rho_w g \frac{\partial s}{\partial x}

\frac{\partial}{\partial y}\left ( 2 \eta 
\left(2\frac{\partial v}{\partial y}+\frac{\partial u}{\partial x}\right)\right)
+\frac{\partial}{\partial x}\left(\eta \left(
\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
+\frac{\partial}{\partial z}\left(\eta \frac{\partial v}{\partial z}\right)
=\rho_w g \frac{\partial s}{\partial y}

There are three differences that you should note

  1. The vertically integrated model includes H, or ice thickness in each term. This is a reflection of the integration and does not appear in the first order equations.
  2. Accounting for the thickness not appearing, the only other difference is the presence of a vertical diffusion of horizontal velocities. This is the the third term on the left in the above equations.
  3. The first order equations must be solved on each of a set of horizontal "layers". Layers communicate with each other through the diffusion term.

Both sets of equations are non-linear elliptical equations. Much of the same "technology" can be used solve them. Naively, we would construct a set of horizontal layers that intersect the ice sheets. What's wrong with that? How can we address this problem?

Higher-order Glimmer/CISM

Balance velocities for the Greenland ice sheet (left) from Jonathan Bamber and modeled velocities from higher-order Glimmer/CISM (right) with basal sliding tuned to match balance velocities.

As mentioned above, there are currently two higher-order dynamics (Blatter-Pattyn) schemes being implemented and tested in Glimmer/CISM. They are discussed in varying amounts of detail here (Univ. of Montana implementation of Pattyn's model by Jesse Johnson and Tim Bocek) and here (joint U.K./U.S. effort by Tony Payne and Steve Price). Boundary conditions for the two models are discussed here and here and the numerical solution is discussed here and here.

Relatively recently, Frank Pattyn organized an extremely useful higher-order model intercomparison, ISMIP-HOM, which compared the results from numerous higher-order and full Stokes models in an attempt to define a set of "benchmark" experiments for higher-order models. The results of that study are published in Pattyn et. al (2008). A nice feature of Glimmer/CISM is the ability to run the higher-order codes through an automated test suite, which generates nice figures for comparing model output to the ISIMIP-HOM benchmarks.

Exercises using higher-order Glimmer/CISM

We will do two exercise using higher-order Glimmer/CISM:

  1. Run the higher-order code through a number of tests from ISMIP-HOM and use the automated test suite to look at output for various model setups.
  2. Use the higher-order solver to obtain a solution for a simplified ice shelf domain. Then, add a new module to Glimmer/CISM in order to calculate the evolution of that domain.


  • Bueler, E. and J. Brown. Shallow shelf approximation as a "sliding law" in a thermomechanically coupled ice sheet model. J. Geophys. Res., F03008, doi:10.1029/2008JF001179, 2009.
  • Pollard, D. and R.M. DeConto. Modelling West Antarctic ice sheet growth and collapse through the past five million years. Nature, 458, doi:10.1038/nature07809, 2009.
  • Pattyn, F. A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes, J. Geophys. Res., 108(B8), 2003.