Motivation

A significant achievement in the scope of glaciology was Nye's (1952) application of Glen's (1952) isotropic relation between stress and strain of polycrystalline ice. The model is versatile, as it can explain glacier behavior in both map and cross sectional views. Because it sparked the modern age of glaciology, every glaciologist should probably construct and use it at least once. Therefore, the topic for the first introduction to glacier modeling is an obvious choice. An understanding of numerical methods and the ability to intelligently solve problems with a computer are important tools in the twenty-first century, regardless of the field of study. The simplicity of Nye's lamellar flow model leads to a algorithm that is easy to implement and perform for a discretized 2-D glacier.

Reading List

Required

• Nye, JF (1952), The Mechanics of Glacier Flow, Journal of Glaciology, 2, 12, 82-93

Suggested

This content expands upon the work in Cuffey and Patterson's (2010) The Physics of Glaciers and van der Veen's (2013) Fundamentals of Glacier Dynamics.

For strong mathematical derivations, Hutter's (1983) Theoretical glaciology and Greeve and Blatter's (2009) Dynamics of Ice Sheets and Glaciers are very comprehensive references.

List of Tasks

1. Read the recommended literature, the section on theory, and the section on the numerical implementation
2. In a chosen programming language, write a script that can solve the time-evolution of the thickness of a glacier
3. Observe the numerical model under growth conditions leading to the Vialov profile
4. Plot geometry as function of time and surface mass balance
5. Compare plots of 'steady state' geometries to those Vialov profiles
6. Model ice sheet relaxation conditions under zero surface mass balance
7. Generate a hypothesis concerning glacier behavior under the model
8. Write a small report on the testing and results of the hypothesis

List of Symbols

Theoretical Background

This follows the basic derivation given by Nye (1952). There are two governing principles behind the model of lamellar flow in glaciers and ice sheets. The first is that the only source of drag is interactions at the glacier boundaries, such as that between the ice and the bed or the ice and the lateral walls. This means for a two-dimensional cross-section, the ice is in contact with the bed but the upper surface of the ice is free of drag. The same model can describe the observed transverse velocity gradients that characterize glacier flow in map-view if both sides are bound to the side walls. The model described here is that of a 2-D cross-sectional flowline with no basal sliding (\( u_b\)), meaning basal drag (\(\tau_{bx}\)) is equal in magnitude to the driving stress (\(\tau_{dx}\)). \[ \mbox{ \( \tau_{dx} = \tau_{xz} = \tau_{bx} \) when \( u_b = 0\) } \]

The second principle is is the principle of continuity, which states that there is no change in volume through time (accumulation and ablation neglected) and the density (\( \rho \)) does not change. In reality, continuity is based on the the conservation of mass, but because the model presented here concerns a perfectly homogeneous glacier, volume can serve as an appropriate proxy for mass. For any infinitessimal-width vertical column of ice, if the upper surface has a slope that is not directly perpendicular to gravity, there is a flux imbalance between the ice entering one side, and exiting the other, which is described by \( \phi = H U \), with \( H \) the height (thickness) of the column, and \( U \) the vertically-averaged velocity for the column. Because of this imbalance, the rate-change of the ice column thickness with time is non-zero, and is the sum of the flux imbalance and the surface mass balance (\(M\)). \[ {\partial H \over \partial t} = -{d \over dx} \phi + M \] In this derivation, the representation of \( U \) is sought to provide closure to the above continuity equation.

Consider a simple flowline model of a glacier with coordinates \( x \) oriented downflow and parallel to the surface and \( z \) oriented perpendicular to \( x \) with positive-\( z \) upwards. The shear stress at any depth \( z = d \) in a glacier is given by \( \tau_{xz} = \rho g d \sin{\alpha} \). Using the isotropic Nye-Glen creep relation \[ \dot{\varepsilon} = \left({\tau \over B}\right)^n \] it is possible to calculate the velocity at any depth (\( d \)) as \[ u_d = {C \over n+1} (\sin{\alpha})^n d^{n+1} + u_0 \] with \( C=\left({-\rho g / B}\right)^n \). Because \( u_b = 0 \), the velocity difference between the surface and bed is simple to calculate. \[ u_0 = {C \over n+1} (\sin{\alpha})^n H^{n+1} \] Integrating the velocity at depth over the full thickness of the ice column gives the flux (\( \phi \)) of ice through the column. \[ \phi = u_b H + {C\over n+2}(\sin{\alpha})^n H^{n+2} \] The average velocity of ice in the column is \[ U = {2 A H \over n+2}{\tau_{dx}}^n. \] This enables the solution to the continuity equation in terms of \( H \), \( \alpha \), and \( \partial H / \partial t \).

