ASIS:Solver flow

Velocity solver overview
(Tony P initial description for distribution to interested parties)

Equations for momentum balance are simplified by neglecting inertial terms so that they reduce to a static force balance. A non-linear constitutive relation for ice (Glen’s flow law) is used to relate stresses to strain rates (and hence velocities). After a fair amount of manipulation involving further simplifications in the vertical force balance (on the basis of aspect ratio, see above) the following coupled 3-d elliptic equations can be derived (similar to a Stokes equation)

$$2\frac{\partial }{\partial x}f\left( 2\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} \right)+\frac{\partial }{\partial y}f\left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)+\frac{\partial }{\partial z}f\frac{\partial u}{\partial z}=-\rho g\frac{\partial s}{\partial x}$$

$$2\frac{\partial }{\partial y}f\left( 2\frac{\partial v}{\partial y}+\frac{\partial u}{\partial x} \right)+\frac{\partial }{\partial x}f\left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)+\frac{\partial }{\partial z}f\frac{\partial v}{\partial z}=-\rho g\frac{\partial s}{\partial y}$$

where u and v are the two components of horizontal velocity, is ice density (a constant) and g gravity, and the effective viscosity f is

$$f=\frac{1}{2}A^{-1/n}\left[ \frac{\partial u}{\partial x}^{2}+\frac{\partial v}{\partial y}^{2}+\left( \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} \right)^{2}+\frac{1}{4}\left( \left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)^{2}+\frac{\partial u}{\partial z}^{2}+\frac{\partial v}{\partial z}^{2} \right) \right]^{\left( 1-n \right)/2n}$$

where A and n are parameters in the constitutive relation (the former is a constant value 3, while the latter varies by several orders of magnitude with ice temperature). These equations are in z, however we use

$$\sigma $$ coordinates in the model. The transformations for the second-order horizontal terms turn out to be fairly involved. The upper boundary condition is zero traction (von Neumann) and the lower (basal) boundary condition can be:

1.	zero traction (floating ice),

2.	zero velocity (Dirichlet, in the case of ice frozen to the substratum), or

3.	mixed (where the base is at melting point and the substratum is deforming; the form of the boundary condition then

represents the rheology of the underlying sediment, see below). For a marine boundary, the lateral boundary condition is derived from a force balance with the weight of displaced ocean water (von Neumann). For a margin on land, the force balance is zero; alternatively velocity can simply be set to zero (Dirichlet).  The force balance can get fairly complicated if the margin is not aligned with the grid, so that in practice the lateral boundary condition is applied only at the boundary of the computational domain and an artificial layer of thin ice is created around the true ice sheet.  This technique simplifies the implementation of the boundary condition and has little effect on the results.  We are, however looking into applying the lateral boundary condition correctly at the true ice sheet margin.

The solver has two loops. The inner one solves the linear equation with f calculated from the previous iteration. This is done separately for u and v, using velocity from the previous iteration for the other component. The code calculates coefficients etc and packs them into a sparse matrix. We currently use the SLAP preconditioned conjugate gradient solver (Greenbaum and Seager, LLNL) to solve the linear system. The outer iteration updates velocity components and calculates the new effective viscosity. The velocities are under-relaxed and we iterate until things converge. As with the temperature equation (see below), we would prefer not to employ multiple grids in the vertical but only for the horizontal dimensions of the problem.

In the event that this 3d system is too demanding, there is a simplification. For slow-flow, ice deformation is dominated by vertical shear and the horizontal terms in (2) and (3) have little effect. In this case, the equations reduce to

$$u(z)=-2\left( \rho g \right)^{n}\left[ \frac{\partial s}{\partial x}^{2}+\frac{\partial s}{\partial y}^{2} \right]^{(n-1)/2}\frac{\partial s}{\partial x}\int\limits_{h}^{z}{A\left( s-z \right)^{n}}dz+u(h)$$

and equivalent expression for v. In the case of fast flow or floating ice, deformation is dominated by extension and compression with little vertical shear so that the vertically-integrated (2d elliptic) equations can be used

$$2\frac{\partial }{\partial x}Hf\left( 2\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} \right)+\frac{\partial }{\partial y}Hf\left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)+\tau _{bx}=-\rho g\frac{\partial s}{\partial x}$$

with equivalent expressions for v and f, where

$$\tau _{bc}$$ is basal traction and u, v and f are now assumed constant with depth.

This second approach has the disadvantage that one has to a priori determine what flow regime is appropriate using some diagnostic such as the presence of water at the bed. In addition, it does not couple fast-flowing ice with adjacent areas of slow flow correctly.