# Analytical solutions for the effect of topography, accumulation rate and lateral flow divergence on isochrone layer geometry

F. Parrenin, R. Hindmarsh, F. Remy

► To cite this version:

F. Parrenin, R. Hindmarsh, F. Remy. Analytical solutions for the effect of topography, accumulation rate and lateral flow divergence on isochrone layer geometry. Journal of Glaciology, 2006, 52 (177), pp.191 à 201. 10.3189/172756506781828728. insu-00375735

HAL Id: insu-00375735

https://insu.hal.science/insu-00375735

Submitted on 23 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analytical solutions for the effect of topography, accumulation rate and lateral flow divergence on isochrone layer geometry

F. PARRENIN,¹,² R.C.A. HINDMARSH,³ F. RÉMY¹

¹Legos, 18 av. Edouard Belin, 31055 Toulouse Cedex, France

²Laboratoire de Glaciologie et Géophysique de l'Environnement (CNRS-UJF), 54 rue Molière, BP 96, 38402 Saint-Martin-d'Hères Cedex, France

E-mail: parrenin@ujf-grenoble.fr

³British Antarctic Survey, Natural Environment Research Council, Madingley Road, Cambridge CB3 0ET, UK

**ABSTRACT.** The effect of spatial variations in ice thickness, accumulation rate and lateral flow divergence on radar-detected isochrone geometry in ice sheets is computed using an analytical method, under assumptions of a steady-state ice-sheet geometry, a steady-state accumulation pattern and a horizontally uniform velocity shape function. By using a new coordinate transform, we show that the slope of the isochrones (with a normalized vertical coordinate) depends on three terms: a principal term which determines the sign of the slope, and two scale factors which can modify only the amplitude of the slope. The principal term depends only on a local characteristic time (ice thickness divided by accumulation rate minus melting rate) between the initial and final positions of the ice particle. For plug flow, only the initial and final values have an influence. Further applications are a demonstration of how the vertical velocity profile can be deduced from sharp changes in isochrone slopes induced by abrupt steps in bedrock or mass balance along the ice flow. We also demonstrate ways the new coordinate system may be used to test the accuracy of numerical flow models.

## NOTATION

Symbols used in this paper, listed in the order in which they are defined:

<table><tr><td>Y(x)</td><td>Flow-tube width</td></tr><tr><td>B(x)</td><td>Bedrock elevation</td></tr><tr><td>S(x)</td><td>Surface elevation</td></tr><tr><td>H(x)</td><td>Total ice thickness</td></tr><tr><td>a(x)</td><td>Surface accumulation rate</td></tr><tr><td>m(x)</td><td>Basal melting rate</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>x</td><td>Horizontal coordinate</td></tr><tr><td>z</td><td>Vertical coordinate</td></tr><tr><td>u<sub>x</sub>(x, z)</td><td>Horizontal velocity</td></tr><tr><td>u<sub>z</sub>(x, z)</td><td>Vertical velocity</td></tr><tr><td>q(x, z)</td><td>Partial horizontal flux (Equation (1))</td></tr><tr><td>Q(x)</td><td>Total horizontal flux</td></tr><tr><td>ω(ζ)</td><td>Flux shape function (Equation (2))</td></tr><tr><td>ζ</td><td>Normalized vertical coordinate</td></tr><tr><td>U&#773;<sub>x</sub>(x)</td><td>Average horizontal velocity</td></tr><tr><td>ω'(ζ)</td><td>Horizontal velocity shape function</td></tr><tr><td>s = ω'(ζ = 0)</td><td>Basal sliding rate</td></tr><tr><td>u<sub>ζ</sub>(x, ζ) = dζ/dt</td><td>Vertical velocity in ζ coordinate (Equations (17) and (19))</td></tr><tr><td>ā(x) = a(x)/H(x)</td><td>Relative surface accumulation rate</td></tr><tr><td>ṁ(x) = m(x)/H(x)</td><td>Relative basal melting rate</td></tr><tr><td>μ(x) = m(x)/(a(x) - m(x))</td><td>Ratio relating basal melting to mass balance (Equation (18))</td></tr><tr><td>T(x)</td><td>Local characteristic time</td></tr><tr><td>π(x) = ln(Q(x)/Q<sub>ref</sub>)</td><td>Logarithmic-flux horizontal coordinate (Equation (20))</td></tr><tr><td>θ(ζ) = ln[(ω(ζ) + μ)/(1 + μ)]</td><td>Logarithmic-flux vertical coordinate (Equation (21))</td></tr><tr><td>u<sub>π</sub>(π, θ) = dπ/dt</td><td>Horizontal velocity in π coordinate (Equation (23))</td></tr><tr><td>ψ(θ)</td><td>Time shape function (Equation (24))</td></tr><tr><td>u<sub>θ</sub>(π, θ) = dθ/dt</td><td>Vertical velocity in θ coordinate (Equation (26))</td></tr><tr><td>χ</td><td>Age</td></tr><tr><td>β</td><td>Non-plug-flow parameter (Equations (42) and (52))</td></tr><tr><td>α</td><td>Slope-sign parameter (Equations (46) and (53))</td></tr><tr><td>λ</td><td>Horizontal amplification term (Equation (51))</td></tr><tr><td>p</td><td>Deformation exponent</td></tr></table>
# 1. INTRODUCTION

Radio-echo sounding of the cold Greenland and Antarctic ice sheets dates back to the 1970s. Such observations not only provide the ice thickness and basal conditions but also show internal layering. Internal layers can be traced over very long distances and are commonly regarded as isochronal layers, i.e. layers with the same age (Bailey and others, 1964; Robin and others, 1969). They are flow-markers, and potentially provide significant information about ice rheology and dynamics, ice boundary conditions such as accumulation rate or basal sliding, ice flow divergence or convergence, the effect of bedrock irregularities, etc. They also provide a strong constraint on, and a test of, ice-sheet models.

An immediate application was the qualitative description of the actual or past ice flow. Robin and others (1977) first pointed out that internal layers often do not closely follow the bedrock, and suggested that ice flow must be diverted by subglacial valleys. The discovery of past local disruption of layers observed in Kamb Ice Stream (Ice Stream C), West Antarctica, allowed Jacobel and others (1996) to conclude that there had been an abrupt stoppage of the ice stream a few centuries ago. Ng and Conway (2004) have deduced that the stagnant Kamb Ice Stream flowed at velocities of over 350 m a⁻¹, from the internal layer stratigraphy.

A second application has been the stratigraphic correlation of ice cores, allowing the comparison of their respective chronologies: Jacobel and Hodge (1995) linked GRIP and GISP2 in Greenland, while Siegert and others (1998) linked Vostok and Dome C in East Antarctica.

A third application has been the quantitative reconstruction of the parameters controlling ice flow. Basal layers observed above Vostok subglacial lake, West Antarctica (Bell and others, 2002) and Ellsworth subglacial lake (Siegert and others, 2006), allowed evaluation of the basal melting and refreezing. The varying depth of shallow isochrones allowed reconstruction of the spatial variation of the accumulation in various areas. For the near-surface isochrones, ice-equivalent layer depth can be viewed as directly proportional to local accumulation rate (e.g. Pinglot and others, 2001). For slightly deeper layers, one needs to correct for total strain rate with a one-dimensional (1-D) model (Vaughan and others, 1999; Fahnestock and others, 2001), and also sometimes for horizontal advection of the ice with two- or three-dimensional models (Nereson and others, 2000; Baldwin and others, 2003; Siegert and others, 2003).

A fourth application has been the extraction of past divide location histories exploiting the Raymond effect (Raymond, 1983), where the non-linear rheology of ice gives rise to slower vertical velocities at the divide. These result in anticlines in the radar echo layers, and are particularly convincing evidence that the layers are isochrones (Conway and others, 1999).

Theoretical investigations into the relation between isochrone characteristics and ice flow boundary conditions or dynamics have been carried out for three decades (Weertman, 1976; Whillans, 1976). Quantitative studies have tended to focus on one aspect of the processes controlling isochrone geometry. This paper is particularly concerned with investigating how the interaction of accumulation rate, ice-thickness changes and deformation