Diffusive Form of Continuity

Combining the continuity equations with those for the driving stress and the average velocity of an ice column, it is possible to represent thickness change as a diffusion equation of the form \[ {\partial H \over \partial t} = {\partial \over \partial x} \left( D {\partial h \over \partial x} \right) + M, \] which has a known solution and is therefore useful for the purposes here. This diffusive form also ensures more stable results during the running of the numerical model that is presented below.

Numerical Solution

While a solution is now possible given the continuous solution above the numerical model presented here is discrete, consisting of only a set of points rather than a continuous surface. Continuous models and equations can be computationally consuming, so models are commonly discretized spatially into a set of points and temporally into time-steps. When the members of the spatial set and the temporal set become infinite in number, the continuous solution is at hand, but even a limited number of points can be employed to good effect. In this section the numerical model implementation of the above solution is given, and certain numerical methods are demonstrated to perform better than others. Here a staggered grid is used, a diagram of which is shown below.

For the diagram of a staggered grid, the thin line shows a continuous ice sheet surface. The dark circles show the discretization of the continuous surface into a set of grid points, here labeled 1 through 6. The smaller blue circles show the location of the off-grid points, which are located at the midpoints of a line connecting adjacent grid points. Notice that at the grid points, ice thickness is known and updated during the model run. The off-grid points are one less in number than the grid points, and have their own set of indices.

Discretization

A staggered grid is such that there are 'grid points' that contain 'true data' on ice thickness and surface elevation and there are 'off-grid points' that contain interpolated data on ice thickness and surface elevation. The off-grid points can occur anywhere between grid points, though it is most simple to consider evenly spaced grid points at steps \( \Delta x \) with off-grid points located at the midpoint between grid points. The grid points are the data set that is updated every time step. Grid points can be thought of as a 1-dimensional column of ice and odff-grid points are used in consideration of surface slope and can be treated as flux gates, where ice enters and exits the 1-D column of ice represented by the grid points. By considering the flux through off-grid points on either side of a grid point, the flux imbalance can be determined, and the rate of thickness change can be calculated just as in the continuous solution.

