A introduction to the hierarchy of ice-sheet models

Why use ice-sheet models? | What is a model? | How reliable are ice-sheet models? | Model inputs and outputs | Prognostic and diagnostic models | Flowline and spatially distributed models | Modelling ice flow | The role of grid size | Summary | References | Comments |

Why use ice-sheet models?

A glacier or ice-sheet model is a building block, used for experiments to make quantitative predictions about the future of glaciers and ice sheets. They are important, because they can help scientists to constrain future contributions to sea level rise and ocean and atmospheric circulation. They can be stand-alone units, or can be coupled to earth-system, ocean and climate models for a comprehensive assessment of past and future changes.

In 2007, the IPCC stated that,

“Dynamical processes related to ice flow not included in current [computer] models but suggested by recent observations could increase the vulnerability of the ice sheets to warming, increasing future sea level rise. Understanding of these processes is limited and there is no consensus on their magnitude”1.

This means that ice-sheet models prior to 2007 could not replicate some of the complex, non-linear changes in ice-sheet flow that might affect rates and magnitudes of sea level rise, such as ice shelf collapse or marine ice sheet instability. Modelling changes in sea-ice extent has proven especially difficult.

Since 2007, extensive work has been conducted improving the physical basis of these models. Calving processes have been a particular area of focus. Improvements in computer power mean that ever more complex models can be formulated, coupled and run over long timescales.

Alongside these computationally expensive models is a series of smaller, lighter, simpler models that can be applied to specific problems and research questions. They are capable of carrying out many simulations rapidly, allowing a fuller exploration of parameter space and the effect of different variables.

What is a model?

A model is a simplification of reality. It is never a completely true representation of the real system, which is far too complex to replicate in its entirety. It is for this reason that George Box (1987) wrote,

essentially, all models are wrong, but some are useful.”

Because ice-sheet and glacier models are an approximation of reality, we cannot expect them to yield perfect results. However, they can provide insights into likely scenarios, and if they can replicate past change well, can be used to provide estimates of future change.

Numerical ice-sheet models represent physical real-world processes by a series of equations. These equations are then solved to understand how the system will respond under different scenarios.

The process of modelling.

The process of modelling.

Glacier and ice-sheet models essentially combine two fundamental processes; first, they must model the change in mass from snow fall (accumulation) and snow melt (ablation) onto the glacier (the glacier’s mass balance). This is usually calculated through degree day or energy balance models2. Secondly, they must compute the flow of ice downslope under its own weight.

Equilibrium line altitudes in a hypothetical glacier

Equilibrium line altitudes in a hypothetical glacier

The video below explains the process of modelling Rhone Glacier, in Switzerland.

Glacier flow is a continuous process, but this is very difficult to simulate mathematically. Instead, modellers use numerical methods, which solve equations in a series of steps. Models are split up into grid cells; the numerical methods calculate the change in ice thickness through time (using a series of partial differential equations) for one grid cell at a time, and use the solution as an input for the next cell. This is then iteratively repeated at each time step. For example, in the figure below, the change in ice thickness since the last time step is calculated for each cell, one at a time.

Grid cells in a flowline model (cross-sectional view)

Grid cells (50 m resolution) in a flowline model (cross-sectional view). Equations are solved step-wise, one grid cell at a time. The solution for the first grid cell is used as input for the second.

How reliable are ice-sheet models?

Can a numerical model be used to make a useful prediction? Perhaps, if the model physics are a suitable approximation of the real system, and the appropriate values for the input variables have been used. Models are therefore usually verified (the model physics tested against good models or good analytical solutions) (e.g., 3-5) and validated (run against a period for which observations are available, and tested against these observations) (e.g., 6).

Model inputs and outputs