style produce isochrone geometry. This is done with the assumptions of a steady-state ice-sheet geometry and accumulation pattern, and of spatially uniform flow conditions (i.e. spatially uniform velocity shape functions). These assumptions permit development of an analytical expression that relates isochrone slope to spatial variations in accumulation rate, ice thickness and lateral flow divergence. We illustrate our model with a few instructive examples. Such an analytical model shows the influence of each parameter (accumulation, ice thickness, divergence) on isochrones, and gives information about the general nature of the influence. This is complementary to numerical models, which are not necessarily limited by such simple physical assumptions, but which give information about isolated points in parameter space.

# 2. THE (π, θ)-COORDINATE SYSTEM

We consider a flow-tube of an ice sheet in steady state. The use of this assumption is discussed in section 5.2. Time is represented by $t$, and we write the equation in $(x, z)$ coordinates, where $x$, the horizontal coordinate, is the distance from the ice divide and $z$ is the vertical coordinate. We suppose that the direction of the flow does not depend on the vertical coordinate, and we represent the horizontal flow divergence by the varying flow-tube width, $Y(x)$. The ice-sheet geometry is given by $B(x)$ the bedrock elevation, $S(x)$ the surface elevation and $H(x) = S(x) - B(x)$, the total ice thickness. We suppose that snow densifies instantaneously. Let $a(x)$ be the surface accumulation of ice and let $m(x)$ be the basal melting rate at the ice-bedrock interface. We denote $u_x(x, z)$ the horizontal velocity and $u_z(x \pm z)$ the vertical velocity of the ice particles. We further define $q(x, z)$, the partial horizontal flux in the flow-tube, as the flux transported below level $z$:

$$q(x, z) = Y(x) \int_{B}^{z} u_x(x, z') \, dz', \quad (1)$$

with $Q(x) = q(x, S)$ being the total horizontal flux.

We now define the flux shape function, $\omega$, as the function defined on [0; 1] such that the partial flux $q$ is given by:

$$q(x, z) = Q(x)\omega(\zeta), \quad (2)$$

where $\zeta = (z - B)/H$ is the normalized vertical coordinate. Such shape functions have been used previously in glaciology (Lliboutry, 1979; Reeh, 1988; Ritz, 1989). We assume here that $\omega$ is a strictly increasing function (i.e. no reverse flow), and that it is the same function all along the flow-tube (i.e. the flux shape function is uniform). The use of this assumption is discussed in section 5.2.

From Equation (1), the horizontal velocity $u_x$ is given by $u_x(x, z) = (1/Y)(\partial q/\partial z)$, leading to:

$$u_x(x, z) = \bar{U}_x(x)\omega'(\zeta), \quad (3)$$

where $\bar{U}_x = Q/YH$ is the average horizontal velocity and $\omega'$ is the derivative of $\omega$ with respect to $\zeta$. In other words, $\omega'$ is the shape function for the horizontal velocity. In the following, we deduce the vertical velocity profile from the ice mass conservation relationship.

Owing to the incompressibility of the ice, the ice mass conservation can be written as:

$$\frac{1}{Y} \frac{\partial (Y u_x)}{\partial x} + \frac{\partial u_z}{\partial z} = 0. \quad (4)$$
Multiplying Equation (3) by Y and differentiating with respect to x leads to:

$$ \frac{\partial Y u_x}{\partial x} = \frac{\partial Y \bar{U}_x}{\partial x} \omega'(\zeta) + Y \bar{U}_x \omega''(\zeta) \frac{\partial \zeta}{\partial x'} \quad (5) $$

and the definition $\zeta = (z - B)/H$ leads to:

$$ \frac{\partial \zeta}{\partial x} = -\frac{1}{H} \left( \frac{\partial B}{\partial x} + \zeta \frac{\partial H}{\partial x} \right). \quad (6) $$

Moreover, the ice mass conservation applied to a whole ice column leads to:

$$ \frac{\partial Q}{\partial x} = Y(a - m). \quad (7) $$

Now

$$ Q = H Y \bar{U}_x, \quad (8) $$

and therefore

$$ \frac{\partial Q}{\partial x} = \frac{\partial H}{\partial x} Y \bar{U}_x + H \frac{\partial Y \bar{U}_x}{\partial x}. \quad (9) $$

Finally, replacing $\partial\zeta/\partial x$ and $\partial Y \bar{U}_x/\partial x$ in Equation (5) using Equations (6), (7) and (9) leads to:

$$ \frac{\partial Y u_x}{\partial x} = \frac{1}{H} \left( a - m - \frac{\partial H}{\partial x} \bar{U}_x \right) Y \omega'(\zeta) - \frac{1}{H} Y \bar{U}_x \left( \frac{\partial B}{\partial x} + \zeta \frac{\partial H}{\partial x} \right) \omega''(\zeta). \quad (10) $$

From Equation (4), $u_z$ can be written as:

$$ u_z = u_z(z = B) - \int_B^z \frac{1}{Y} \frac{\partial (Y u_x)}{\partial x} dz'. \quad (11) $$

Now $u_z(z = B)$ is given by:

$$ u_z(z = B) = -m + \frac{\partial B}{\partial x} \bar{U}_x s, \quad (12) $$

where $s = \omega'(\zeta = 0)$ is the sliding rate ($s = 0$ if no sliding and $s = 1$ if full sliding, i.e. plug flow). Therefore