The model specifications are as follows: This model is composed of \( N \) points evenly distributed in a one-dimensional space oriented parallel to the x-axis. Each point has an associated value for bed elevation (\( b \)) and surface elevation (\( h \)) (or thickness, (\( H \)). The model is time-evolving, with a finite time step \( \Delta t \) over which the model progresses under only the forcings of a uniform surface mass balance rate \( M \) and internal deformation related to weight and gravity.

Algorithm

Depending on the programming language, there may be list-like or vector-like objects that would be well-suited for the one-dimensional coordinate system used in the numerical model. The algorithm presented below assumes that vectors are used to store and access the data in an organized way.

0. Create a vector \( \mathbf{H} \) of length \( N \), which contains real numbers describing the ice thickness in meters, and a vector \( \mathbf{G} \) of length \( N \), which contains real numbers describing the bed elevation at grid-points spaced evenly every \( \Delta x \) meters. \[ \mathbf{H} = \{H_1,H_2,H_3...H_N\} \] \[ \mathbf{G} = \{G_1,G_2,G_3...G_N\} \]
1. For every grid-point, calculate the diffusivity \( D \) and store in a vector \( \mathbf{D} \) of length \( N \). \[ D_{\#} = {2 A\over n+2} (\rho g)^n H_{\#}^{n+2} \left | {\partial h \over \partial x}\right | ^{n-1} \]
2. Calculate the flux of ice \( F \) between adjacent gridpoints and store in a vector \( \mathbf{F} \) of length \( N-1 \). \[ F_{\#'} = {1\over2} (D_{\#} + D_{\#+1}) ({h_{\#+1} - h_{\#} \over \Delta x}) \]
3. Calculate the change in thickness \( C \) of the grid points due to flux imbalance and store in a vector \( \mathbf{C} \) of length \( N \) \[ C_{\#} = {\partial H_{\#} \over \partial t} = {\Delta t \over \Delta x} (F_{\#'} - F_{\#'-1}) \]
4. Update the ice thickness at every grid-point. \[ H_{\#} = H_{\#} + C_{\#} + M \Delta t \]
5. Increment time. \[ t = t + \Delta t \]
6. Loop through steps 1-5 until model end conditions are met.

Model Parameterization

Having established the shell of the program, it is necessary to parameterize the model with numerical values for different rate factors, exponents, and dimensions.

Ice Rheology

Assuming the flow exponent \( n = 3 \), here is a quick reference to some potential values for the prefactor \( A \) in units of kPa^-3 yr^-1 given by Hooke (1981).

Ice Sheet Dimensions

The maximum width of the modern Antarctic ice sheet is around 4600 km, reaching from coast to coast. The maximum width of the modern Greenland ice sheet in the East-West direction is around 1200 km and in the North-South direction is around 2400 km. The maximum width of the Laurentide ice sheet was probably around 5000 km.

Surface Mass Balance

Though there is clearly spatial heterogeneity in accumulation and ablation rates across both modern ice sheets, some first-order values for average surface mass balance can be given. Modeled Antarctic surface mass balance for the grounded portion of the Antarctic ice sheet is an average of around 18 cm ice-equivalent per year (Lenaerts et al, 2012). Modeled Greenlandic surface mass balance for the Greenland ice sheet is an average of around 30 cm ice-equivalent per year (Ettema et al, 2009).

Model Validation

There are certain parameterizations of the model that have unique analytical solutions. By parameterizing the model to the conditions, the correct implementation of the numerical model can be verified. These analytical solutions, their parameterizations, and ways to modify the general algorithm in compliance are presented here.

Vialov Steady-state Solution

The Vialov profile is an analytical solution for the surface profile of an ice sheet with a perfectly horizontal bed subject to a spatially and temporally uniform surface mass balance, and pinned laterally at positions \( x = \pm L \) where the thickness is zero. Under these conditions, there is a steady-state profile defined by Vialov (1958). The continuous function relating ice sheet thickness \( H \), which is equivalent to surface elevation \( h \) for a horizontal bed at elevation 0, at a scalar distance \( x \) from the ice divide to the half-width \( L \) of the ice sheet and the uniform surface mass balance \( M \) is \[ H(x) = \left({2 (n+2)^{1/n} \over \rho g} \left({M \over 2 A}\right)^{1/n} \left(L^{1+1/n} - x^{1+1/n}\right)\right)^{-\left(2+2/n\right)} \]

To compare the numerical solution to the analytical solution, the numerical model must be parameterized with a horizontal bed, and be grown from nothing, meaning that at the initialization of the model, the ice thickness is zero. Therefore, in step 0 of the algorithm, set \( H_\# = 0 \), and \( G_\# = 0 \). Because the Vialov profile requires the edges of the model be pinned at zero, after step 4 of the algorithm, ensure that the \( H_1 = H_N = 0 \). Using these simple modifications, and running the model for long periods to approach the steady state, it is possible to compare near-steady state condition in the numerical model to the analytical solution. With regards to plotting the solutions, it is probably best to center the model domain at \( x=0 \), so that the two solutions are horizontally registered.

Halfar Ice Sheet Relaxation

While the Vialov profile concerns an ice sheet in steady-state, an understanding of the collapse or relaxation of an ice sheet can also be analytically derived. Given an ice sheet like that described by the Vialov profile, with a horizontal bed, if the surface mass balance is turned to zero and the model domain is increased to allow for lateral expansion, then the ice sheet will be in disequilibrium, and the ice sheet will both decrease thickness at the ice divide and widen. Analytical solutions for the condition of ice sheet collapse under nil surface mass balance are given by Halfar (1983) and Nye (2000) for axisymmetric ice sheets. The basic principle of these analytical solutions is that there is a persistent similarity between ice sheet elevation profiles for any time, and that due to the zero mass balance, the ice sheet has a constant volume. In a way, the solutions hinge on using normalized ice sheet thickness \( \hat{H} = H/H_0 \) and radius \(\hat{r} = r/R \) to describe the profile, with the time-dependency describing the mapping between normalized parameter and the steady-state solution.

For an ice sheet with initial ice divide thickness \( \tilde{H}_0 \) and radius \( \tilde{R} \), the thickness at a time \( t \) after the nullification of surface mass balance is given by \[ H_0 (t) = \tilde{H}_0 \left(1+ {t \over T} \right)^{-2/(5n+3)} \] and the corresponding ice sheet radius is given by \[ R(t) = \tilde{R} \left(1+ {t \over T} \right)^{1/(5n+3)} \] with \[ T = {1 \over (5n+3)A_0}\left({2n+1 \over n+1}\right)^n { \tilde{R}^{n+1} \over \tilde{H}_0^{2n+1}} \] \[ A_0 = {2A\over n+2}(\rho g)^n \]

Bueler Thickness-dependent Accumulation

Bueler et al (2005) give an exact, continuous solution for an ice sheet elevation profile under very specific conditions of non-zero surface mass balance. \[ M = {\lambda \over T} H \] with \( \lambda \) a real constant and the starting conditions is not \( t = 0 \), but rather \( t = T \). There are three basic cases for how the choice of \( \lambda \) affects the ice sheet growth or collapse. If \( \lambda = 0 \), then \( M = 0 \), and surface mass balance is nil, which is the condition described in the previous section on the Halfar (1983) solution. If \( \lambda > 0 \), then surface mass balance is positive, and there is accumulation that increases with thickness. If \( \lambda < 0 \), then surface mass balance is negative, and there is ablation that increases with thickness. The reality of the relation between surface mass balance and ice thickness can be questioned, but regardless, the equation for the elevation profiles of ice sheet is \[ H(r,t) = \tilde{H}_0 \left({t \over T}\right)^{-a}\left[1-\left(\left({t\over T}\right)^{-b}{r\over\tilde{R}}\right)^{n/(2n+1)}\right]^{n/(2n+1)} \] with \[ a = {2-(n+1)\lambda \over 5n + 3} \] \[ b = {1+(2n+1) \lambda \over 5n+3} \] The corresponding solution for the thickness at the ice divide can be obtained by setting \( r \rightarrow 0 \) \[ H_0(t) = \tilde{H}_0 \left( {t \over T} \right)^{-a} \] The solution for the ice sheet radius is given by \[ R(t) = \tilde{r} \left( {t\over T} \right)^{b} \]

Future Work

Having tested the model on simple geometries, and simple conditions of surface mass balance, exploration into more dynamic and complex conditions is warranted. There are a number of hypotheses to test with this model, which could concern the effect of bed topography or non-uniform surface mass balance on the surface elevation profile. Because the model is time-evolving, the effects of perturbation at one location can be tracked towards remote reaches of the glacier. It is advised that original hypotheses be created and tested through the use of the model.

References

• Bueler, E, and 4 others (2005), Exact solutions and verification of numerical models for isothermal ice sheets, Journal of Glaciology, 51, 173, 291-306
• Ettema, J, and 6 others (2009), Higher surface mass balance of the Greenland ice sheet revealed by high-resolution climate modeling, Geophysical Research Letters, 36, L12501, doi:10.1029/2009GL038110
• Halfar, P (1983), On the Dynamics of the Ice Sheets 2, Journal of Geophysical Research, 88, C10, 6043-6051
• Hooke, RLB (1981), Flow Law for Polycrystalline Ice in Glaciers: Comparison of Theoretical Predictions, Laboratory Data, and Field Measurements, Reviews of Geophysics and Space Physics, 19, 4, 664-672
• Cuffey, KM & Patterson, WSB (2010), The Physics of Glaciers, Elsevier, isbn:9780123694614
• Glen, JW (1952), Experiments on the Deformation of Ice, Journal of Glaciology, 2, 12, 111-114
• Lenaerts, JTM, and 4 others (2012), A new, high-resolution surface mass balance map of Antarctica (1979-2010) based on regional atmospheric climate modeling, Geophysical Research Letters, 39, L04501, doi:10.1029/2011GL050713
• Nye, JF (1952), The Mechanics of Glacier Flow, Journal of Glaciology, 2, 12, 82-93
• Nye, JF (2000), A flow model for the polar caps of Mars, Journal of Glaciology, 46, 438-444
• van der Veen, CJ (2013), Fundamentals of Glacier Dynamics, CRC Press, isbn:9781439835661
• Vialov, SS (1958), Regularities of glacial shields movement and the theory of plastic viscous flow, International Association of Hydrological Science Publication 47, 266-275