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.

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 models^{2}. Secondly, they must compute the flow of ice downslope under its own weight.

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.

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

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* models^{9}. Flowline models simulate ice flow along a single line, normally the centreline of the glacier. Spatially distributed models use the whole domain.

### 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 models have been used very successfully with valley glaciers^{6,10} and ice streams within ice sheets^{11-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 classifications
^{9}, 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).

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.

### 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}.

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 length^{4} (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 surface^{16}. 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 line^{4}.

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 divide^{17}. 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 flow^{18}.

The SSA has also been integrated into SIA models as a sliding law^{19} 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 line^{18}.

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 Sheet^{20,21} and New Zealand^{15} 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 domain^{8}.

### 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 depth^{4,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 coast^{17}.

For this reason, there have been a number of attempts to ‘tile’ or integrate Full Stokes and SIA in different regions of ice sheet models^{20}, 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 use^{22}. 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 interior^{17}. 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 km^{8}. Modellers are also limited by a lack of suitable ice-thickness data at that resolution; for example, BEDMAP2 has 5 km resolution^{23}.

A number of experiments have demonstrated that grid size is an important control on the accuracy of the model^{22,24}. This is especially important when modelling the grounding line. Grounding lines require fine-scale modelling^{25}, 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 conditions^{3}.

### 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 reversible^{3,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 line^{17,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.

## Summary

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

- Ice Sheet Modelling (Huybrechts, 2007) – a highly readable introduction
- Numerical modelling glossary
- Ice sheet models
- Mountain glacier models

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?