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

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:

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 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: 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±, 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:

(1)

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:

(2)

where ζ =(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 ω 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:

(3)

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:

(4)

Multiplying Equation (3) by Y and differentiating with respect to x leads to:

(5)

and the definition ζ = (z – B)/H leads to:

(6)

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

(7)

Now

(8)

and therfore

(9)

Finally, replacing ∂ζ/∂x and ∂YŪ_x/∂x in Equation (5) using Equations (6), (7) and (9) leads to:

(10)

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

(11)

Now u_z(z = B) is given by:

(12)

where s = ω′(ζ = 0) is the sliding rate (s = 0 if no sliding and s = 1 if full sliding, i.e. plug flow). Therefore

(13)

An integration by parts gives:

(14)

leading finally, with two cancellations, to

(15)

The vertical velocity in the ζ coordinate, uζ = dζ/dt, can be written in terms of the vertical velocity in the z coordinate, u_z , using the chain rule, as:

(16)

leading to

(17)

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:

(18)

we can write

(19)

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:

(20)

(21)

where Q_ref is a reference flux at a reference position x_ref. 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,

(22)

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

(23)

where ψ(θ) is the inverse of the horizontal-velocity shape function ω′ written in the θ-coordinate,

(24)

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, ζ,

(25)

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

(26)

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 (π, θ)

(27)

which is equivalent to:

(28)

where π₀ = ln(Q₀/Q_ref) and Q₀ is the total flux at the initial position, x ₀, 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:

(29)

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)

(30)

Similarly, the age, χ, of a particle at a given position (π, θ) is readily computed from

(31)

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:

(32)

(33)

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

(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 (π, θ) 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 ₀, x ₀ + δx₀ (corresponding to π₀, π₀ + δπ₀). 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 = δπ₀, and we have δπ = δθ_traj — δθ (see Fig. 1). This leads to the following expression for the slope of the isochrones (assuming that δπ0 → 0):

Fig. 1. Scheme for the derivation of the formula for isochrone slope.

(35)

The notation |_χ is used throughout this paper to define constant age, i.e. an isochrone layer. To proceed further, we must relate δπ to δπ₀. Firstly, we need to evaluate the time, δt _π0, required for the first particle to travel from π₀ to ҡ₀ + δҡ₀ and the time, δt _π, for the second particle to travel from π to π + δπ (see Fig. 1). When δπ, δπ₀ are small, these quantities can be derived from Equation (23) as:

(36)

(37)

with T₀ = T(π₀ ) and ψ ₀ = (θ = 0).

Secondly, we need to evaluate the difference in duration between the two journeys of the particles from π₀ + δπ₀ 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:

(38)

or, if we differentiate the ψ function,

(39)

Using the fact that δθ_traj is constant and equal to –δπ₀ as illustrated in Figure 1, we find:

(40)

or, replacing θ by π₀ — π by Equation (28) and taking the limit as δπ ₀ → 0,

(41)

where

(42)

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:

(43)

and replacing δt _π0 and δt _π by Equations (36) and (37), we get:

(44)

or

(45)

with

(46)

Finally we replace δπ by (1 – α)δπ₀ in Equation (35) to reach the first goal of deriving an expression for the (π, θ) isochrone slope:

(47)

The (π, θ) slope depends on the value of the characteristic time T = H/(a – m) at the initial and final positions of the particle, π₀ and π, and depends also on the value of T between π₀ 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,

(48)

In this formula, note that Q/Y represents the local flux

(49)

Now, since we can write the isochrone slope as:

(50)

with γ(x) defined by:

(51)

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:

(52)

leading, with Equation (46), to

(53)

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

(54)

and the (π, θ) and (x, ζ) slopes, respectively, become

(55)

and

(56)

or

(57)

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 ₁ = 4000 m from the divide to 30 km downstream (position x = x1, corresponding to π = π1; see Fig. 2), then decreasing abruptly to H₂ = 2000 m, and keeping this value downstream. The ice surface is flat.

Fig. 2. Isochrones over a step change in bedrock elevation at x = x ₁ (corresponding to π = π₁), 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ļ/H₂ ) gives the (π, θ) slope of the isochrones as a constant of–1/2. Equation (57) gives the (x,ζ) slope of the isochrones as:

(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 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 = x₂, corresponding to π = π₂; see Fig. 3), creating a plateau in the bed.

Fig. 3. Isochrones over a bedrock plateau delimited at x = x ₁ and x = x₂ (corresponding to π = π₁ and π = π₂). 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 – (H₂/H₁), 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:

(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, ω, which permits internal deformation,

(60)

(61)

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:

(62)

leading to

(63)

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:

(64)

The slope of the isochrones in the logarithmic-flux coordinate system (π, θ) is given by Equation (47), with α given by Equation (53),

(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, ∂T/∂π = 0 between initial and final horizontal positions (π0 and π), except for the bedrock step ∆B at π = π₁ which leads to a step ∆₁ T = T ₂ – T ₁ in the T parameter. Thus, replacing π₀ using Equation (28), we obtain:

Fig. 4. Isochrones over a bedrock plateau delimited at x = x ₁ and x = x₂ (corresponding to π = π-i and π = π₂). See section 4.3 for a complete description of the experiment. Ice is experiencing internal deformation. Details as for Figure 2.

(66)

and

(67)

giving the following expression for the slope of the isochrones:

(68)

Immediately downstream of the bedrock step at π = π₁, 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:

(69)

where π2 is the horizontal coordinate of the second bedrock step (see Fig. 4). Again, the (π, θ) slope is 1 just after the bedrock step π₂ and decreases toward 1/2 as π increases. For area A4, one can also show, from Equation (65), that:

(70)

with ∆₂ T = T ₃ – T ₂. This leads to

(71)

giving the following formula for the (π, θ) slope of the isochrones:

(72)

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 (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:

(73)

(74)

(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:

(75)

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 (Schøtt 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, µ = 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 (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, !, 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,

(76)

we obtain

(77)

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:

(78)

(79)

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, ω:

(80)

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 (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.