Higher order velocity schemes

From Interactive System for Ice sheet Simulation
Revision as of 22:24, 30 July 2009 by Sprice (Talk | contribs)

Jump to: navigation, search



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:

  • 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, we need higher-order models.
  • Even in some regions of very slow flow horizontal-stress gradients may be 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 strectching dominates).
  • In shallow ice models, we can solve for the velocity field point by point (in map view) because the velocity solution only depends on local variables (slope, thickness) at that point. In general, higher-order models require the simultaneous solution of all velocities in the model domain because the velocity at any one point is coupled (horizontally as well as vertically) to the velocities elsewhere in the domain (the equations are elliptical). This makes higher-order models computationally much more expensive than shallow ice models, requiring more efficient (and complicated) solution techniques.

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).

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).

  • 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 them model. For the purposes of the summer school, we will be able to do some simple model runs 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 over the other. David Pollard of Penn State and Ed Bueler of Univ. of Alaska Fairbanks currently run large-scale implentation of this type of model.

Practical differences between shallow ice models and higher-order models

By "practical differences", we mean how do we deal with solving the momentum equations in each case (the dynamics) and then how do we use that information (the velocity fields) to evolve the ice sheet geometry in time (i.e. how do we solve dH/dt = ...). There are very large differences for how this is done 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 tridiagonal coefficient matrix that is relatively easy to invert (to solve for the velocities). This process is also readily parallelized, as (if necessary) each column of velocities could be solved for on a separate processor. For higher-order models, we cannot do this, since the solution is not only dependent on vertical gradients but on horizontal gradients as well. The velocity at any point depends on non-local information, leading to an elliptic system of equations. Every velocity must be solved for simultaneously with every other velocity. This is a more difficult numerical problem to solve on one processor, and much more difficult to solve on multiple processors. Nevertheless, for large scale applications (e.g. whole-ice sheet models) and for coupling with climate models, efficient solution and parallelization will be essential. 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. For higher-order models we cannot do this. We must solve the dynamic equations for the velocity field first. After this, we must use some other scheme to calculate thickness evolution.

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 being done on Glimmer/CISM has to do with "upgrading" the model so that it can be used effectively and efficiently with higher-order dynamics schemes.


Higher-order Glimmer/CISM

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 one exercise that involves running the higher-order code through a number of tests from ISMIP-HOM and using the automated test suite to look at output for various model setups.

Another exercise will be to use the higher-order solver to obtain a solution for a simplified ice shelf domain, and then to add a new module to Glimmer/CISM for calculating evolution of that domain.