Higher Order Physics

From Interactive System for Ice sheet Simulation
Revision as of 09:18, 4 August 2009 by Sprice (Talk | contribs)

Jump to: navigation, search


Pattyn-Bocek Diagnostic Model

We have integrated into CISM an improved version of Frank Pattyn's incomplete 2nd order model (Pattyn 2003). This model modifies the full Stokes equations by neglecting vertical resistive stress (Pattyn). This assumption that the pressure at any point in the ice is due only to the weight of the ice above it and not due to resistance to motion is called the hydrostatic approximation (cite?). The model also assumes that the vertical gradients in the horizontal velocity field are much greater in magnitude than the horizontal gradients in the vertical velocity field (Pattyn). These assumptions allow the model to be simplified by solving for only two rather than three components of the velocity vector field.

Stress Field

Main Article: Stress Field Equations for Pattyn 2003 Model

Boundary Conditions

Main Article: Boundary Conditions for Pattyn 2003 Model

Software Integration

Main Article: Integration of Pattyn 2003 Model

Payne-Price Diagnostic Model

Stress Field

Main Article: a;slkdjf;lksdjfl;skf

Boundary Conditions

Main Article: asldkjfas;ldkfj

Software Integration

Enabling and using higher order computation

The integration has been carried out in such a way that it should not be difficult to place alternate solver options in CISM. In fact, this use has been anticipated and the structure of options in the configuration file input has been designed accordingly. To enable the computation of higher-order velocities using Pattyn's model, add to the configuration file

diagnostic_scheme=[1 or 2]

Both 1 and 2 call Pattyn's model. The difference is that 1 computes velocities on the unstaggered (ice) grid and averages them onto the staggered (velocity) grid before returning them to CISM. Setting the option to 2 averages the ice geometry onto the velocity grid, and computes velocities directly at the centroids. This was done because there was some question regarding the accuracy of the computation when the geometry gradients had been subject to averaging. This question has since been resolved, but the distinct options remain in case they are of use later.

Because many higher-order intercomparison experiments require only that velocities be computed and not that the ice geometry be altered, an additional option has been introduced to support this use case. Adding


to the configuration file will force CISM to output and exit as soon as velocities have been computed.

Specifying the basal boundary condition

Basal boundary conditions are specified using two options in the [ho_options] section: basal_stress_type and basal_stress_input. If basal_stress_type is not specified or is specified as 0, a linear bed model is used; a value of 1 specifies a plastic bed. The basal_stress_input option specifies how the \beta^2 or \tau_0 field is specified. Possible values are:

0: No basal sliding
1: Reciprocal of existing "soft" field from CISM (for easy SIA integration)
2: Reciprocal of existing "btrc" field from CISM (for easy SIA integration)
3: Read 'beta' field from the input NetCDF file. Sliding will be turned off anywhere this field contains a NaN.
4: Use a slip ratio as described in the ISMIP-HOM-F experiment.

Specifying an initial velocity guess

The success of an iterative scheme depends on the initial guess used to begin the iteration. For most purposes, using the computed shallow-ice (0th order) velocity solution as an initial guess is close enough to the actual solution to lead to convergence. In some cases, however, convergence can be too slow or impossible with this guess; we suspect that this is the case for some ice shelf experiments.

By default, the shallow ice velocity is used as an initial guess. To use a different velocity field as the initial guess, set

guess_specified = 1

and include the initial guess you want to use in the uvelhom and vvelhom fields in the NetCDF input. Note that these fields must be 3D and on the staggered (velocity) grid rather than on the ice grid.

(This hasn't been implemented yet but I'm documenting the requirements here).

Internal Software Concerns

A design goal for the new higher-order modules was that they should be as separate as possible from the older shallow ice velocity computation. To this end, two new modules were introduced. The first, ice3d_lib, contains an analogue to Pattyn's original set of F77 routines (though these have been heavily modified both to better use F90 idioms and to implement additional model requirements). The second file, glide_velo_higher, acts as a facade between the higher-order routines and the rest of CISM. This extra facade module was created for two reasons. First, it separates the higher order code from the rest of CISM and allows the series of calls necessary to perform the computation (and transform data between formats required by Pattyn and formats used by CISM) to appear as a single call to the higher-order solver in the thickness evolution code. Second, it anticipates future expansion by allowing glide_velo_higher to include facades to other higher-order models as they are introduced while ice3d_lib remains strictly a module for Pattyn's model.

The glide_global_type has been modified in two ways by the introduction of higher-order code. First, options relevant to the solver have been introduced into the glide_options substructure. Second, a new substructure, glide_velocity_hom, has been introduced to manage fields specifically related to higher-order computation. This was done both because many of the fields are analogous to fields already present in glide_velocity (such as uvel and vvel) and because doing so avoided polluting existing data structures with the large number of fields required to support higher-order computation. The Glide NetCDF variables list has been modified to expose these variables such that they can be written to the output without affecting the output of shallow ice-related variables (for example, glide_velocty_hom%uvel and glide_velocity_hom%vvel output as "uvelhom" and "vvelhom" respectively).

Prognostic Model

An experimental prognostic scheme using the existing thickness evolution model from Glimmer has been developed. The scheme, similar to the one that Pattyn describes, simply converts the higher-order velocities to vector diffusivities and passes them to the thickness evolution routine. However, it has not been extensively tested on prognostic experiments, and is likely to lack stability.

Alternative prognostic schemes have been suggested, such as the use of higher-order velocities as a sliding law for shallow ice (Bueler) or an incremental remapping scheme (Lipscomb) that has been developed in a stand-alone model and is awaiting integration with CISM.


  • Hooke, R, Principles of Glacier Mechanics. Upper Saddle River, NJ: Prentice Hall, 1998.
  • MacAyeal, et. al., "An ice-shelf model test based on the Ross Ice Shelf, Antarctica", Annals of Glaciology, vol. 23, 1996.
  • Pattyn, F, "A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes", Journal of Geophysical Research, Vol. 108, No. B8, 2003.
  • Versteed, H and Malalasekera, W, An introduction to Computational Fluid Dynamics: The Finite Volume Method, Edinburgh Gate, Essex, England: Pearson Education, 1995.