Numerical glacier and ice-sheet models need input data. This includes model constants (such as ice and sea water density, the thermal conductivity of ice, acceleration due to gravity), and environmental variables (such as temperature, precipitation, geothermal heat flux, temperature and precipitation lapse rates, temperature ranges. Other input data may include sea level, initial ice surface elevation and bed topography.

Some model parameters may be unknown, and must therefore be tuned. This means that the model parameters are systematically varied until the model outputs (ice geometry, velocity, length, volume and so on) match observations. This is generally the most time consuming part of any modelling programme, and there are generally only a few solutions that produce a realistic model output.

Model constants and environmental variables.

Example of some model constants and environmental variables.

Once the model has been parameterised, it may be dynamically calibrated. This means that a series of short (100-200 years), time-dependent runs are carried out, forced by observed temperature and precipitation. The transient response of the model is recorded, and parameters are fine-tuned until model outputs match the record of glacier length or volume, including the rate of change over the time period. If the model can replicate the recent glacier fluctuations, you may have more confidence that it is validated and is therefore likely to produce reliable results.

Prognostic and diagnostic models

There are many different kinds of model, which can use a variety of different physical solutions. The most fundamental division in kinds of model is the difference between prognostic and diagnostic models.

Diagnostic models focus on a specific process. They may isolate a part of a glacier or ice sheet (for example, the grounding line or calving front) and explore the workings of that particular system. They consider the physics of a specific process in a schematic manner.

Prognostic models are used for understanding the past, and for predicting the future. They are time-marching, simulating evolution through time. Steady-state prognostic models have stable input parameters and are used to find a specific solution (such as, the geometry of an ice field at the LGM 7). Transient simulations have time-dependent input variables; parameters like temperature and precipitation may change through time. For example, proxy records of air or ocean temperature may be used to drive the model. The scientist is interested in the transient response of the ice mass to these changes (e.g., Ref. 8).

Flowline and spatially distributed models

The next major distinction in numerical models is whether they are spatially distributed (also known as plan view) or flowline models9. Flowline models simulate ice flow along a single line, normally the centreline of the glacier. Spatially distributed models use the whole domain.

Plan view and spatially distributed numerical glacier and ice-sheet models

Plan view and spatially distributed numerical glacier and ice-sheet models

Flowline models

Flowline models simulate ice-flow along a single, one-cell wide centreline. They usually include a width-dependent shape factor that accounts for the valley width, and which accounts for conservation of mass as the valley widens and narrows. For example, compressional flow will occur as the valley narrows, resulting in accelerating ice velocities. The incorporation of a width-dependent variable ensures that this is included in the model’s physics.

Flowline models are often called ‘one-dimensional (1D)’ models as they have one dimension. If they average the processes of ice flow vertically in the ice column, they are ‘depth-integrated’. Some flowline models may be depth-dependent, and the processes of ice flow may vary with temperature in the ice column.

Flowline model domain. After Oerlemans, 1997

Flowline model domain. After Oerlemans, 1997

Flowline models have been used very successfully with valley glaciers6,10 and ice streams within ice sheets11-13, where ice-flow is topographically constrained.

Spatially distributed models

Spatially distributed (also known as plan view) models use the whole domain, and may be preferable when the influence of regional topography may be important (for example, in controlling ice accumulation)14.

There are two main kinds of spatially distributed models: 2D and 3D.

  • 2D models are depth-integrated; they vertically average the processes of ice flow and temperature through the ice column. These kinds of models may also be called ‘isothermal’. These kinds of models have been used very effectively to model ice dynamics, such as in the Southern Alps in New Zealand 7. They are more computationally efficient than the more complex 3D models and often there is little difference in results; this means that they are often used in longer time simulations.
  • 3D models are depth-dependent, and ice temperature is allowed to vary in the ice column, which affects the ease with which the ice deforms. These models are at the top end of model classifications9, and are often used to simulate complex ice dynamics and ice streaming.

Modelling ice flow

Numerical models can use a variety of different solutions, of differing complexities, to simulate ice flow. They have different approximations of the nine principle stress components. The more complex models will include all nine, including longitudinal (stretching and compressional) and transverse stresses (such as drag against the valley sides).

Driving and resisting stresses operating on a block of ice on an inclined slope.

Driving and resisting stresses operating on a block of ice on an inclined slope.

Determining the full stress field and its evolution through time involves using the Full Stokes equations (differential equations that form the cornerstone of fluid dynamics). Most models solve simplified (or reduced) forms of the Stokes equations to approximate the stress field. There is therefore a hierarchy of model complexities, depends on which approximation of the full stokes equations are used.

A hierarchy of ice-sheet models

A hierarchy of ice-sheet models

Shallow Ice Approximation (SIA)

The Shallow Ice Approximation (SIA) neglects longitudinal (along flow stretching and compression) and transverse stresses (lateral drag against slower ice for an ice stream or valley walls for a valley glacier)15, and vertical stress gradients. It is computationally efficient and works well in non-ice streaming regions of an ice sheet, and works well for valley glaciers. It is sometimes called a ‘zero-order’ model.

However, basal sliding is difficult to implement properly, and so the accuracy of the SIA decreases as the contribution of basal slip to ice velocity increases. It is therefore not very good at simulating ice streams and other key regions of an ice sheet, such as ice divides, grounding lines and floating ice shelves (with zero basal drag)16.

Stresses included in Full Stokes versus the Shallow Ice Approximation

Stresses included in Full Stokes versus the Shallow Ice Approximation

The SIA is the simplest approximation of the Full Stokes equations. It assumes that basal shear stress of the grounded ice is completely balanced by the gravitational driving stress. It also assumes a small aspect ratio of vertical to horizontal length4 (which is appropriate for an ice sheet) and a large ratio of vertical to horizontal stresses. Vertical shearing is concentrated close to the bedrock, with almost no vertical shearing near the surface16. It also assumes that the longitudinal derivatives of stress, velocity and temperature are small. This is typical of slow ‘sheet flow’ in the interior of ice sheets, and the SIA has been shown to model this well.

“Shallowness” refers to the small depth to width ratio of ice sheets, which simplifies the equations and reduces computational cost. The SIA is therefore very efficient at modelling ice sheets over long simulations.

Advantages of the SIA include that it is computationally cheap, and so well suited to long runs, and it works well on ‘sheet flow’ in the interior of ice sheets. It is well understood, with a long history of use. Disadvantages of the SIA include that it cannot properly transition between grounded and floating ice, and does not represent ice shelves or ice streams. It is also not appropriate for complex, local changes, such as at the grounding line; it excludes membrane stresses cross the grounding line4.

There are several important models that use the SIA, including:

Shallow-Shelf Approximation (SSA)

The Shallow-Shelf Approximation (SSA) was developed to model ice shelves, where basal shear stress is zero and where longitudinal stresses dominate. The SSA is a vertically integrated (2D) model, with a depth-averaged ice velocity. It is therefore difficult to implement in regions where vertical variations in speed are important, such as across grounding lines or complex ice flow near the ice divide17. The SSA is also a zero-order model.

Hybrid models (HySSA)

Some models combine the SIA and SSA, to couple ice-sheet to ice-shelf flow18.

Hybrid models: SSA and SIA coupling across the grounding line

Hybrid models: SSA and SIA coupling across the grounding line. After Kirchner et al., 2011.

The SSA has also been integrated into SIA models as a sliding law19 for ice streams. The proportion of movement from SIA or SSA varies across the domain, from fully SIA (sheet flow in interior) to a combination (ice streams) to fully SSA (ice shelves). A parameter is adjusted to balance the amount of ice motion due to SSA and SIA, from pure SSA on floating parts of the ice sheet domain, to pure SIA when ice is frozen at the bedrock interface. This enables accurate simulation of ice streaming, where most of the displacement is due to basal slip, and across the grounding line18.

Sliding and deformation in HySSA models. After Kirchner et al., 2011

Sliding and deformation in HySSA models. After Kirchner et al., 2011

The HySSA is computationally efficient, and the combined SIA/SSA stress balance allows higher ice velocities to be simulated. However, the physical assumptions do still break down near the grounding line in these kinds of models. This kind of model has been successfully used to model the Antarctic Ice Sheet20,21 and New Zealand15 ice sheet over several timescales. The Parallel Ice Sheet Model (PISM) is an example of a HySSA model, which superposes the SIA and SSA across the entire model domain8.

Higher-order models

Higher order models include more stress terms, including longitudinal stress gradients, which is significant for ice streams, and vertical shear stresses, which is significant for slower-flowing regions. There are two main kinds; Blatter-Pattyn type models and Second-Order models. Higher-Order models are typically ‘three-dimensional thermomechanical’ models, that account for variations in ice temperature and viscosity with depth4,17. They are an improvement over hybrid models, but are also computationally more expensive and so are less suitable for whole ice-sheet or long simulations.

Full Stokes

Full Stokes models are the gold standard of numerical glacier and ice-sheet models; they account for all nine stress tensors. They are particularly useful over grounding lines, but may not be needed for the interior of the ice sheet, where improvements in modelling of ice flow is minimal but computational cost is far higher. They are challenging to use at continental scale with a fine grid size (resolution). The largest differences between Full Stokes and approximate models are therefore typically found near the coast17.

For this reason, there have been a number of attempts to ‘tile’ or integrate Full Stokes and SIA in different regions of ice sheet models20, with the aim of trying to limit the spatial extent where Full Stokes is applied, speeding up model run times. These hybrid models have been shown to compare well with Full Stokes simulations.

The role of grid size

After choosing a mass balance scheme and method of modelling ice flow, the next most important decision a modeller must make is the grid size to use22. Ideally, the grid size would be approximately one ice-thickness, so around 600 m around the margins of the Greenland ice sheet and around 3 km in the interior17. Unfortunately, we are limited by computer capacity and it is much more common to see grid sizes of 20-30 km for whole ice-sheet modelling. For example, a recent model of the Antarctic Ice Sheet evolution through the last glacial cycle had a grid size of 15 km8. Modellers are also limited by a lack of suitable ice-thickness data at that resolution; for example, BEDMAP2 has 5 km resolution23.

A number of experiments have demonstrated that grid size is an important control on the accuracy of the model22,24. This is especially important when modelling the grounding line. Grounding lines require fine-scale modelling25, but this is very computationally expensive and would result in very long model runs.

A number of approaches have been implemented to solve this problem, including: nested grids, moving grids and adaptive mesh grids.

Nested grids

Nested grids are a fixed-grid approach, but they have a nested finer-resolution grid near the grounding line. They do not track the grounding line explicitly, and the grounding must fall between grid cells. It uses a fixed, Cartesian coordinate system (a square grid). Nested grid approaches are good at simulating steady-state ice sheets, but are less successful at simulating transient conditions3.

Nested grid models

Nested grid models

Moving Grids

Moving grids track the grounding line explicitly, with a stretched coordinate system where the grounding line always coincides with a point. Moving grids are good at capturing grounding line migration, and perturbations are reversible3,22.

Adaptive Mesh Refinement

A relatively new development has been adaptive mesh grids. These grids, based on polygons rather than squares, can be resized as needed, allowing for finer resolution near the grounding line17,25. The grid size is varied in different parts of the model domain, and the mesh evolves as the grounding line moves. Adaptive mesh grids do not necessarily conform to a coordinate system.

Ice sheet modules

Aside from the mechanics of ice flow, many ice-sheet models include different modules designed to deal with specific issues, such as glacio-isostatic rebound, sea ice, ocean and geothermal heat fluxes, basal sliding and calving. Ice sheet models may also be dynamically coupled to earth-system, climate or ocean models. Two-way coupling means that changes in one affect changes in the other, so changes in ice-surface elevation may affect surface air temperatures and atmospheric circulation.


In summary, ice-sheet and glacier models may incorporate a variety of different ways of calculating surface mass balance and ice flow, and deal with complex ice flow, for example in ice streams, at ice divides or at the grounding line, in different ways.

Glacier and ice-sheet models are very powerful tools, essential to the study of our icey world. They allow us to quantitatively link climate change to ice-mass change, and to investigate the effects of surface mass balance changes to ice flow. However, they are a simplification of reality and should not be held as ‘truth’. They all approximate the real world to varying degrees, and all have different problems.

As computer power continues to grow, we can expect to see ever more complex models, running at a finer resolution at a continental scale for ever longer simulations. By coupling with climate models, we will be able to understand dynamic feedbacks between ice-sheets and climate. Overall, the use of ice-sheet and glacier models promises to be vital in glaciology for coming decades.

Further reading


1              IPCC core writing team. Climate change 2007: synthesis report. (IPCC, 2007).

2              Hock, R. Glacier melt: a review of processes and their modelling. Progress in Physical Geography 29, 362-391, (2005).

3              Pattyn, F. et al. Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP. The Cryosphere 6, 573-588, (2012). (open access)

4              Pattyn, F. et al. Grounding-line migration in plan-view marine ice-sheet models: results of the ice2sea MISMIP3d intercomparison. Journal of Glaciology 59, 410-422, (2013).

5              Payne, A. J. et al. Results from the EISMINT model intercomparison: the effects of thermomechanical coupling. Journal of Glaciology 46, 227-238, (2000).

6              Oerlemans, J. A flowline model for Nigardsbreen, Norway: projection of future glacier length based on dynamic calibration with the historic record. Journal of Glaciology 24, 382-389, (1997). (free access)

7              Golledge, N. R. et al. Last Glacial Maximum climate in New Zealand inferred from a modelled Southern Alps icefield. Quaternary Science Reviews 46, 30-45, (2012).

8              Golledge, N. R. et al. Antarctic contribution to meltwater pulse 1A from reduced Southern Ocean overturning. Nat Communications 5, 5107, (2014).

9              Huybrechts, P. Ice sheet modeling. in: B. Riffenburgh (ed): Encyclopedia of the Antarctic, Routledge, New York and London, 514-517, (2007).

10           Davies, B. J. et al. Modelled glacier response to centennial temperature and precipitation trends on the Antarctic Peninsula. Nature Climate Change 4, 993–998, (2014).

11           Gladstone, R. M. et al. Calibrated prediction of Pine Island Glacier retreat during the 21st and 22nd centuries with a coupled flowline model. Earth and Planetary Science Letters 333–334, 191-199, (2012).

12           Jamieson, S. S. R. et al. Understanding controls on rapid ice-stream retreat during the last deglaciation of Marguerite Bay, Antarctica, using a numerical model. Journal of Geophysical Research: Earth Surface 119, 1-17, (2014).

13           Jamieson, S. S. R. et al. Ice-stream stability on a reverse bed slope. Nature Geosci 5, 799-802, (2012).

14           Plummer, M. A. & Phillips, F. M. A 2-D numerical model of snow/ice energy balance and ice flow for paleoclimatic interpretation of glacial geomorphic features. Quaternary Science Reviews 22, 1389-1406, (2003).

15           Adhikari, S. & Marshall, S. J. Parameterization of lateral drag in flowline models of glacier dynamics. Journal of Glaciology 58, 1119-1132, (2012).

16           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: Solid Earth 108, 2382, (2003).

17           Larour, E., Seroussi, H., Morlighem, M. & Rignot, E. Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM). Journal of Geophysical Research: Earth Surface 117, F01022, (2012).

