Notation
Symbols used in this paper, listed in the order in which they are defined:
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 (Reference Bailey, Evans and RobinBailey and others, 1964; Reference Robin, Evans and BaileyRobin 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. Reference Robin, Drewryand and MeldrumRobin 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 Reference Jacobel, Scambos, Raymond and GadesJacobel and others (1996) to conclude that there had been an abrupt stoppage of the ice stream a few centuries ago. Reference Ng and ConwayNg and Conway (2004) have deduced that the stagnant Kamb Ice Stream flowed at velocities of over 350 ma–1, from the internal layer stratigraphy.
A second application has been the stratigraphic correlation of ice cores, allowing the comparison of their respective chronologies: Reference Jacobel and HodgeJacobel and Hodge (1995) linked GRIP and GISP2 in Greenland, while Reference Siegert, Hodgkins and DowdeswellSiegert 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 (Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002) and Ellsworth subglacial lake (Reference SiegertSiegert 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. Reference Pinglot, Hagen, Melvold, Eiken and VincentPinglot and others, 2001). For slightly deeper layers, one needs to correct for total strain rate with a one-dimensional (1-D) model (Reference Vaughan, Doake and WaddingtonVaughan and others, 1999; Reference Fahnestock, Abdalati, Luo and GogineniFahnestock and others, 2001), and also sometimes for horizontal advection of the ice with two- or three-dimensional models (Reference Nereson, Raymond, Jacobel and WaddingtonNereson and others, 2000; Reference Baldwin, Bamber, Payne and LayberryBaldwin and others, 2003; Reference Siegert, Hindmarsh and HamiltonSiegert and others, 2003).
A fourth application has been the extraction of past divide location histories exploiting the Raymond effect (Reference RaymondRaymond, 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 (Reference Conway, Hall, Denton, Gades and WaddingtonConway 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 (Reference WeertmanWeertman, 1976; Reference WhillansWhillans, 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 uz ((x±, 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:
with Q(x) = q(x, S) being the total horizontal flux.
We now define the flux shape function, ω, as the function defined on [0;1] such that the partial flux q is given by:
where ζ =(z — B)/H is the normalized vertical coordinate. Such shape functions have been used previously in glaciology (Reference LliboutryLliboutry, 1979; Reference ReehReeh, 1988; Reference RitzRitz, 1989). We assume here that ω 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)(∂q/∂z), leading to:
where Ūx = Q/YH is the average horizontal velocity and ω′ is the derivative of ω with respect to ζ. In other words, ω′ 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:
Multiplying Equation (3) by Y and differentiating with respect to x leads to:
and the definition ζ = (z – B)/H leads to:
Moreover, the ice mass conservation applied to a whole ice column leads to:
Now
and therfore
Finally, replacing ∂ζ/∂x and ∂YŪx/∂x in Equation (5) using Equations (6), (7) and (9) leads to:
From Equation (4), uz can be written as:
Now uz(z = B) is given by:
where s = ω′(ζ = 0) is the sliding rate (s = 0 if no sliding and s = 1 if full sliding, i.e. plug flow). Therefore
An integration by parts gives:
leading finally, with two cancellations, to
The vertical velocity in the ζ coordinate, uζ = dζ/dt, can be written in terms of the vertical velocity in the z coordinate, uz , using the chain rule, as:
leading to
with ã = a/H and 
the normalized accumulation and melting rates. Notice that this form removes the dependence of the vertical velocity in ζ coordinate on the basal and surface slopes.
Defining the ratio μ, relating basal melting to mass balance, as:
we can write
where T(x) = H(x)/(a(x) — m(x)) is called the local characteristic time. We assume here that μ is uniform all along the flow-tube, which results in a spatially uniform shape function for the vertical velocity, uζ. To further the analysis we change the coordinate system from (x, ζ) to a new system (π, θ) defined by:
where Qref is a reference flux at a reference position xref. In the following, we will call π the logarithmic-flux horizontal coordinate and θ the logarithmic-flux vertical coordinate. Note that the change of variable from x- to π-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 ζ- to θ-coordinate requires ω(ζ) to be a strictly increasing function, which is the case when ω′ is a strictly positive function of ζ.
We then define the horizontal and vertical velocity components in this new coordinate system as u π = dπ/dt and uθ = dθ/dt. The logarithmic-flux horizontal velocity u π can be derived from the horizontal velocity in the standard x-coordinate,
Replacing Q using Equation (8), ∂Q/∂x using Equation (7) and dx/dt using Equation (3) leads to:
where ψ(θ) is the inverse of the horizontal-velocity shape function ω′ written in the θ-coordinate,
and is called the time shape function. We have used a different notation here for 1/ω' to emphasize that it is the derivative of 1/ω′ with respect to ζ, and not with respect to θ.
Similarly, we can derive the vertical velocity in the logarithmic-flux vertical coordinate, θ, from the corresponding quantity in the normalized vertical coordinate, ζ,
and replacing dζ/dt by its expression in Equation (17) leads to the following expression:
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 (π, θ)
which is equivalent to:
where π0 = ln(Q0/Qref) and Q0 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π/dθ)|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:
where x 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 (π, θ) 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)
Similarly, the age, χ, of a particle at a given position (π, θ) is readily computed from
For plug flow, by definition ω′ ≡ 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:
with Equation (27) also applying to plug flow. The formula for the age simplifies to:
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 (π, θ) slope, the (x, ζ) slope or the (x, z) slope. We will focus in particular on the (x, ζ) 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 + δx0 (corresponding to π0, π0 + δπ0). After time ∆t they will have positions x, x + δx (corresponding to π, π + δπ) with elevations ζ, ζ + δζ (corresponding to θ, θ + δθ). Our goal is to derive an expression for the isochrone (x, ζ) slope, but we will first derive the (π, θ) slope. As indicated by the notation, we will assume that δπ, etc. are small quantities.
As a consequence of Equation (27), the trajectories in the (π, θ) space of both particles are two parallel lines with a slope of –1. The vertical shift of the two lines is equal to δθtraj = δπ0, and we have δπ = δθtraj — δθ (see Fig. 1). This leads to the following expression for the slope of the isochrones (assuming that δπ0 → 0):
Scheme for the derivation of the formula for isochrone slope.
The notation |χ is used throughout this paper to define constant age, i.e. an isochrone layer. To proceed further, we must relate δπ to δπ0. Firstly, we need to evaluate the time, δt π0, required for the first particle to travel from π0 to ҡ0 + δҡ0 and the time, δt π, for the second particle to travel from π to π + δπ (see Fig. 1). When δπ, δπ0 are small, these quantities can be derived from Equation (23) as:
with T0 = T(π0 ) and ψ 0 = (θ = 0).
Secondly, we need to evaluate the difference in duration between the two journeys of the particles from π0 + δπ0 to π, represented by δt π0→π. In plug flow, δt π0→π = 0 because the horizontal velocity does not depend on the vertical position of the particle. In the general case, we can determine δt π0→π using Equation (23) to be:
or, if we differentiate the ψ function,
Using the fact that δθtraj is constant and equal to –δπ0 as illustrated in Figure 1, we find:
or, replacing θ by π0 — π by Equation (28) and taking the limit as δπ 0 → 0,
where
Note for plug flow that ß = 0 as ψ ≡ 1. This is why we call ß the non-plug-flow parameter.
Now, by construction, both particles have the same age:
and replacing δt π0 and δt π by Equations (36) and (37), we get:
or
with
Finally we replace δπ by (1 – α)δπ0 in Equation (35) to reach the first goal of deriving an expression for the (π, θ) isochrone slope:
The (π, θ) slope depends on the value of the characteristic time T = H/(a – m) at the initial and final positions of the particle, π0 and π, and depends also on the value of T between π0 and π through the function ß. 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 → π and ζ → θ are monotonic, the signs of the (x, ζ) and (π, θ) slopes are the same. Moreover, assuming that ω” is positive, ß is negative and thus, from Equation (46), (1 – α) is positive. We will thus call α the slope-sign parameter, because its sign determines the sign of the isochrone slope.
From Equations (20) and (21), we see that dπ = dQ/Q, dθ = dω/(ω + µ) and we also use the identity dω = ω′dζ. With dQ = (a – m) Y dx from Equation (7), we deduce the (x, ζ) slope,
In this formula, note that Q/Y represents the local flux 
Now, since 
we can write the isochrone slope as:
with γ(x) defined by:
Note that γ(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 α, which illustrates different aspects of the problem, can be derived as follows. Integrating Equation (42) by parts, ß can be expressed as:
leading, with Equation (46), to
The interest of this formula is that the new integral term in α contains the derivative of T with respect to π, instead of the derivative of ψ with respect to θ. It can now be seen immediately that the (x, ζ) 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 ζ-flat to qualify these isochrones, because it is equivalent to having a zero slope in the (x, ζ)-, (Q, ζ)- or (π, ζ)-coordinate systems. Similarly, we will use z-flat for isochrones with a zero slope in the (x, z)-, (Q, z)- or (π, z)-coordinate systems.
3.2. Plug flow as a special case
In plug flow, ω = ζ, ω′ ≡ 1 and ψ ≡ 1, lead to ß = 0. Thus
and the (π, θ) and (x, ζ) slopes, respectively, become
and
or
The remarkable aspect of Equation (56) is that the (x, ζ) slope depends only upon the elevation, the initial (through α) and final positions of the particle, and the local flux 
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(π) 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, ζ) 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 ζ. Indeed, we can see that Equation (17) still holds when, at some points, ∂B/∂x and ∂H/∂x tend to infinity and the a and 
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, ω, 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 π (directly deduced from Q) and the x-dependent scale factor (also by an integration from the ice divide). Second, for every ζ level, we can compute the logarithmic-flux vertical variable, θ, the ω′ value and the depth-dependent scale factor. Third, for every (x, ζ) node, we can compute the spatial origin of the particle from Equation (28) or (30), the ß term from the convolution-type Equation (42) (by a simple numeric integration against the horizontal π variable) and the principal term α/(α — 1) (deduced directly from ß). 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 = constant), (2) accumulation rate to be constant and equal to 3 cm a–1 of ice, (3) no basal melting all along the flow-tube (µ = 0) and (4) ice thickness equal to H 1 = 4000 m from the divide to 30 km downstream (position x = x1, corresponding to π = π1; see Fig. 2), then decreasing abruptly to H2 = 2000 m, and keeping this value downstream. The ice surface is flat.
Isochrones over a step change in bedrock elevation at x = x 1 (corresponding to π = π1), marked by a thick grey vertical line. Ice is in plug flow. (a) (x, z)-coordinate system. (b) (π, θ)-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).
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 α = 1 – (Hļ/H2 ) gives the (π, θ) slope of the isochrones as a constant of–1/2. Equation (57) gives the (x,ζ) slope of the isochrones as:
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 (Reference Rubinstein and RubinsteinRubinstein 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 nonzero (x, ζ) 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 0m at 60 km (position x = x2, corresponding to π = π2; see Fig. 3), creating a plateau in the bed.
Isochrones over a bedrock plateau delimited at x = x 1 and x = x2 (corresponding to π = π1 and π = π2). See section 4.2 for a complete description of the experiment. Ice is in plug flow. Details as for Figure 2.
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), α = 0 and the isochrones are consequently ζ-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 ζ-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 α = 1 – (H2/H1), we find that the slope of the isochrones in the logarithmic-flux coordinate system (π,θ) is constant and equal to 1, and Equation (57) gives the slope of the isochrones in the (x, ζ)-coordinate system as:
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, ω, which permits internal deformation,
with p = 1.5 chosen in our example. We call p the deformation exponent. For 1 ≤ p ≤ 2, these simulate typical velocity profiles for simple shear. This shape function is qualitatively similar to standard glaciological forms:
leading to
but is analytically more tractable. A difference that requires comment is the fact that ω″ (which is a multiple of the shear rate) presents a singularity at ζ = 0 (it tends to +∞ at ζ = 0), while ω″ was used during the analytical development in section 2. This singularity is not important, because only the non-singular integrated forms 
or 
appear in the analytical development.
Using our additional assumption that µ = 0, it is straightforward to show that ψ is given by:
The slope of the isochrones in the logarithmic-flux coordinate system (π, θ) is given by Equation (47), with α given by Equation (53),
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, ∂T/∂π = 0 between initial and final horizontal positions (π0 and π), except for the bedrock step ∆B at π = π1 which leads to a step ∆1 T = T 2 – T 1 in the T parameter. Thus, replacing π0 using Equation (28), we obtain:
Isochrones over a bedrock plateau delimited at x = x 1 and x = x2 (corresponding to π = π-i and π = π2). See section 4.3 for a complete description of the experiment. Ice is experiencing internal deformation. Details as for Figure 2.
and
giving the following expression for the slope of the isochrones:
Immediately downstream of the bedrock step at π = π1, the (π, θ) slope is –1/2 and decreases with π (with a minimum value of –1).
Similarly, one can show that in area A5 the (π, θ) slope of the isochrones is given by:
where π2 is the horizontal coordinate of the second bedrock step (see Fig. 4). Again, the (π, θ) slope is 1 just after the bedrock step π2 and decreases toward 1/2 as π increases. For area A4, one can also show, from Equation (65), that:
with ∆2 T = T 3 – T 2. This leads to
giving the following formula for the (π, θ) slope of the isochrones:
Unlike the previous example, isochrones are not ζ-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, ω, 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, ζ, eliminates this primary effect. In other words, the (x, ζ) slope informs us about the processes governing ice flow.
From Equation (50), the (x, ζ) slope is the product of three terms: one principal term, and two scale factors. The scale factors can modify the amplitude of the (x, ζ) slope but not its sign.
One scale factor, γ(x), depends upon the horizontal position of the particle, and the other, (ω(ζ)+ µ)/ω′(ζ), depends upon the normalized elevation, ζ, 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, !. For plug flow, it simply increases linearly from µ at the bedrock to 1 + µ at the surface. When there is no melting, the depth-dependent scale factor tends to zero as ζ 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, α/(1 – α), may tend to +∞. When there is melting but no sliding, the depth-dependent scale factor tends to +∞. That does not mean, however, that the isochrones approach a vertical orientation as they approach the bedrock interface, because the slope-sign term, α, tends to zero.
The principal term α/(1 – α) 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, ζ) 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, ζ) isochrone slope is positive and, conversely, if it is smaller the (x, ζ) slope is negative. For general flow, the basic pattern remains the same, but the values of the characteristic time in between enter in the α parameter in Equation (46) through ß, 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 (Reference Ritz, Rommelaere and DumasRitz and others, 2001), i.e. less than 10%. There are hints from numerical calculations (Reference Nereson, Raymond, Jacobel and WaddingtonNereson 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:
(with R(t) a non-negative multiplication factor) we can come back to the steady-state problem with melting rate 
and accumulation rate ā(x) by changing the time variable from t to 
defined by:
Because of the bijection between t and 
(for every value of t there is one unique value of 
) isochrones in t are equivalent to isochrones in 
. The true age, χ, can simply be deduced from the normalized age, χ, with Equation (73). This technique to include transient accumulation rates has been used by previous workers (Reference Schøtt, Waddington and RaymondSchøtt and others, 1992; Reference Nereson, Waddington, Raymond and JacobsonNereson 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, µ = 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 µ 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, !. 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 (Reference LliboutryLliboutry, 1979; Reference HutterHutter, 1983). More discussions on this can be found in Reference Hindmarsh, Straughan, Greve, Ehrentraut and WangHindmarsh (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, !, 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 (π, θ) slope of –1. The relationship between the logarithmic-flux vertical coordinate, θ, and the normalized vertical coordinate, ζ, 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 π, 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 ω(ζ = 1). However, if we suppose that the flux shape function, ω, is uniform, gradients in surface velocity still give spatial variations in accumulation, which is sufficient to compute the isochrone slope and the x → π change of coordinate.
To transform the vertical coordinate from ζ to θ, we need to know the flux shape function, ω, and the ratio (relating basal melting to mass balance), µ. Assuming that we restrict our consideration to an area where there is no melting, and taking the flux shape function, ω, from section 4.3,
we obtain
Because the flux shape function ω′ is unknown, the measured isochrone layers cannot be plotted in the (π, θ)-coordinate system. As an alternative, we can use a trial coordinate 
as the vertical coordinate. In (π, 
) space, particle trajectories now have a slope of 1 /p, giving an estimate of the deformation parameter p.
More generally, i.e. with an arbitary flux shape function, ω, we can write:
and thus, if we define S(ζ) as the slope of the trajectories in (π, θ) space, writing dω/dζ = (dω/dθ)(dθ/dπ)(dπ/dθ)(dθ/ dζ) leads to the following differential equation for the flux shape function, ω:
In this case, it is therefore, in principle, possible to deduce the flux shape function, ω, from a discontinuity in the slope of the isochrones in the (π, 
) 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 (π, θ) 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, ζ)- coordinate system (x being a horizontal, along-flow coordinate and ζ being the normalized vertical coordinate). The (x, ζ) 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, ζ) 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, ζ) 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, ζ) 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, ζ) 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 (Reference WeertmanWeertman, 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’.
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, Reference HindmarshHindmarsh (1999) shows that the Lagrangian derivation operator D(·)/Dt in ζ coordinates can be written as:
where ∇H is the horizontal gradient in ζ coordinate and (·) represents the field being advected. In steady-state conditions and with a uniform flux shape function, ω,
The age equation Dχ/Dt = 1 may therefore be written as:
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
and
where
This hyperbolic equation has characteristic equations
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:
which gives linear trajectories in the (π, ớ)-coordinate system, as obtained in section 2.
The second equality in Equation (A8) gives the rate of change of age with the vertical coordinate, ω,
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).
