Model dynamics

The model dynamics allow the glacier thickness and width to evolve in response to the climate perturbations applied.

function calcthck

This function returns ice thickness along the glacier at the new time including the effects of local mass balance and ice flow. This is done using an implicit technique and is complicated by the need to incorporate the effect of possible branches in the glacier system. The equation to be solved is,

[math]\displaystyle{ \frac{\delta\,s} {\delta\,t}= b - \frac{1} {w} \frac{\delta\,wq} {\delta\,x} }[/math]

[math]\displaystyle{ = b- \frac{\delta\,q} {\delta\,x}- \frac{q}{w} \frac{\delta\,w} {\delta\,x} }[/math]

[math]\displaystyle{ = b+ \frac{\delta} {\delta\,x}D \frac{\delta\,s} {\delta\,x}+ \frac{D}{w} \frac{\delta\,w} {\delta\,x}\frac{\delta\,s} {\delta\,x} }[/math]

[math]\displaystyle{ = b+ \frac{\delta} {\delta\,x}D \frac{\delta\,s} {\delta\,x}+ c \frac{\delta\,s} {\delta\,x} }[/math].

Note that we are not incorporating changes in bedrock elevation so that [math]\displaystyle{ \frac{\delta\,s} {\delta\,t}= \frac{\delta\,H} {\delta\,t} }[/math]. The numerical technique involved is based on Crank-Nicolson time differencing:

[math]\displaystyle{ s_{t+\Delta\,t}=s_{t}+ \Delta\,t \left [{b}\right ]_{t}+ \Delta\,t\alpha\left [{\frac{\delta}{\delta\,x}D\frac{\delta\,s} {\delta\,x}+c\frac{\delta\,s} {\delta\,x}}\right ]_{t}+\Delta\,t \left({1-\alpha}\right)\left [{\frac{\delta}{\delta\,x}D\frac{\delta\,s} {\delta\,x}+c\frac{\delta\,s} {\delta\,x}}\right ]_{t+\Delta\,t} }[/math]

where we follow the suggestion of Hindmarsh () and set [math]\displaystyle{ \alpha=1.5 }[/math]. The calculation is best explained in terms of the matrix equation Ax=b, where the vector x (in vector answ) contains the unknown thickness at the new time level; the vector b (in vector rhsd) contains the known thickness at the old time level, the known mass balance and the known flow at the old time level; and the matrix A contains coefficients relating to the thickness and surface elevation at the new time level (i.e., the finite difference approximations to the two terms within the bracketed quantities on the far right hand side). The algorithm proceeds in five stages.

1. Determine the diffusivities (D) from available thickness and surface elevation at the old time level; we approximate diffusivity at the new time level by the value at the old level. Diffusivity is calculated on a staggered grid halfway between adjacent thickness points. This is done by calling function getdiff and checking to see if the particular point under consideration is affected by a branch.
2. Determine the flow convergence (c) from available width and diffusivity at the old time level; again we approximate convergence at the new time level by the value at the old level. Diffusivity is calculated at the original thickness grid points. This is done by calling function getconv and checking to see if the particular point under consideration is affected by a branch.
3. We are now in a position to fill matrix A. We enter a value into the sparse form of A if the mask variable shows that this is appropriate (see function calcmask). This is done by repeated calls to putpcgc. In the case of a non-branching system, this results in a simple tridiagonal matrix. A great deal of care is required to enter this information correctly for branched systems. The known vector b is filled at the same time using the diffusive and convergence calculated above.
4. With the matrix A and vector b filled, we can obtain x by standard matrix inversion algorithms (in this case SLAP, see slapsolve). We only do this if there is a significant length of glacier etc.
5. Finally, unpack the new surface elevations from the matrix inversion output to the appropriate positions in the surface elevation array (usrf) using information in mask, and then find new thicknesses using bedrock elevation (topg) and glacier length.

subroutine putpcgc

This simple subroutine fills a sparse matrix using the coefficients supplied. The sparse matrix has a three column format (row, column and value) and can be used with the SLAP matrix algebra library. The subroutine checks to ensure no zero values are entered. The sparse matrix itself is held in module pcgdwk. This subroutine is called from calcthck and is called once per matrix element.

subroutine slapsolv

This subroutine employs the SLAP matrix algebra library to solve the matrix inversion problem Ax=b, where A and b are known and x is required (in this case, ice thickness). It expects a sparse form of the matrix to exist in module pcgdwk, and will first convert the anticipated three-column format into SLAP’s own storage format (ds2y). The inversion itself is performed using one of a number of preconditioned conjugate-gradient methods available in SLAP (dslucs). Finally, the convergence of the inversion is checked and the program is stopped if convergence has not been achieved. This subroutine is called from calcthck.

function calcflux

This function returns the total horizontal flux of ice along the glacier. This is determined from,

[math]\displaystyle{ Q^{'}=q^{'}w^{'}=-D^{'}w{'}\frac{\delta\,s^{'}}{\delta\,x^{'}} }[/math]

where standard symbols are defined in Table 1. The flux is scaled and is determined on a staggered grid that has points positioned halfway between adjacent thickness grid points. The algorithm takes into account any branching in the system. If both up and downstream points have no ice then the flux is set to zero.
Note that this calculation is not needed for the main thickness-evolution calculation and is performed solely for diagnostic reasons. This function is called directly from valley and the returned information is held in array, flux.

function calcwdth

This function returns the scaled width of the glacier. This is determined using,
[math]\displaystyle{ w^{'}=w^{'}_{r}\left({\frac{H^{'}}{H^{'}_{r}}}\right)^{p}+ w^{'}_{l} }[/math]

where standard symbols are defined in , and [math]\displaystyle{ w_{r} }[/math] and [math]\displaystyle{ H_{r} }[/math] are (scaled) reference widths and thickness. The parameter p is typically 0.5 and is hard coded in the function, while [math]\displaystyle{ w_{l} }[/math] is a small number incorporated to stop width becoming zero in the event of zero ice thickness (which is important in the ice thickness solver). This function is called from valley and the returned information is held in array wdth.