18           Kirchner, N., Hutter, K., Jakobsson, M. & Gyllencreutz, R. Capabilities and limitations of numerical ice sheet models: a discussion for Earth-scientists and modelers. Quaternary Science Reviews 30, 3691-3704, (2011).

19           Bueler, E. & Brown, J. Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model. Journal of Geophysical Research 114, F03008, (2009).

20           Seroussi, H. et al. Coupling ice flow models of varying orders of complexity with the Tiling method. Journal of Glaciology 58, 776-786, (2012).

21           Pollard, D. & DeConto, R. M. Modelling West Antarctic ice sheet growth and collapse through the past five million years. Nature 458, 329-U389, (2009). (Open access PDF)

22           Vieli, A. & Payne, A. J. Assessing the ability of numerical ice sheet models to simulate grounding line migration. Journal of Geophysical Research: Earth Surface 110, F01003, (2005).

23           Fretwell, L. O. et al. Bedmap2: improved ice bed, surface and thickness datasets for Antarctica. The Cryosphere 7, 375-393, (2013). (open access)

24           Durand, G. et al. Full Stokes modeling of marine ice sheets: influence of the grid size. Annals of Glaciology 50, 109-114, (2009).

25           Cornford, S. L. et al. Adaptive mesh, finite volume modeling of marine ice sheets. Journal of Computational Physics 232, 529-549, (2013).

Share this

If you enjoyed this post, please consider subscribing to the RSS feed to have future articles delivered to your feed reader.

One thought on “A introduction to the hierarchy of ice-sheet models

  1. hello, Bethan. the ‘pure’ SIA is based on the assumption that ice is frozen at the bedrock interface. does that means that there is no sliding at base at all for ice sheet flow? I notice in the figure that SIA illustrates both sliding and deformation, although the latter is dominant contribution to motion. If the basal ice adhere closely to the bedrock, how can it flow?

Leave a Reply

Your email address will not be published.