CISM: Solving the equation for thickness evolution

From Interactive System for Ice sheet Simulation
Jump to: navigation, search


Conservation of mass for an ice sheet

For an ice sheet (or any depth-integrated mass flow for that matter), conservation of mass can be expressed by

\frac{\partial H}{\partial t}=-\nabla \cdot \left( UH \right)+\dot{b}

The change in thickness per unit time is given by minus the flux divergence plus a source term (accumulation or ablation). We can expand this out a bit as

\begin{align}\frac{\partial H}{\partial t}=-\left( \frac{\partial \left( UH \right)}{\partial x}+\frac{\partial \left( VH \right)}{\partial y} \right)+\dot{b} \\ 

where {\vec{U}} is the depth-averaged velocity vector (in map plane), and U and V are the depth integrated x and y velocity fields, H is the ice thickness, and {\dot{b}} is the source term, which is the surface mass balance (>0 for accumulation and <0 for ablation). Note that the negative sign in front of the divergence terms (terms in parentheses on the right-hand side) insures sensible behavior. Consider a section of the ice sheet where the thickness is nearly constant and there is no accumulation or ablation. If the velocity gradient along-flow is negative (the ice is slowing down), we expect that to lead to thickening locally (left-hand side of the equation > 0) and vice versa (if the ice is speeding up along flow, that should lead to thinning locally). This is the equation that needs to be solved to calculate ice sheet evolution.

A diffusive approach

For the case of a SIA model, the values of U and V are recast in terms of ice thickness and elevation gradients, in which case the whole problem can be recast as a diffusion equation in ice thickness. In 1d, the equation becomes

\frac{\partial H}{\partial t}=\frac{\partial }{\partial x}\left( D\frac{\partial h}{\partial x} \right)+\dot{b},\quad D=\frac{2A}{n=2}\left( \rho g \right)^{n}H^{n+2}\left| \frac{\partial h}{\partial x} \right|^{n-1}

where D is the non-linear diffusivity (because it depends on the solution to the equation, H), A is the rate factor in Glen's flow lay and n is the power-law exponent, h is the ice surface elevation, and ρ and g are the ice density and the acceleration due to gravity. Importantly, we need only local information in order to solve the above equation. If our velocity solution can not be solved locally, as in the case of higher-order models, we cannot easily use the above formulation to solve ice sheet evolution. In an attempt to use this form and retain a diffusion-solution approach to the problem (diffusion problems generally have nice numerical properties), we could try the following approach (again, in 1d only),

\frac{\partial H}{\partial t}=\frac{\partial }{\partial x}\left( D\frac{\partial h}{\partial x} \right)+\dot{b},\quad D=UH\left( \frac{\partial h}{\partial x} \right)^{-1},

where the U in the expression for D is the depth-integrated velocity field from the higher-order model. This is the approach initially taken by Pattyn (2003). Notice, however, that ice sheet surface slope is in the denominator of the diffusivity here and, as the slopes get smaller and smaller, as they tend to do in regions of fast flow like ice streams and ice shelves, the diffusivity will get larger and larger (approaching infinity as the slope goes to zero). The faster velocities in these regions appear in the numerator of the diffusivity, also making it larger. This is a severe restriction on this approach because the diffusive CFL condition says that

\Delta t<\frac{\left( \Delta x \right)^{2}}{2D},

where Δt is the time step required for stability and Δx is the grid spacing. As the diffusivity goes to infinity, the stable time step goes to zero. In practice, this approach has proven very difficult to use for calculating ice sheet evolution in most of the areas we care about (fast flowing areas with shallow slopes). Thus, it seems that we need an alternate approach.

Advection schemes

The alternate approach is to solve the evolution equation using an advection scheme. Numerically, advection schemes are more problematic than diffusion schemes, but in some cases like this one, they are difficult to get around. The most general advection scheme, and one that we will implement in a related exercise, is a first-order, upwind advection scheme (first-order refers to first-order accurate, as opposed to second-order accurate, which we might prefer).