$$ u_z = -m + \frac{\partial B}{\partial x} \bar{U}_x s - \left( a - m - \frac{\partial H}{\partial x} \bar{U}_x \right) \omega(\zeta) + \bar{U}_x \frac{\partial B}{\partial x} (\omega'(\zeta) - s) + \bar{U}_x \frac{\partial H}{\partial x} \int_0^\zeta \zeta' \omega''(\zeta') d\zeta'. \quad (13) $$

An integration by parts gives:

$$ \int_0^\zeta \zeta' \omega''(\zeta') d\zeta' = \zeta \omega'(\zeta) - \omega(\zeta), \quad (14) $$

leading finally, with two cancellations, to

$$ u_z = -m - (a - m)\omega(\zeta) + \bar{U}_x \left( \frac{\partial B}{\partial x} + \zeta \frac{\partial H}{\partial x} \right) \omega'(\zeta). \quad (15) $$

The vertical velocity in the $\zeta$ coordinate, $u_\zeta = d\zeta/dt$, can be written in terms of the vertical velocity in the $z$ coordinate, $u_z$, using the chain rule, as:

$$ u_\zeta = \frac{1}{H} \left[ u_z - \left( \frac{\partial B}{\partial x} + \zeta \frac{\partial H}{\partial x} \right) u_x \right], \quad (16) $$

leading to

$$ u_\zeta = -\tilde{m} - (\tilde{a} - \tilde{m})\omega(\zeta), \quad (17) $$

with $\tilde{a} = a/H$ and $\tilde{m} = m/H$ the normalized accumulation and melting rates. Notice that this form removes the dependence of the vertical velocity in $\zeta$ coordinate on the basal and surface slopes.

Defining the ratio $\mu$, relating basal melting to mass balance, as:

$$ \mu(x) = \frac{m(x)}{a(x) - m(x)}, \quad (18) $$

we can write

$$ u_\zeta = - \frac{\omega(\zeta) + \mu(x)}{T(x)}, \quad (19) $$

where $T(x) = H(x)/(a(x) - m(x))$ is called the local characteristic time. We assume here that $\mu$ is uniform all along the flow-tube, which results in a spatially uniform shape function for the vertical velocity, $u_\zeta$. To further the analysis we change the coordinate system from $(x, \zeta)$ to a new system $(\pi, \theta)$ defined by:

$$ \pi(x) = \ln \left( \frac{Q(x)}{Q_{\text{ref}}} \right), \quad (20) $$

$$ \theta(\zeta) = \ln \left( \frac{\omega(\zeta) + \mu}{1 + \mu} \right), \quad (21) $$

where $Q_{\text{ref}}$ is a reference flux at a reference position $x_{\text{ref}}$. In the following, we will call $\pi$ the logarithmic-flux horizontal coordinate and $\theta$ the logarithmic-flux vertical coordinate. Note that the change of variable from $x$- to $\pi$-coordinate requires $Q(x)$ to be an increasing function, which is the case when the accumulation rate, $a$, is strictly positive all along the flow-tube. Similarly, the change of variable from $\zeta$- to $\theta$-coordinate requires $\omega(\zeta)$ to be a strictly increasing function, which is the case when $\omega'$ is a strictly positive function of $\zeta$.

We then define the horizontal and vertical velocity components in this new coordinate system as $u_\pi = d\pi/dt$ and $u_\theta = d\theta/dt$. The logarithmic-flux horizontal velocity $u_\pi$ can be derived from the horizontal velocity in the standard $x$-coordinate,

$$ u_\pi = \frac{1}{Q} \frac{\partial Q}{\partial x} \frac{dx}{dt}. \quad (22) $$

Replacing $Q$ using Equation (8), $\partial Q/\partial x$ using Equation (7) and $dx/dt$ using Equation (3) leads to:

$$ u_\pi = \frac{1}{T(\pi)\psi(\theta)}, \quad (23) $$

where $\psi(\theta)$ is the inverse of the horizontal-velocity shape function $\omega'$ written in the $\theta$-coordinate,

$$ \psi(\theta) = \psi \left[ \ln \left( \frac{\omega(\zeta) + \mu}{1 + \mu} \right) \right] = \frac{1}{\omega'(\zeta)}, \quad (24) $$

and is called the time shape function. We have used a different notation here for $1/\omega'$ to emphasize that it is the derivative of $\omega$ with respect to $\zeta$, and not with respect to $\theta$.

Similarly, we can derive the vertical velocity in the logarithmic-flux vertical coordinate, $\theta$, from the corresponding quantity in the normalized vertical coordinate, $\zeta$,

$$ u_\theta = \frac{1}{\omega + \mu} \frac{\partial \omega}{\partial \zeta} \frac{d\zeta}{dt}, \quad (25) $$

and replacing $d\zeta/dt$ by its expression in Equation (17) leads to the following expression:

$$ u_\theta = - \frac{1}{T(\pi)\psi(\theta)}. \quad (26) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Scheme for the derivation of the formula for isochrone slope.</figcaption></figure>

We immediately remark that the expressions in Equations (23) and (26) sum to zero. This gives simple linear trajectories for ice particles in this logarithmic-flux coordinate system $(\pi, \theta)$

$$ \left. \frac{d\pi}{d\theta} \right|_c = -1, \qquad (27) $$

which is equivalent to:

$$ \theta = \pi_0 - \pi, \qquad (28) $$

where $\pi_0 = \ln(Q_0/Q_{\text{ref}})$ and $Q_0$ is the total flux at the initial position, $x_0$, of the particle at the surface. Mathematically minded readers will guess that we have used the notation $(d\pi/d\theta)|_c$ because trajectories are the characteristics of the (hyperbolic) age equation. Indeed, as shown in the Appendix, we could have reached the same results for the trajectories using the age equation:

$$ \frac{D\chi}{Dt} = 1, \qquad (29) $$

where $\chi$ is the age, and subsequently using standard theory for solution of hyperbolic equations.

The fact that trajectories are linear in the logarithmic-flux coordinate system $(\pi, \theta)$ is useful for solving the age equation. We will see in the next sections that use of this coordinate system allows analytical determination of the slope of the isochrones everywhere in the ice sheet. It also permits computation of the initial position of a particle from its current position, by taking the exponential of Equation (28)

$$ Q_0 = Q \frac{\omega + \mu}{1 + \mu}. \qquad (30) $$

Similarly, the age, $\chi$, of a particle at a given position $(\pi, \theta)$ is readily computed from

$$ \chi(\pi, \theta) = \int_{\pi+\theta}^{\pi} T(\pi') \psi(\pi - \pi' + \theta) d\pi'. \qquad (31) $$

For plug flow, by definition $\omega' \equiv 1$ (i.e. uniform horizontal velocity along a vertical line). This situation corresponds to 100% of the motion arising from basal sliding or,

equivalently, to ice deformation concentrated at the extreme base of the ice sheet. In this case, the velocity field becomes:

$$ u_{\pi} = \frac{1}{T(\pi)}, \qquad (32) $$

$$ u_{\theta} = -\frac{1}{T(\pi)}, \qquad (33) $$

with Equation (27) also applying to plug flow. The formula for the age simplifies to:

$$ \chi(\pi, \theta) = \int_{\pi+\theta}^{\pi} T(\pi') d\pi'. \qquad (34) $$

# 3. SLOPE OF ISOCHRONES

## 3.1. General case

We will determine here an analytical formula for the slope of the isochrones. The isochrone slope is dependent on the coordinate system, and in the future it will be qualified according to the coordinate system, and called, for example, the $(\pi, \theta)$ slope, the $(x, \zeta)$ slope or the $(x, z)$ slope. We will focus in particular on the $(x, \zeta)$ slope, which is an observable, as the $(x, z)$ slope is strongly affected by the local bedrock and surface slope.

Consider two particles initially at the upper surface of the ice sheet at arbitrary neighbouring positions $x_0, x_0 + \delta x_0$ (corresponding to $\pi_0, \pi_0 + \delta \pi_0$). After time $\Delta t$ they will have positions $x, x + \delta x$ (corresponding to $\pi, \pi + \delta \pi$) with elevations $\zeta, \zeta + \delta \zeta$ (corresponding to $\theta, \theta + \delta \theta$). Our goal is to derive an expression for the isochrone $(x, \zeta)$ slope, but we will first derive the $(\pi, \theta)$ slope. As indicated by the notation, we will assume that $\delta \pi$, etc. are small quantities.

As a consequence of Equation (27), the trajectories in the $(\pi, \theta)$ space of both particles are two parallel lines with a slope of -1. The vertical shift of the two lines is equal to $\delta \theta_{\text{traj}} = \delta \pi_0$, and we have $\delta \pi = \delta \theta_{\text{traj}} - \delta \theta$ (see Fig. 1). This leads to the following expression for the slope of the isochrones (assuming that $\delta \pi_0 \to 0$):

$$ \left. \frac{d\theta}{d\pi} \right|_{\chi} = \frac{\delta\theta}{\delta\pi} = \frac{\delta\pi_0 - \delta\pi}{\delta\pi}. \qquad (35) $$

The notation $|_{\chi}$ is used throughout this paper to define constant age, i.e. an isochrone layer. To proceed further, we must relate $\delta \pi$ to $\delta \pi_0$. Firstly, we need to evaluate the time,
$\delta t_{\pi_0}$, required for the first particle to travel from $\pi_0$ to $\pi_0 + \delta\pi_0$ and the time, $\delta t_{\pi}$, for the second particle to travel from $\pi$ to $\pi + \delta\pi$ (see Fig. 1). When $\delta\pi, \delta\pi_0$ are small, these quantities can be derived from Equation (23) as:

$$ \delta t_{\pi_0} = T_0 \psi_0 \delta \pi_0, \quad (36) $$

$$ \delta t_{\pi} = T(\pi)\psi(\theta)\delta\pi, \quad (37) $$

with $T_0 = T(\pi_0)$ and $\psi_0 = \psi(\theta = 0)$.

Secondly, we need to evaluate the difference in duration between the two journeys of the particles from $\pi_0 + \delta\pi_0$ to $\pi$, represented by $\delta t_{\pi_0-\pi}$. In plug flow, $\delta t_{\pi_0-\pi} = 0$ because the horizontal velocity does not depend on the vertical position of the particle. In the general case, we can determine $\delta t_{\pi_0-\pi}$ using Equation (23) to be:

$$ \delta t_{\pi_0 \to \pi} = \int_{\pi_0 + \delta\pi_0}^{\pi} T(\pi') [\psi(\theta) - \psi(\theta - \delta\theta_{\text{traj}})] d\pi', \quad (38) $$

or, if we differentiate the $\psi$ function,

$$ \delta t_{\pi_0 \to \pi} = - \int_{\pi_0 + \delta\pi_0}^{\pi} T(\pi') \frac{\partial\psi}{\partial\theta}(\theta) \delta\theta_{\text{traj}} d\pi'. \quad (39) $$

Using the fact that $\delta\theta_{\text{traj}}$ is constant and equal to $-\delta\pi_0$ as illustrated in Figure 1, we find:

$$ \delta t_{\pi_0 \to \pi} = \delta \pi_0 \int_{\pi_0 + \delta \pi_0}^{\pi} T(\pi') \frac{\partial \psi}{\partial \theta}(\theta) d\pi', \quad (40) $$

or, replacing $\theta$ by $\pi_0 - \pi$ by Equation (28) and taking the limit as $\delta\pi_0 \to 0$,

$$ \delta t_{\pi_0 \to \pi} = \beta \delta \pi_0, \quad (41) $$

where

$$ \beta(\pi, \pi_0) = \int_{\pi_0}^{\pi} \left[ T(\pi') \frac{\partial \psi}{\partial \theta}(\pi_0 - \pi') \right] d\pi'. \quad (42) $$

Note for plug flow that $\beta = 0$ as $\psi \equiv 1$. This is why we call $\beta$ the non-plug-flow parameter.

Now, by construction, both particles have the same age:

$$ \delta t_{\pi} = \delta t_{\pi_0} - \delta t_{\pi_0-\pi}, \quad (43) $$

and replacing $\delta t_{\pi_0}$ and $\delta t_{\pi}$ by Equations (36) and (37), we get:

$$ T\psi\delta\pi = T_0\psi_0\delta\pi_0 - \beta(\pi, \pi_0)\delta\pi_0, \quad (44) $$

or

$$ \delta\pi = (1 - \alpha)\delta\pi_0, \quad (45) $$

with

$$ \alpha = 1 - \frac{1}{T\psi} [T_0\psi_0 - \beta(\pi, \pi_0)]. \quad (46) $$

Finally we replace $\delta\pi$ by $(1 - \alpha)\delta\pi_0$ in Equation (35) to reach the first goal of deriving an expression for the $(\pi, \theta)$ isochrone slope:

$$ \left. \frac{d\theta}{d\pi} \right|_{\chi} = \frac{\alpha}{1 - \alpha}. \quad (47) $$

The $(\pi, \theta)$ slope depends on the value of the characteristic time $T = H/(a - m)$ at the initial and final positions of the particle, $\pi_0$ and $\pi$, and depends also on the value of $T$ between $\pi_0$ and $\pi$ through the function $\beta$. The initial position of the particle can be evaluated from Equation (30) and the slope of the isochrone can be estimated directly by integrating Equation (42). Because the changes of coordinates $x \to \pi$ and $\zeta \to \theta$ are monotonic, the signs of the $(x, \zeta)$ and $(\pi, \theta)$ slopes are the same. Moreover, assuming that $\omega''$ is positive, $\beta$ is negative and thus, from Equation (46), $(1 - \alpha)$

is positive. We will thus call $\alpha$ the slope-sign parameter, because its sign determines the sign of the isochrone slope.

From Equations (20) and (21), we see that $d\pi = dQ/Q$, $d\theta = d\omega/(\omega + \mu)$ and we also use the identity $d\omega = \omega'd\zeta$. With $dQ = (a - m)Y dx$ from Equation (7), we deduce the $(x, \zeta)$ slope,

$$ \left. \frac{d\zeta}{dx} \right|_{\chi} = \frac{\alpha}{1 - \alpha} \frac{a(x) - m(x)}{Q(x)} Y(x) \frac{\omega + \mu}{\omega'}. \quad (48) $$

In this formula, note that $Q/Y$ represents the local flux $\bar{Q}$,

$$ \bar{Q} = \frac{Q}{Y} = \int_{B}^{S} u_{x}(z) dz. \quad (49) $$

Now, since $Q = \int_{0}^{x} (a - m)Ydx$ we can write the isochrone slope as:

$$ \left. \frac{d\zeta}{dx} \right|_{\chi} = \frac{\alpha}{1 - \alpha} \lambda(x) \frac{\omega(\zeta) + \mu}{\omega'(\zeta)}, \quad (50) $$

with $\lambda(x)$ defined by:

$$ \lambda(x) = \left[ \int_{0}^{x} \frac{(a(x') - m(x'))Y(x')}{(a(x) - m(x))Y(x)} dx' \right]^{-1}. \quad (51) $$

Note that $\lambda(x)$ is simply $1/x$ when accumulation rate, $a$, melting rate, $m$, and flow-tube width, $Y$, are all independent of $x$.

A second, related, expression for $\alpha$, which illustrates different aspects of the problem, can be derived as follows. Integrating Equation (42) by parts, $\beta$ can be expressed as:

$$ \beta = T_0 \psi_0 - T \psi + \int_{\pi_0}^{\pi} \frac{\partial T}{\partial \pi} \psi(\pi_0 - \pi') d\pi', \quad (52) $$

leading, with Equation (46), to

$$ \alpha = \frac{1}{T\psi} \int_{\pi_0}^{\pi} \frac{\partial T}{\partial \pi} \psi(\pi_0 - \pi') d\pi'. \quad (53) $$

The interest of this formula is that the new integral term in $\alpha$ contains the derivative of $T$ with respect to $\pi$, instead of the derivative of $\psi$ with respect to $\theta$. It can now be seen immediately that the $(x, \zeta)$ slope of the isochrones is zero when the particle has experienced constant $T$ conditions along its journey, which is not obvious from the previous formulation. We will use the term $\zeta$-flat to qualify these isochrones, because it is equivalent to having a zero slope in the $(x, \zeta)$-, $(Q, \zeta)$- or $(\pi, \zeta)$-coordinate systems. Similarly, we will use $z$-flat for isochrones with a zero slope in the $(x, z)$-, $(Q, z)$- or $(\pi, z)$-coordinate systems.

## 3.2. Plug flow as a special case

In plug flow, $\omega = \zeta, \omega' \equiv 1$ and $\psi \equiv 1$, lead to $\beta = 0$. Thus

$$ \alpha = 1 - \frac{T_0}{T}, \quad (54) $$

and the $(\pi, \theta)$ and $(x, \zeta)$ slopes, respectively, become

$$ \left. \frac{d\theta}{d\pi} \right|_{\chi} = \frac{\alpha}{1 - \alpha}, \quad (55) $$

and

$$ \left. \frac{d\zeta}{dx} \right|_{\chi} = \frac{\alpha}{1 - \alpha} \frac{a(x) - m(x)}{Q(x)} Y(x)(\zeta + \mu), \quad (56) $$

or

$$ \left. \frac{d\zeta}{dx} \right|_{\chi} = \frac{\alpha}{1 - \alpha} \lambda(x)(\zeta + \mu). \quad (57) $$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Isochrones over a step change in bedrock elevation at $x = x_1$ (corresponding to $\pi = \pi_1$), marked by a thick grey vertical line. Ice is in plug flow. (a) $(x, z)$-coordinate system. (b) $(\pi, \theta)$-coordinate system. In both panels, thin black lines (plain and dot-dashed) are isochronous layers, and the thick dashed line is the trajectory of a particle deposited immediately upstream of the bedrock step (see section 4.1 for a complete description of the zones).</figcaption></figure>

The remarkable aspect of Equation (56) is that the $(x, \zeta)$ slope depends only upon the elevation, the initial (through $\alpha$) and final positions of the particle, and the local flux $\bar{Q} = Q/Y = \int_B^S u_x(z) dz$. The intervening bedrock geometry does not affect the isochrone slope.

# 4. ILLUSTRATIONS

We illustrate the effect of various parameters on the slope of the isochrones with three different cases, which are combinations of steps and plateaux associated with plug flows and flows with distributed internal deformation. These examples are presented because they are simple. For example, the $T(\pi)$ parameter takes only two different values, because the bedrock elevation takes only two different values, separated by steps. As a consequence, the solutions are entirely analytical.

It is clear that at the steps ice particles are moving upwards at infinite speed. In reality, horizontal deviator stress gradients would prevent this, but act over length scales comparable with the ice thickness. Without these stress gradients, the predicted $(x, z)$ slope is infinite, but the $(x, \zeta)$ slope remains well behaved, as it depends upon the integral of the parameters it experiences through its journey. Thus, after any step, the isochrones emerge with unchanged $\zeta$. Indeed, we can see that Equation (17) still holds when, at some points, $\partial B/\partial x$ and $\partial H/\partial x$ tend to infinity and the $\tilde{a}$ and $\tilde{m}$ functions are discontinuous. The derivation for the slope of the isochrones is consequently not affected. This example can be seen as the limit of cases when the bedrock slope tends to infinity.

It should also be noted that, at short horizontal length scales, the assumptions regarding the spatially constant velocity shape function do not hold in cases with abrupt bedrock reliefs. These illustrations are only intended to show the use of our analytical solution to understand the geometry of isochrone layers on long horizontal length scales.

However, we do not want to give the impression that our approach only holds for steps in bedrock topography or accumulation rate. The results for plug flow – that the slope depends only on initial and final conditions – holds for any intermediate variation. For a given flux shape function, $\omega$, the isochrone slopes can be computed as follows. First, for

every x column, we can compute the corresponding total flux $Q$ (by an integration of the mass balance from the ice divide), the logarithmic-flux horizontal variable $\pi$ (directly deduced from $Q$) and the x-dependent scale factor (also by an integration from the ice divide). Second, for every $\zeta$ level, we can compute the logarithmic-flux vertical variable, $\theta$, the $\omega'$ value and the depth-dependent scale factor. Third, for every $(x, \zeta)$ node, we can compute the spatial origin of the particle from Equation (28) or (30), the $\beta$ term from the convolution-type Equation (42) (by a simple numeric integration against the horizontal $\pi$ variable) and the principal term $\alpha/(\alpha - 1)$ (deduced directly from $\beta$). We expect this numerical method to be significantly faster and more accurate than a brute-force method of back-tracking the trajectories of the particles by a finite-difference approach.

## 4.1. Step in $H/a$, plug flow

In this first experiment, we assume (1) plug flow and no lateral divergence ($Y = \text{constant}$), (2) accumulation rate to be constant and equal to $3 \text{ cm a}^{-1}$ of ice, (3) no basal melting all along the flow-tube ($\mu = 0$) and (4) ice thickness equal to $H_1 = 4000 \text{ m}$ from the divide to $30 \text{ km}$ downstream (position $x = x_1$, corresponding to $\pi = \pi_1$; see Fig. 2), then decreasing abruptly to $H_2 = 2000 \text{ m}$, and keeping this value downstream. The ice surface is flat.

Results are plotted in Figure 2, where we can distinguish three different areas, depending on the initial and final horizontal positions of the ice particles:

*   **Area A1**, before the bedrock step. The ice particles in this zone have trajectories entirely in this area, and have thus experienced constant accumulation, bedrock and surface elevation conditions during their motion. Hence, the isochrones are flat from Equations (54) and (56) (or even simply by the integration of Equation (17), because the rate of elevation change is independent of the horizontal coordinate).

*   **Area A2** represents particles whose trajectory starts upstream of the bedrock step but which have experienced the bedrock step. In this case, applying Equation (55) with $\alpha = 1 - (H_1/H_2)$ gives the $(\pi, \theta)$ slope of the isochrones as a constant of $-1/2$. Equation (57) gives
<figure><img src="image_4.png" /><figcaption>Fig. 3. Isochrones over a bedrock plateau delimited at $x = x_1$ and $x = x_2$ (corresponding to $\pi = \pi_1$ and $\pi = \pi_2$). See section 4.2 for a complete description of the experiment. Ice is in plug flow. Details as for Figure 2.</figcaption></figure>

the $(x, \zeta)$ slope of the isochrones as:

$$ \left. \frac{d\zeta}{dx} \right|_{\chi} = -\frac{\zeta}{2x}. \qquad (58) $$

The position of the step does not appear as might be expected because $x$ is proportional to the total flux $Q$.

Area A3 represents particles whose trajectory starts downstream of the bedrock step. As with area A1, these particles have experienced constant accumulation and ice-thickness conditions, and thus the slope of the isochrones is zero.

The dashed line separating areas A2 and A3 is a trajectory of ice particles (or characteristic, if we follow the formalism of the hyperbolic equations; see Appendix). It is clear from the illustrations that when considering a step in bedrock elevation the age solution is non-differentiable across these trajectories or characteristics. Such non-differentiability does not indicate the invalidity of the solutions, since such a property is expressly known to be a possible feature of solution fields across characteristics of hyperbolic equations (Rubinstein and Rubinstein, 1998). This is also true for the following examples.

Finally, the solution shows where the 1-D approximation is invalid. This approximation, where horizontal advection is neglected, is often used in the study of isochrones. The non-zero $(x, \zeta)$ slopes in this example are a consequence of horizontal advection; the equivalent 1-D solution would have computed the isochrones in area A3 to have filled the whole zone downstream of the bedrock step, and the isochrones would have been flat; in other words, area A2 would not have existed.

## 4.2. Plateau in $H/a$, plug flow

In this example, the assumptions are the same as in section 4.1, except that the bedrock elevation decreases abruptly from 2000 m to 0 m at 60 km (position $x = x_2$, corresponding to $\pi = \pi_2$; see Fig. 3), creating a plateau in the bed.

We can distinguish several different areas:

Areas A1–A3 are identical to the corresponding cases in the example of section 4.1.

Area A4 represents particles whose trajectory starts before the first bedrock step and finishes after the second

step. From Equation (54), $\alpha = 0$ and the isochrones are consequently $\zeta$-flat and $z$-flat. This is not a priori trivial to understand, and we see here a concrete application of the formula of isochrone slope derived in section 3. In the next example we will see that isochrones are not $\zeta$-flat in area A4 when there is internal deformation.

Area A5 represents particles whose trajectory starts between the two bedrock steps but continues to an area downstream of the second step. If we apply Equation (47) with $\alpha = 1 - (H_2/H_1)$, we find that the slope of the isochrones in the logarithmic-flux coordinate system $(\pi, \theta)$ is constant and equal to 1, and Equation (57) gives the slope of the isochrones in the $(x, \zeta)$-coordinate system as:

$$ \left. \frac{d\zeta}{dx} \right|_{\chi} = \frac{\zeta}{x}. \qquad (59) $$

Area A6 represents particles whose trajectory begins downstream of the second bedrock step. As with area A1, these particles experience constant accumulation and ice-thickness conditions, and thus the slope of the isochrones is zero and the depth of the isochrones is the same as in A1.

Similarly to the previous example, areas A2, A4 and A5 would not be computed correctly by the 1-D approximation.

## 4.3. Plateau in $H/a$, distributed deformation

In this example we use the same profiles for accumulation, bedrock and surface elevations and lateral divergence as in the previous example, but a different flux shape function, $\omega$, which permits internal deformation,

$$ \omega = \zeta^p, \qquad (60) $$

$$ \omega' = p\zeta^{p-1}, \qquad (61) $$

with $p = 1.5$ chosen in our example. We call $p$ the deformation exponent. For $1 \le p \le 2$, these simulate typical velocity profiles for simple shear. This shape function is qualitatively similar to standard glaciological forms:

$$ \omega = \frac{1}{n+1} \left[ (1 - \zeta)^{n+2} + (n+2)\zeta - 1 \right], \qquad (62) $$

leading to

$$ \omega' = \frac{n+2}{n+1} \left[ 1 - (1 - \zeta)^{n+1} \right], \qquad (63) $$
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Isochrones over a bedrock plateau delimited at $x = x_1$ and $x = x_2$ (corresponding to $\pi = \pi_1$ and $\pi = \pi_2$). See section 4.3 for a complete description of the experiment. Ice is experiencing internal deformation. Details as for Figure 2.</figcaption></figure>

but is analytically more tractable. A difference that requires comment is the fact that $\omega''$ (which is a multiple of the shear rate) presents a singularity at $\zeta = 0$ (it tends to $+\infty$ at $\zeta = 0$), while $\omega''$ was used during the analytical development in section 2. This singularity is not important, because only the non-singular integrated forms $\int \omega''$ or $\int \zeta \omega''$ appear in the analytical development.

Using our additional assumption that $\mu = 0$, it is straightforward to show that $\psi$ is given by:

$$ \psi(\theta) = \frac{1}{p} \exp\left(\frac{1-p}{p}\theta\right). \qquad (64) $$

The slope of the isochrones in the logarithmic-flux coordinate system ($\pi, \theta$) is given by Equation (47), with $\alpha$ given by Equation (53),

$$ \alpha = \frac{1}{\psi T} \int_{\pi_0}^{\pi} \frac{\partial T}{\partial \pi} \psi(\pi_0 - \pi') d\pi'. \qquad (65) $$

As in section 4.2, we can distinguish the same six areas, A1–A6, with respect to the initial and final horizontal positions of the ice particles (see Fig. 4). Isochrones are still flat in areas A1, A3 and A6 because these ice particles have experienced constant accumulation and ice-thickness conditions during their trajectories. For area A2, $\partial T / \partial \pi = 0$ between initial and final horizontal positions ($\pi_0$ and $\pi$), except for the bedrock step $\Delta B$ at $\pi = \pi_1$ which leads to a step $\Delta_1 T = T_2 - T_1$ in the $T$ parameter. Thus, replacing $\pi_0$ using Equation (28), we obtain:

$$ \alpha = \frac{\Delta_1 T \psi[\theta + (\pi - \pi_1)]}{T_2 \psi(\theta)}, \qquad (66) $$

and

$$ \alpha = -\exp\left[\frac{1-p}{p}(\pi - \pi_1)\right], \qquad (67) $$

giving the following expression for the slope of the isochrones:

$$ \left. \frac{d\theta}{d\pi} \right|_{\chi} = - \frac{\exp\left[\frac{1-p}{p}(\pi - \pi_1)\right]}{1 + \exp\left[\frac{1-p}{p}(\pi - \pi_1)\right]}. \qquad (68) $$

Immediately downstream of the bedrock step at $\pi = \pi_1$, the ($\pi, \theta$) slope is $-1/2$ and decreases with $\pi$ (with a minimum value of $-1$).

Similarly, one can show that in area A5 the ($\pi, \theta$) slope of the isochrones is given by:

$$ \left. \frac{d\theta}{d\pi} \right|_{\chi} = \frac{\exp\left[\frac{1-p}{p}(\pi - \pi_2)\right]}{2 - \exp\left[\frac{1-p}{p}(\pi - \pi_2)\right]}, \qquad (69) $$

where $\pi_2$ is the horizontal coordinate of the second bedrock step (see Fig. 4). Again, the ($\pi, \theta$) slope is 1 just after the bedrock step $\pi_2$ and decreases toward 1/2 as $\pi$ increases.

For area A4, one can also show, from Equation (65), that:

$$ \alpha = \frac{\Delta_2 T \psi[\theta + (\pi - \pi_2)] - \psi[\theta + (\pi - \pi_1)]}{T_1 \psi(\theta)}, \qquad (70) $$

with $\Delta_2 T = T_3 - T_2$. This leads to

$$ \alpha = \frac{1}{2} \exp\left[\frac{1-p}{p}(\pi - \pi_1)\right] - \frac{1}{2} \exp\left[\frac{1-p}{p}(\pi - \pi_2)\right], \qquad (71) $$

giving the following formula for the ($\pi, \theta$) slope of the isochrones:

$$ \left. \frac{d\theta}{d\pi} \right|_{\chi} = \frac{\exp\left[\frac{1-p}{p}(\pi - \pi_2)\right] - \exp\left[\frac{1-p}{p}(\pi - \pi_1)\right]}{1 - \exp\left[\frac{1-p}{p}(\pi - \pi_2)\right] - \exp\left[\frac{1-p}{p}(\pi - \pi_1)\right]}. \qquad (72) $$

Unlike the previous example, isochrones are not $\zeta$-flat in this area (see Fig. 4), as a result of the distributed deformation in the flow. It suggests that, in this case, we might be able to reconstruct the flux shape function, $\omega$, from the isochrone geometry.

# 5. DISCUSSION

## 5.1. Influence of accumulation, ice thickness and lateral flow-tube divergence

The analytical solution for the slope of the isochrones where the velocity shape functions are spatially uniform allows us to study the influence of the spatial variation of accumulation, ice thickness and lateral flow-tube divergence on the slope of the isochrones. Two primary controls on isochrones are the bedrock and surface elevation variations: an isochronic layer has an approximately constant relative depth in the ice sheet, so that plotting the isochrones with the normalized coordinate, $\zeta$, eliminates this primary effect. In other words, the $(x, \zeta)$ slope informs us about the processes governing ice flow.
From Equation (50), the $(x, \zeta)$ slope is the product of three terms: one principal term, and two scale factors. The scale factors can modify the amplitude of the $(x, \zeta)$ slope but not its sign.

One scale factor, $\lambda(x)$, depends upon the horizontal position of the particle, and the other, $(\omega(\zeta) + \mu)/\omega'(\zeta)$, depends upon the normalized elevation, $\zeta$, of the particle. The x-dependent scale factor tends to zero as the distance from the ice divide increases. It also depends on the flow-tube mass-balance term (accumulation rate minus melting rate multiplied by the flow-tube width) at the local position compared to its values upstream. If the accumulation increases downstream, or if the flow-tube width increases downstream, the scale factor decreases more slowly with distance, $x$, from the ice divide.

The depth-dependent scale factor depends on the flux shape function, $\omega$. For plug flow, it simply increases linearly from $\mu$ at the bedrock to $1 + \mu$ at the surface. When there is no melting, the depth-dependent scale factor tends to zero as $\zeta$ tends to zero. That does not mean, however, that the isochrones approach a horizontal orientation as they approach the bedrock interface, because the principal term, $\alpha/(1 - \alpha)$, may tend to $+\infty$. When there is melting but no sliding, the depth-dependent scale factor tends to $+\infty$. That does not mean, however, that the isochrones approach a vertical orientation as they approach the bedrock interface, because the slope-sign term, $\alpha$, tends to zero.

The principal term $\alpha/(1 - \alpha)$ depends on the characteristic time $T = H/(a - m)$ (ice thickness divided by accumulation rate minus melting rate) between the initial and final horizontal positions of the particle. Changes in mass balance (accumulation minus melting) or in ice thickness have exactly the same influence on this term. Lateral flow-tube divergence does not impact this principal term and thus the sign of the $(x, \zeta)$ slope, but only the x-dependent scale factor. For plug flow, the situation is even simpler, because this term depends only upon the ratio between the initial and the final values of this local characteristic time (and not the values at intermediate locations along the particle trajectory); if it is larger at the final position, the $(x, \zeta)$ isochrone slope is positive and, conversely, if it is smaller the $(x, \zeta)$ slope is negative. For general flow, the basic pattern remains the same, but the values of the characteristic time in between enter in the $\alpha$ parameter in Equation (46) through $\beta$, the latter being defined in Equation (42). It should be emphasized here that through our analytical solution we can now assess the influence of the bedrock geometry and spatial accumulation rate variations on isochrone layers. Such a formulation would have been rather more difficult to detect with numerical experiments with a finite number of bedrock geometry and accumulation rate scenarios.

## 5.2. Limitations of assumptions

In the derivation of our analytical model, we made several assumptions.

Firstly, we assumed that the ice sheet is in steady state. This assumption implies that our results should apply with least adjustment to the central parts of the ice sheets, where thickness changes are subdued compared with the margins. For Greenland or for West Antarctica, the stable Holocene period occupies a large part of the ice-sheet thickness, and thus this assumption should not be too constraining

providing we do not consider isochrones too close to the bedrock. However, for central East Antarctica, the Holocene period occupies only a small fraction of the total ice thickness. Counteracting this is the fact the ice thickness conditions have been relatively stable, with changes of the order of 200 m (Ritz and others, 2001), i.e. less than 10%. There are hints from numerical calculations (Nereson and others, 2000) that even severe changes in geometry do not affect the isochrone geometry enormously, at least where the flow is horizontally uniform.

The fact that the accumulation rate changed between interglacial and glacial periods is not so important, because the isochrone slope depends only on the relative spatial variations of the accumulation. Indeed, if we assume that spatial and temporal variations of accumulation and melting can be separated as:

$$a(x, t) = \bar{a}(x)R(t), \qquad (73)$$

$$m(x, t) = \bar{m}(x)R(t), \qquad (74)$$

(with $R(t)$ a non-negative multiplication factor) we can come back to the steady-state problem with melting rate $\bar{m}(x)$ and accumulation rate $\bar{a}(x)$ by changing the time variable from $t$ to $\bar{t}$ defined by:

$$\bar{t} = \int_{0}^{t} R(t) \, dt. \qquad (75)$$

Because of the bijection between $t$ and $\bar{t}$ (for every value of $t$ there is one unique value of $\bar{t}$) isochrones in $t$ are equivalent to isochrones in $\bar{t}$. The true age, $\chi$, can simply be deduced from the normalized age, $\bar{\chi}$, with Equation (73). This technique to include transient accumulation rates has been used by previous workers (Schott and others, 1992; Nereson and others, 1996). With this hypothesis, the particles follow the same paths regardless of the time-dependent term $R(t)$ and they go faster or slower as $R(t)$ increases or decreases, respectively. There is no surface adjustment and no delay in velocity response to mass-balance changes.

Secondly, we assumed a uniform ratio, $\mu = m/(a - m)$, relating basal melting to mass balance. This assumption applies in particular if there is no basal melting. If there is significant basal melting, it is very unlikely that $\mu$ is uniform and thus our analytical model is not valid. However, at shallow depths, the velocity field is almost unaffected by basal melting, and thus the error implied by our assumption should be small.

Thirdly, we assumed a uniform flux function, $\omega$. Although we have not explicitly made any mechanical assumptions, relying purely on kinematics, we are clearly motivated by the shallow-ice approximation (SIA), because the assumption of the velocity shape function is consistent with the SIA (Lliboutry, 1979; Hutter, 1983). More discussions on this can be found in Hindmarsh (2001). Near bedrock steps, the SIA will not hold over short distances, much shorter than those of an appreciable fraction of the ice-sheet extent that we are concerned with here. In these narrow zones, the shape function will not be uniform. Over longer distances, the assumption of a horizontally uniform flux shape function, $\omega$, will hold as long as (1) the basal conditions are constant (e.g. no subglacial lakes which permit increased sliding) and (2) the temperature profile does not change too much (e.g. no change in geothermal heat flux, only small changes in ice thickness).
## 5.3. Inference of velocity shape functions from observed layers

Bedrock steps lead to traceable, propagating discontinuities in isochrone slopes which correspond to ice particle trajectories and which must have $(\pi, \theta)$ slope of -1. The relationship between the logarithmic-flux vertical coordinate, $\theta$, and the normalized vertical coordinate, $\zeta$, depends on the velocity shape function. The shape function does enter into the formal definition of the coordinate transform so, if all the other relevant parameters are known, we can compute or invert for the shape function.

To change the horizontal coordinate from $x$ to $\pi$, we need to know the mass-balance term $a - m$ (or just the accumulation rate $a$ if basal melting is negligible). This could be done by reconstructing the accumulation rate from the shallow isochrones, which have experienced negligible strain. An alternative way would be to measure the surface velocity; the gradient in surface strain rate gives the accumulation rate. To do that we need, in principle, to know the ratio between the surface velocity and the average velocity $\omega'(\zeta = 1)$. However, if we suppose that the flux shape function, $\omega$, is uniform, gradients in surface velocity still give spatial variations in accumulation, which is sufficient to compute the isochrone slope and the $x \to \pi$ change of coordinate.

To transform the vertical coordinate from $\zeta$ to $\theta$, we need to know the flux shape function, $\omega$, and the ratio (relating basal melting to mass balance), $\mu$. Assuming that we restrict our consideration to an area where there is no melting, and taking the flux shape function, $\omega$, from section 4.3,

$$ \omega = \zeta^p, \qquad (76) $$

we obtain

$$ \theta = p \ln(\zeta). \qquad (77) $$

Because the flux shape function $\omega$ is unknown, the measured isochrone layers cannot be plotted in the $(\pi, \theta)$-coordinate system. As an alternative, we can use a trial coordinate $\tilde{\theta} = \ln(\zeta)$ as the vertical coordinate. In $(\pi, \tilde{\theta})$ space, particle trajectories now have a slope of $1/p$, giving an estimate of the deformation parameter $p$.

More generally, i.e. with an arbitrary flux shape function, $\omega$, we can write:

$$ d\theta = \frac{d\omega}{\omega + \mu}, \qquad (78) $$

$$ d\tilde{\theta} = \frac{d\zeta}{\zeta}, \qquad (79) $$

and thus, if we define $S(\zeta)$ as the slope of the trajectories in $(\pi, \tilde{\theta})$ space, writing $d\omega/d\zeta = (d\omega/d\theta)(d\theta/d\pi)(d\pi/d\tilde{\theta})(d\tilde{\theta}/d\zeta)$ leads to the following differential equation for the flux shape function, $\omega$:

$$ \frac{d\omega}{d\zeta} = -\frac{1}{S(\zeta)} \frac{\omega + \mu}{\zeta}. \qquad (80) $$

In this case, it is therefore, in principle, possible to deduce the flux shape function, $\omega$, from a discontinuity in the slope of the isochrones in the $(\pi, \tilde{\theta})$ space characteristic of a particle trajectory.

## 5.4. Test for numerical age models

With any numerical model, it is essential to test how the discretization of the physical equations affects the results. Our analytical expression may be used to test numerical

methods for the age equation in various ways. In particular, the trajectories of ice particles in the logarithmic coordinate system $(\pi, \theta)$ should be as close as possible to lines of slope equal to -1. These logarithmic coordinates give more weight to the particles close to the bedrock, which tests the accuracy of the solutions at depths where there may be large gradients in age between, for example, finite-difference points. Also, the extent to which abrupt horizontal changes of bedrock or accumulation rate affect the numerical solutions can be readily investigated.

# 6. CONCLUSIONS AND PERSPECTIVES

Our analytical model allows us to understand the first-order variations of isochrone slope: where the velocity shape function is uniform we can characterize how downstream changes in accumulation, ice thickness and lateral flow-tube divergence impact the isochrone slope in the $(x, \zeta)$-coordinate system ($x$ being a horizontal, along-flow coordinate and $\zeta$ being the normalized vertical coordinate). The $(x, \zeta)$ slope of the isochrone at a given point depends on the whole mechanical history of the particle, and thus not only on local but also on upstream accumulation and flow conditions.

By using a new 'logarithmic-flux' coordinate system, we show that the $(x, \zeta)$ slope of the isochrones depends on three terms: one principal term, and two non-negative scale factors which modify the amplitude of the slope but not its sign. The principal term depends only on the characteristic time (ice thickness divided by accumulation rate minus melting rate) between the initial and final positions of the ice particle. For plug flow, only the initial and final values of the characteristic time have an influence. For example, immediately after an upward step in the bedrock, isochrones dip downstream. Increases in accumulation rate have the same effect. Downward steps and decreases in the accumulation rate have the opposite effect on isochrone $(x, \zeta)$ slope. The effects of these changes can be detected over distances that are an appreciable fraction of the ice-sheet extent.

The downstream regions of influence of changes in the bed topography or accumulation rate are clearly bounded. The boundaries, indicated by discontinuities of isochrone $(x, \zeta)$ slope, ultimately relate to the bed topography or accumulation rate regime experienced by a trajectory. We sketched how these slope discontinuities may be used to estimate the vertical profile of velocity, a parameter of primary interest when studying ice flow, but difficult to measure.

We used a steady-state assumption to derive the expression for isochrone $(x, \zeta)$ slope, but in fact it can be relaxed for the accumulation rate as long as the spatial pattern of accumulation remains constant. This is of primary importance when looking at isochrone profiles for Greenland and Antarctica, because these ice sheets have experienced large changes in accumulation rate between glacial and interglacial periods.

We have used an assumption of spatially uniform velocity shape function, but we expect that similar analyses can be applied to the case when the velocity shape functions vary with position, $x$. For example, how do horizontal variations in basal sliding affect isochrone slope (Weertman, 1976)? Or changes in the ratio relating basal melting to surface accumulation? An interesting further step would be to extend the analysis to examples that are not in steady state.
We hope that such analytical expression will help glaciologists qualitatively understand the geometry of isochrones in ice sheets. It may also help define the more appropriate data to invert in a quantitative inverse method based on a physically more sophisticated ice-flow model. In such a refined inverse method, it may also be used as a first guess.

## ACKNOWLEDGEMENTS

We thank C. Ritz for helpful discussions. We warmly thank E. Waddington for his very thorough and constructive review, as well as another anonymous reviewer. Some of this work was funded by UK Natural Environment Research Council grant NER/A/S/2001/01011 'Calculation of ice accumulation rates in Antarctica over the last glacial cycle'.

## REFERENCES

Bailey, J.T., S. Evans and G.de Q. Robin. 1964. Radio echo sounding of polar ice sheets. *Nature*, **204**(4957), 420–421.

Baldwin, D.J., J.L. Bamber, A.J. Payne and R.L. Layberry. 2003. Using internal layers from the Greenland ice sheet, identified from radio-echo sounding data, with numerical models. *Ann. Glaciol.*, **37**, 325–330.

Bell, R.E., M. Studinger, A.A. Tikku, G.K.C. Clarke, M.M. Gutner and C. Meertens. 2002. Origin and fate of Lake Vostok water frozen to the base of the East Antarctic ice sheet. *Nature*, **416**(6878), 307–310.

Conway, H., B.L. Hall, G.H. Denton, A.M. Gades and E.D. Waddington. 1999. Past and future grounding-line retreat of the West Antarctic ice sheet. *Science*, **286**(5438), 280–283.

Fahnestock, M.A., W. Abdalati, S. Luo and S. Gogineni. 2001. Internal layer tracing and age–depth–accumulation relationships for the northern Greenland ice sheet. *J. Geophys. Res.*, **106**(D24), 33,789–33,797.

Hindmarsh, R.C.A. 1999. On the numerical computation of temperature in an ice sheet. *J. Glaciol.*, **45**(151), 568–574.

Hindmarsh, R.C.A. 2001. Notes on basic glaciological computational methods and algorithms. In Straughan, B., R. Greve, H. Ehrentraut and Y. Wang, eds. *Continuum mechanics and applications in geophysics and the environment*. Berlin, etc., Springer-Verlag, 222–249.

Hutter, K. 1983. *Theoretical glaciology; material science of ice and the mechanics of glaciers and ice sheets*. Dordrecht, etc., Reidel Publishing Co./Tokyo, Terra Scientific Publishing.

Jacobel, R.W. and S.M. Hodge. 1995. Radar internal layers from the Greenland summit. *Geophys. Res. Lett.*, **22**(5), 587–590.

Jacobel, R.W., T.A. Scambos, C.F. Raymond and A.M. Gades. 1996. Changes in the configuration of ice stream flow from the West Antarctic ice sheet. *J. Geophys. Res.*, **101**(B3), 5499–5504.

Lliboutry, L. 1979. A critical review of analytical approximate solutions for steady state velocities and temperatures in cold ice sheets. *Z. Gletscherkd. Glazialgeol.*, **15**(2), 135–148.

Nereson, N.A., E.D. Waddington, C.F. Raymond and H.P. Jacobson. 1996. Predicted age–depth scales for Siple Dome and inland WAIS ice cores in West Antarctica. *Geophys. Res. Lett.*, **23**(22), 3163–3166.

Nereson, N.A., C.F. Raymond, R.W. Jacobel and E.D. Waddington. 2000. The accumulation pattern across Siple Dome, West Antarctica, inferred from radar-detected internal layers. *J. Glaciol.*, **46**(152), 75–87.

Ng, F. and H. Conway. 2004. Fast-flow signature in the stagnated Kamb Ice Stream, West Antarctica. *Geology*, **32**(6), 481–484.

Pinglot, J.F., J.O. Hagen, K. Melvold, T. Eiken and C. Vincent. 2001. A mean net accumulation pattern derived from radioactive layers and radar soundings on Austfonna, Nordaustlandet, Svalbard. *J. Glaciol.*, **47**(159), 555–566.

Raymond, C.F. 1983. Deformation in the vicinity of ice divides. *J. Glaciol.*, **29**(103), 357–373.

Reeh, N. 1988. A flow-line model for calculating the surface profile and the velocity, strain-rate, and stress fields in an ice sheet. *J. Glaciol.*, **34**(116), 46–54.

Ritz, C. 1989. Interpretation of the temperature profile measured at Vostok, East Antarctica. *Ann. Glaciol.*, **12**, 138–144.

Ritz, C., V. Rommelaere and C. Dumas. 2001. Modeling the evolution of Antarctic ice sheet over the last 420,000 years: implications for the altitude changes in the Vostok region. *J. Geophys. Res.*, **106**(D23), 31,943–31,964.

Robin, G.de Q., S. Evans and J.T. Bailey. 1969. Interpretation of radio echo sounding in polar ice sheets. *Philos. T. Roy. Soc. London, Ser. A*, **265**(1166), 437–505.

Robin, G.de Q., D.J. Drewry and D.T. Meldrum. 1977. International studies of ice sheet and bedrock. *Philos. T. Roy. Soc. London, Ser. B*, **279**(963), 185–196.

Rubinstein, I. and L. Rubinstein. 1998. *Partial differential equations in classical mathematical physics*, Cambridge, Cambridge University Press.

Schott, C., E.D. Waddington and C.F. Raymond. 1992. Predicted time-scales for GISP2 and GRIP boreholes at Summit, Greenland. *J. Glaciol.*, **38**(128), 162–168.

Siegert, M.J., R. Hodgkins and J.A. Dowdeswell. 1998. A chronology for the Dome C deep ice-core site through radio-echo layer correlation with the Vostok ice core, Antarctica. *Geophys. Res. Lett.*, **25**(7), 1019–1022.

Siegert, M.J., R.C.A. Hindmarsh and G.S. Hamilton. 2003. Evidence for a large surface ablation zone in central East Antarctica during the last Ice Age. *Quat. Res.*, **59**, 114–121.

Siegert, M.J. and 7 others. 2006. Subglacial Lake Ellsworth: A candidate for in-situ exploration in West Antarctica. *Geophys. Res. Lett.*, **31**(23), L23403. (10.1029/2004GL021477.)

Vaughan, D.G., H.F.J. Corr, C.S.M. Doake and E.D. Waddington. 1999. Distortion of isochronous layers in ice revealed by ground-penetrating radar. *Nature*, **398**(6725), 323–326.

Weertman, J. 1976. Sliding–no-sliding zone effect and age determination of ice cores. *Quat. Res.*, **6**(2), 203–207.

Whillans, I.M. 1976. Radio-echo layers and the recent stability of the West Antarctic ice sheet. *Nature*, **264**(5582), 152–155.

# APPENDIX

## THE HYPERBOLIC AGE EQUATION

The characteristic Equations (23), (26) and (27) could have been obtained using the formalism of the hyperbolic equations instead of calculating the horizontal and vertical velocities. Indeed, Hindmarsh (1999) shows that the Lagrangian derivation operator $D(\cdot)/Dt$ in $\zeta$ coordinates can be written as:

$$ \frac{D(\cdot)}{Dt} = \frac{\partial(\cdot)}{\partial t} + u\nabla_H(\cdot) + \frac{\zeta\nabla_H Q - \nabla_H q - a\zeta - m(1 - \zeta)\frac{\partial(\cdot)}{\partial t}}{H} \quad (A1) $$

where $\nabla_H$ is the horizontal gradient in $\zeta$ coordinate and $(\cdot)$ represents the field being advected. In steady-state conditions and with a uniform flux shape function, $\omega$,

$$ \nabla_H Q = a - m, \quad (A2) $$

$$ \nabla_H q = \omega \nabla_H Q = \omega(a - m). \quad (A3) $$

The age equation $D\chi/Dt = 1$ may therefore be written as:

$$ u \nabla_H \chi - \frac{\omega(a - m) + m \frac{\partial \chi}{\partial \zeta}}{H} = 1. \quad (A4) $$

The solution to the age equation gives the age of the ice, and its isolines give isochrones. In particular, note that if the age
solutions are non-differentiable, there will be discontinuities in the slope of the isochrones.

In plane flow

$$u_x \frac{\partial \chi}{\partial x} - \frac{\omega(a-m) + m \partial \chi}{H} \frac{\partial \chi}{\partial \zeta} = 1, \quad (A5)$$

and

$$Q \frac{\partial \chi}{\partial Q} - (\omega + \mu) \frac{\partial \chi}{\partial \omega} = \frac{H}{\nu(a-m)}, \quad (A6)$$

where

$$\nu(\omega(\zeta)) = \omega'(\zeta). \quad (A7)$$

This hyperbolic equation has characteristic equations

$$\frac{d\omega}{\omega + \mu} = -\frac{dQ}{Q} = \nu \frac{a-m}{H} d\chi. \quad (A8)$$

The fact that the characteristics are the trajectories can be

demonstrated by showing that they are identical to the characteristics of the pure advection equation, which are known to be particle trajectories.

The first equality in Equation (A8) leads to:

$$d \ln (\omega + \mu) = -d \ln Q, \quad (A9)$$

which gives linear trajectories in the $(\pi, \theta)$-coordinate system, as obtained in section 2.

The second equality in Equation (A8) gives the rate of change of age with the vertical coordinate, $\omega$,

$$\left. \frac{d\chi}{d\omega} \right|_c = \frac{T}{(\omega + \mu)\nu(\omega)} = \psi \frac{T}{(\omega + \mu)}. \quad (A10)$$

The theory of hyperbolic equations tells us that discontinuities in the boundary condition (e.g. accumulation rate) propagate into discontinuities along these characteristics (in this case in the slope of the isochrones).

MS received 22 December 2004 and accepted in revised form 14 March 2006