Most ice sheet models (and fluid dynamic models in general) perform calculations on a "staggered" grid of the type shown below, where velocity components live on one grid and scalar components (e.g. temperatures, pressures, thickness, etc.) live on a grid that is staggered 1/2 grid space in the horizontal dimensions (this leads to numerical advantages - like stability - that we won't get in to much here).

Figure 1: Staggered grid in two dimensions, showing scalars (like thickness, H) at grid cell centers and velocity components, U and V, at grid cell faces. This is a "C" grid. Another staggered-grid possibility is a "B" grid, for which both velocity components live at the grid cell corners. While Glimmer/CISM uses a B grid, averaging of corner velocities to face centers allows one to express them on a C grid if necessary.

Stag grid.png

A "control volume" (or finite volume) approach to solving the equation

\begin{align}\frac{\partial H}{\partial t}=-\left( \frac{\partial \left( UH \right)}{\partial x}+\frac{\partial \left( VH \right)}{\partial y} \right)+\dot{b} \\ 

would be to integrate our equation over the control volume centered on the scalar values. Ignoring source terms for now, and assuming flow only along the x direction (that is, assuming that V~0) we have

\frac{\partial H}{\partial t}=-\frac{1}{\Delta y\Delta x}\int_{s}^{n}{\int_{w}^{e}{\frac{\partial \left( UH \right)}{\partial x}}}dxdy=-\frac{1}{\Delta y\Delta x}\left( HU_{e}-HU_{w} \right)\Delta y.

The "east" and "west" (subscripts e and w) faces of the control volume are shown in the figure below.

Figure 2: Staggered grid in two dimensions, showing locations of interfaces and control volume dimensions. Interfaces e, w, n and s connect the volume centered at P with those volumes to the east, west, north, and south (E, W, N, and S).

Stag gridb.png

In the above equation, we've deliberately left it vague as to which value of H is being advected across the east and west interfaces, into or out of the volume. This is where the term "upwinding" comes in - we choose the scalar value to advect across the interface based on an upwinding criterion. If, for example, the flow at interface e is from left to right (U>0), we would advect the value of H at P OUT of the volume centered at P and INTO the volume centered at E. If, on the other hand, the flow at e was from right to left (U<0), we would advect the value of H at E INTO the volume at P and OUT of the volume at E.

The product of velocity, thickness, and the grid spacing at each interface gives a volume flux in units of cubic meters per year. The sum of the volume fluxes over the total number of faces being considered (two in this case, but four for the 2d case) gives the total volume flux in (total>0) or out (total<0) of the volume. When this number is divided by the area of the volume, the result is the mean thickness change in that volume per unit time required to maintain conservation of mass.

Explicit or Implicit?

Is this an explicit or implicit scheme? If we discretize the right-hand side of the primary equation in time we get

  & \frac{\partial H}{\partial t}=-\nabla \cdot \left( UH \right)+\dot{b} \\
 & \frac{H_{t=1}-H_{t=0}}{\Delta t}\approx -\nabla \cdot \left( UH \right)+\dot{b} \\

This can be rearranged to

H_{t_{t=1}}\approx H_{t=0}+\left[ -\nabla \cdot \left( UH \right)+\dot{b} \right]_{t=0}\Delta t.

The thickness at the new time stop is a function of only variables evaluated at the previous time step. Thus, the scheme is fully explicit and, as such, is subject to the CFL condition for stability.

That's more-or-less all there is to a first-order upwinding scheme, other than extending it to two dimensions for non-zero V. The finite volume formulation guarantees that it will be conservative globally (i.e. all the thickness that's get's moved around on the grid is accounted for at each time step). There are many other "higher-order" advection schemes available, but they are mainly based on the principles outlined here, and mostly rely on fancy corrections to the simple upwind assumption in order to gain more accuracy.

For step-by-step instructions on adding a first-order advection scheme (identical to the one detailed above) to CISM and using it to evolve some of the previous diagnostic test cases, proceed to the third exercise.