1. Introduction
Radio-echo sounding of the cold Greenland and Antarctic ice sheets has been taking place since the 1960s. Such observations not only measure the ice thickness and basal conditions but also determine internal layering. Such layers are commonly regarded as isochronal layers, i.e. layers of the same age (Reference FujitaFujita and others, 1999; Reference Hempel, Thyssen, Gundestrup and ClausenHempel and others, 2000). Since these layers represent former ice surfaces, they potentially contain significant information about the ice flow.
The connection between isochronal layer geometry and flow fields was originally pointed out by Reference WeertmanWeertman (1976), who described the effect of no-sliding to sliding transitions on these layers. Along such a transition, the increase in basal horizontal velocity induces an increase in the vertical velocity at a given depth. As a consequence, isochrones in the lower part slope towards the bed.
Various other studies established and then employed the connection between ice flow and isochronal layers. At an ice divide, Reference RaymondRaymond (1983) first showed that the non-linear rheology of ice gives rise to lower vertical velocities compared to the flank flow. He suggested that this would result in anticlines in the radar echo layers, which were later observed (Reference Vaughan, Corr and DoakeVaughan and others, 1999). The ‘Raymond effect’ has been exploited to infer past divide-flow changes (Reference Conway, Hall, Denton and GadesConway and others, 1999; Reference Martín and HindmarshMartín and others, 2006).
A third application has been the inference of surface and basal mass-balance conditions. The varying depth of shallow isochrones permits reconstruction of the spatial variation of the accumulation in various areas. For the near-surface isochrones, ice-equivalent layer depth can be viewed as being directly proportional to local accumulation rate (e.g. Reference Pinglot, Hagen, Melvold and EikenPinglot and others, 2001). For slightly deeper layers, one needs to correct for total strain rate using a one-dimensional model (Reference Vaughan, Corr and DoakeVaughan and others, 1999; Reference Fahnestock, Abdalati and LuoFahnestock and others, 2001), and also sometimes for horizontal advection of the ice with two-dimensional or three-dimensional (3-D) models (Nere-son and others, 2000; Reference Baldwin, Bamber and PayneBaldwin and others, 2003; Reference Siegert and HindmarshSiegert and others, 2003). Basal melting rates were inferred from basal isochrones above Lake Ellsworth, West Antarctica (Reference SiegertSiegert and others, 2004).
Despite these uses of isochrone geometry for ice-flow reconstruction, the theoretical relationship between isochrone geometry and ice flow remains generally unclear. In particular, it is unclear how much quantitative information is contained in isochrones: an area where theoretical studies are of use. The steady-state case has been most intensively studied. It is intuitive that isochrones and streamlines should have a fairly intimate relationship. For example, in the case of negligible basal melting, Reference Hindmarsh, Leysinger Vieli and RaymondHindmarsh and others (2006) showed that over short wavelengths, basal isochrones tend to track streamlines. However, we need to extend these restrictions in order to describe the relationship between isochrone geometry and flowlines in the interior of a steady ice sheet. Reference Parrenin and HindmarshParrenin and others (2006) derived an analytical solution for the slope of isochrones, which helps to understand the influence of bedrock and surface geometry, spatial accumulation pattern and lateral flow divergence. However, their analysis only applies when the velocity profiles are spatially uniform as a function of a normalized vertical coordinate, and therefore excludes in practice certain cases where the velocity fields are strictly non-uniform as a result of, for example, changes in basal sliding or ice mechanical behaviour.
The aim of the present study is to derive an analytical expression for the slope of the isochrones along the steady vertical flow tube of an ice sheet. The analysis here applies to cases where the velocity profiles (shape functions) are nonuniform.
To do this, we derive in section 2 the well-known stream function of the flow in the general unsteady case, and introduce the normalized stream function (NSF) Ω which is a vertically normalized form of the stream function. In section 3, we study the steady case with a new coordinate system (π, θ) referred to as the logarithmic flux coordinate system (Reference Parrenin and HindmarshParrenin and others, 2006), in which particles have a linear trajectory. The isochrone slope formula is expressed in section 4, and we show that the slope of the isochrones is closely related to the slope of the iso-Ω lines. It is used to interpret the isochrone geometry in the case of varying basal melting in section 5.1, of no-sliding to sliding transitions (Weertman effect) in section 5.2, or close to ice divides (Raymond effect) in section 5.3.
2. General Considerations in the Unsteady Case
2.1. Basic notations
We consider a flowline of an ice sheet. Time is represented by t. We write the equations in the (x, z) coordinates where x, the horizontal coordinate along the flow, is the distance from the ice divide and z is the vertical coordinate. We suppose that the horizontal direction of the flow does not depend on the vertical position and is time-independent, and we represent the lateral flow divergence by the varying flow tube width Y(x). In practice, for grounded ice, the direction of the flow can be determined from the surface elevation gradient. In this case, the steady assumption for the flow tube means that the shape of surface elevation contours does not change with time. However, the analytical development presented in this paper also applies to other cases (e.g. the flow of ice in an ice shelf).
The ice-sheet geometry is given by the bedrock elevation B(x, t), the surface elevation S(x, t) and the total ice thickness H(x, t) = S(x, t) — B(x, t). We suppose that snow densifies instantaneously, and let the surface accumulation of ice and basal melting rate at the ice–bedrock interface be a(x, t) and m(x, t), respectively. In other words, we implicitly transform the vertical real coordinate into a ‘pure ice equivalent’ vertical coordinate. We let ux (x, z, t) denote the horizontal velocity and uz (x, z, t) the vertical velocity of the ice particles. We also let χ(x, z, t) represent the age of the ice particle.
We now define several fluxes that will be used in section 2.2 to derive the stream function. The partial horizontal flux qH (x, z, t) is defined as the horizontal flux passing below depth z:
with QH(X, t) = qH(x, S, t) being the total horizontal flux at position x and time t. We further define the basal melting flux as
and the bedrock uplift flux as
See Table 1 for a list of notation used throughout the paper.
List of symbols
2.2. The stream function
For clarity, in the following formulae we will neglect the dependencies of the functions on the spatial and temporal coordinates. They will be shown only when it is of special interest.
It is well known from the fluid mechanics theory (Reference MasseyMassey, 1998) that a plane flow with a stationary density ρ can be represented by a stream function ψ. The stream function has an important physical meaning. For example, the flux through any curve linking two points A and B is independent of the chosen curve and is given by the difference ψ B — ψ A. A corollary of this property is that along a streamline, ψ is constant.
We can extend this idea to a flowline of varying width as follows. Here, we can subsume the flow tube width Y(x) into a density, changing in space but not in time, and our ice flow thus derives in principle from a stream function. The aim of this subsection is to derive the stream function of the flow, i.e. a function q verifying
To start with, the horizontal velocity can easily be deduced from the vertical derivative of the partial flux:
We now deduce the vertical velocity profile from the ice mass conservation relationship.
We consider a domain (see Fig. 1) delimited by: (1) two vertical lines at x and x+dx; (2) the bedrock; (3) a horizontal surface at z; and of course (4) the lateral flow tube walls. The statement of integrated mass balance for the domain can be written
Illustration of the mass conservation equation on a vertical column truncated at vertical level z.
Taking the limit dx → 0, we obtain
We can now easily check that
verifies Equations (4) and (5), defining the stream function.
2.3. The (x, Ω) coordinate system
We define the total flux at position x and time t by Q(x, t) = q(x, S, t). Considering mass conservation over an entire ice column leads to
Except at the divide where Q(x, t) = 0, we define the normalized stream function Ω(x, z, t) such that the stream function q is given by
This new variable can be related to the horizontal flux shape function ω through
with ω defined (Reference Reeh and PatersonReeh and Paterson, 1988; Reference Hindmarsh, Leysinger Vieli and RaymondHindmarsh and others, 2006; Reference Parrenin and HindmarshParrenin and others, 2006)
and with μ being the ratio of the basal flux to the total horizontal flux, i.e.
Let us now write the velocity field in terms of Ω. From Equation (4), the horizontal velocity ux is given by
From Equations (5) and (10), the vertical velocity uz can be derived as
We now assume Ω to be a strictly increasing function of z for every (x, t). This corresponds to assuming that there is no reverse flow. It permits the use of Ω as the vertical coordinate, and we will call it the stream vertical coordinate. In the following, we will denote by ZΩ(X,Ω, t) the vertical coordinate z written as a function of (x, Ω), i.e. the elevation of the iso-Ω lines.
The temporal variation of Ω following an ice particle, called here the Ω velocity, is obtained using the chain rule:
which simplifies to
3. The Steady State and the (π, θ) Coordinate System
In this section, we assume steady state and introduce a new logarithmic flux coordinate system (π, θ) in which particles have linear trajectories. This is a similar technique to Reference Parrenin and HindmarshParrenin and others (2006), but with a slightly different definition for π and θ. Apart from being aesthetically pleasant for mathematically minded readers, this coordinate system is considerably easier to work with and to derive the analytical formula for the isochrone slope in section 4.
3.1. The (π, θ) Coordinate system
Since we assume that the ice sheet is in steady state, B(x, t) = B(x), S(x, t) = S(x), H(x, t) = H(x), a(x, t) = a(x), m(x, t) = m(x), Ω(x, z, t) = Ω(x, z), etc. From Equation (18), the Ω velocity becomes:
We now transform from (x, z) to a new system (π, θ) defined by
where Qref = Q(x ref ) and where x ref is a reference position along the flowline. In the following, we refer to π as the logarithmic flux horizontal coordinate and θ as the logarithmic flux vertical coordinate. Note that the change of variable from x to π requires Q(x) to be an increasing function, and we will therefore assume that the accumulation a is strictly positive all along the flowline.
We then define the horizontal and vertical velocity components in this new coordinate system: 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:
Using Equations (10) and (15) leads to
Similarly, we can derive the vertical velocity uθ from u Ω :
Replacing u Ω by its expression in Equation (19) leads to:
The parameter κ, defined as either 1/uπ or –1/uθ , will play an important role in the following:
3.2. Trajectories in (π, θ)
We immediately remark that Equations (23) and (25) sum to zero. This gives simple linear trajectories for ice particles in this logarithmic flux coordinate system (π, θ), i.e.
which is equivalent to
where π 0 = ln(Q0/Q ref ) and Q0 is the total flux at the initial position x0 of the particle at the surface. We have used the notation | c to design the trajectories, reflecting the fact that they are the characteristics of the hyperbolic age equation. We can derive a similar relationship by taking the exponential of Equation (28) to obtain
This equation simply represents the fact that the stream function q is constant along a trajectory in steady state.
4. Slope of Isochrones
We derive in the Appendix an analytical formula for the slope of the isochrones. We first derive the slope in the (π, θ) coordinate system and then transform it to the (x, Ω) and (x, z) coordinate systems.
The (π, θ) isochrone slope is given by:
where α, the path-sign parameter, determines the sign of the (π, θ) slope and is given by a simple integration along the particle trajectory
The slope in the (x, Ω) coordinate is
Returning to the physical (x, z) coordinate system, the slope is
and we remind readers that dz/dx|Ωis the slope of the Ω lines in physical space. Consequently, from Equation (33), there is a close relationship between the slope of the isochrones and the slope of the iso-Ω lines. Isochrone slope is equal to iso-Ω line slope, plus a term:
By comparison to dz/dx|Ω which depends only on the local variation of the Ω field, this P term will be called the path term. Indeed, it depends upon the whole history of the ice particle; it is defined by what happened along its path. This new formulation of the slope of isochrones shows that every spatial variation of the normalized stream function Ω will directly impact the isochrone geometry. In the case of no basal melting, the more restrictive statement of Reference Hindmarsh, Leysinger Vieli and RaymondHindmarsh and others (2006) that over short wavelengths basal isochrones tend to track changes in the velocity profiles is easily interpreted in this framework. This means that near the bed the path term is nearly zero and that dz/dx|Ω is approximately the isochrone slope.
We note here that a change of ice thickness or a change of velocity profile have the same impact on isochrone geometry through the z Ω function. This is a useful way of considering situations where the flow does not follow the bedrock geometry (e.g. at short spatial wavelengths (Reference Hindmarsh, Leysinger Vieli and RaymondHindmarsh and others, 2006) or when there is dead ice near the bed).
Let us look more closely at the path term. It has been already discussed in Reference Parrenin and HindmarshParrenin and others (2006) in the case when the velocity profiles are horizontally uniform. In this study, the impact of spatial variations of surface accumulation rates, ice thickness and flow tube width was discussed. The first multiplicative factor of the path term is the (π, θ) slope α/(1 –α) where α is related to the integration of ∂κ/∂x along the particle trajectory from Equation (31). Along the particle trajectory, a decreasing accumulation a or a positive ∂ 2 z Ω/∂x∂ Ω will rotate the isochrone counterclockwise, and this historical imprint will remain attached to the particle.
The second multiplicative factor of the path term is positive and will have no influence on its sign. It is the ratio of a reference vertical velocity u Ω/(∂Ω/∂z) = aΩ and of the horizontal velocity ux . For this reason, it will be called the velocity scale factor. The greater the horizontal velocity ux and the smaller the reference vertical velocity u Ω/(∂Ω/∂z), the more attenuated the path term. Along a typical particle journey, this velocity scale factor generally decreases, since the sinking of the particle mean Ω decreases and the total flux also increases. Generally, historical imprints are thus diluted with time.
We may remark that the slope formula Equation (33) presents a singularity for x = 0: in this case, both α and ux tend to zero. This uncertainty may be removed by an asymptotic analysis around the divide, but this is beyond the scope of the current study.
How may one visualize the path effect? A convenient way is to look at the isochrones in the (x, Ω) coordinate system, where the spatial variations of Ω disappear (see Equation (32)). One may also look at the isochrones in the (π, θ) coordinate system, which illustrates the path term corrected for the velocity scale factor. The reader will become more familiar with these different visualizations by studying the illustrations in the following section.
5. Illustrations
We illustrate with three different examples, with changes in the horizontal and/or velocity profile, the effect of various parameters on the slope of the isochrones. These examples are not meant to be completely realistic, but rather simplified and pedagogical. They will familiarize the reader with the isochrone slope formula obtained in section 4.
In our examples, the ice thickness is constant along the flow tube and it is convenient to use the normalized vertical coordinate ζ defined by
We will also denote by ωζ the ω variable as a function of (x, ζ) as opposed to (x, z).
5.1. Influence of varying basal melting
In this first experiment, we assume: (1) plug flow; (2) no lateral divergence (Y = constant); (3) ice thickness to be constant and equal to 1000 m; (4) accumulation rate to be constant and equal to 0.2 m a–1 of ice; and (5) a varying basal melting of ice (in m a–1)
with x 1 = 6 km.
Results are plotted in Figure 2. Isochrones sink towards the bed, mainly because iso-Ω lines also sink towards the bed. Indeed, while Ω varies linearly on each vertical (because of plug flow), the value of Ω at bedrock, μ/(1 + μ), increases with the ratio of basal melting flux μ.
Effect of varying basal melting. Particle trajectories are thick black dashed curves. (a) Isochrones in the (x, z) coordinates; (b) iso-Ω lines in the (x, z) coordinates; (c) isochrones in the (x, Ω) coordinates (this diagram is a good illustration of the path effect because iso-Ω lines are flat); and (d) isochrones in the (π,θ) coordinates.
We now consider: what is the path effect associated with this flow? As seen in Figure 2c, this path effect is positive, counteracting the sinking of the iso-Ω lines. Indeed, the vertical space between iso-Ω lines increases when x increases, leading to a positive ∂ 2 z Ω/∂x∂ Ω and thus a positive ∂κ/∂x (since a is constant).
5.2. The Weertman effect revisited
The Weertman effect (Reference WeertmanWeertman, 1976) which occurs at nosliding to sliding transitions and can be observed at the boundaries of subglacial lakes, refers to the fact that in these circumstances the isochrones sink towards the bed. For simplicity, we assume that the ice sheet is in steady state, that the ice thickness is constant and that there is no basal melting. As a consequence, the Ω and ω variables are equivalent. A no-sliding to sliding transition corresponds to a transition from a flow with internal deformation defined by its shape function 
to a plug flow (defined by 
If the change is abrupt, then from Equation (33) we can say that the jump of the isochrones is equal to the jump of the iso-Ω lines. Because 
and 
are equal at Ω = 0 and Ω = 1, the maximum jump is obtained at a coordinate Ω
crit
such that
and we note z crit the corresponding depth before the jump. Following Reference Leysinger Vieli, Hindmarsh and Siegert.Leysinger Vieli and others (2007), we refer to the particle trajectory passing through z crit as the trajectory intersecting elevation of maximum change (TIEMC).
We illustrate the Weertman effect with a simple ice flow model with prescribed velocity profile. We assume: (1) no lateral divergence (Y = constant); (2) accumulation rate to be constant and equal to 0.03 m a–1 ; (3) no basal melting all along the flow tube (μ = 0); (4) bed and surface topography constant and equal to B = 0m and S = 4000 m, respectively; (5) full sliding flow (plug flow) from 40 to 80 km downstream of the ice divide (position x1 and x2, corresponding to π 1 and π 2; see Fig. 3); and (6) flow with internal deformation elsewhere, defined by the horizontal flux shape function
Effect of increased sliding (Weertman effect). The sliding area is situated between 40 km and 80 km. Particle trajectories are thick black dashed curves. (a) Isochrones in the (x, z) coordinates; (b) iso-Ω lines in the (x, z) coordinates; (c) isochrones in the (x, Ω) coordinates; and (d) isochrones in the (π, θ) coordinates.
where p is the exponent of the Glen flow law. Here, we take p = 3. Such a shape function can be derived from the shallow-ice approximation for the case of isothermal and isotropic ice (Reference LliboutryLliboutry, 1979; Reference RaymondRaymond, 1983; Reference HindmarshHindmarsh, 1999).
In reality, and as modelled by Reference Hindmarsh, Leysinger Vieli and RaymondHindmarsh and others (2006, e.g. fig. 7), horizontal stress gradients mean that the shape function changes over distances comparable to the ice thickness. This can correspond to a large or small change in π, depending upon the distance from the divide. In most cases of interest (e.g. the onset of an ice stream) the flux will be large. This implies that the π-transition distance from no-sliding to sliding will be short and it will be appropriate to model the jump as abrupt. At first glance, a is no longer defined since such an abrupt transition in Ω presents a singularity and the derivative ∂κ/∂x is infinite. This is, however, a known property of hyperbolic equations which allows for discontinuities in boundary conditions and in internal data. Here, the integral of ∂κ/∂x through the transition is simply equal to the value of κ after the transition minus its value before. In other words, a is bounded everywhere in the domain. The places where a is discontinuous correspond to the places where the (π, θ) slope of the isochrones is discontinuous.
We now investigate what happens just downstream of the transition, i.e. the path effect. As we can see in Figure 3a, and already deduced by Reference Leysinger Vieli, Hindmarsh and Siegert.Leysinger Vieli and others (2007) from modelling experiments, the near-basal isochrones above the sliding area slope towards the bedrock in the downstream direction, while in the upper part of the ice sheet the isochrones slope towards the surface. This is a typical path effect which can be explained by examination of the path term P. In our modelling experiment, we considered spatially constant conditions of accumulation rate a, melting rate m and surface and bedrock elevations S and B. Assuming that we are in the sliding area, iso-Ω lines are flat and the slope of the isochrones is defined by the path term P, whose sign is that of α. From Equation (31), we find that the sign of α is determined by that of ∂κ/∂x (which is the sign of ∂ 2 z Ω/∂x∂Ω) when the particle crosses the sliding transition. For the particles that cross the transition above the crit ical depth z crit, ∂Z Ω/∂Ω increases (∂ 2ZΩ/∂X∂Ω > 0) and thus the slope of the isochrones is positive. Conversely, for the particles that cross the transition below the crit ical depth, the slope of the isochrones is negative. The two areas of positive and negative path term are delimited by the TIEMC.
5.3. The Raymond effect revisited
‘Raymond bumps’ (Reference RaymondRaymond, 1983) are anticlines in the isochrones found only at ice divides, where the non-linear rheology of ice gives rise to slower vertical velocities. Using the notation of this paper, and assuming no basal melting (Ω = ω), Raymond bumps are characterized by a significantly larger NSF Ω in the vicinity of the dome compared to the flanks. In other words, iso-Ω lines have a negative slope in (x, z). This is the main reason why isochrones have a negative slope, following our analytical Equation (33).
We illustrate the Raymond effect with a simple ice flow model with prescribed velocity profile, based on the Roosevelt Island (Antarctica) example (Reference Martín and HindmarshMartín and others, 2006). We assumed: (1 ) no lateral divergence (Y = constant); (2) accumulation rate to be constant and equal to 0.18 m a–1; (3) no basal melting all along the flow tube; (4) bed and surface topography constant and equal to B = 0 m and S = 750 m, respectively; and (5) flow with internal deformation defined by the following z Ω function:
where 
corresponds to the 
horizontal flux shape function as defined in section 5.2, using an exponent p = 6.5. 
corresponds to the following horizontal flux shape function:
and k is defined
We set x0 = 450 m in our modelling experiment. The values of 
and k are consistent with the more realistic full-Stokes thermomechanical simulation of Reference Martín and HindmarshMartín and others (2006). The results are illustrated in Figure 4.
Effect of reduced vertical velocity at a dome (Raymond effect), with an exponential transition from dome to flank velocity profiles. Particle trajectories are thick black dashed curves. Thick white curve is the crit ical line. (a) Isochrones in the (x, z) coordinates (note the depressions surrounding the Raymond bump in the middle of the ice sheet); (b) iso-Ω lines in the (x, z) coordinates; (c) isochrones in the (x, Ω) coordinates; and (d) isochrones in the (π, θ) coordinates. The abscissa axis is discontinuous for this last panel.
The path effect associated with a Raymond bump is of interest. The sign of the path term 
is the sign of α. In our modelling experiment, the accumulation conditions are constant in the area around the divide and α is given by
As ∂Z
Ω/∂X < 0, there is a
crit
ical Ω
crit
(x) value (corresponding to a
crit
ical depth z
crit
(x)) for which ∂
2ZΩ/∂x∂ω = 0. Above Ω
crit, ∂
2
Z
Ω/∂X∂Ω > 0 and the (π, θ) slope increases. Below Ω
crit, ∂
2ZΩ/∂x∂Ω < 0 and the (π, θ) slope decreases. In our example, as a consequence of the separable form (39), Ω
crit
is constant along the flow tube and corresponds to the value where the difference 
is maximum. With our choices of 
and 
a numerical approximation for Ω
crit
is 0.305.
It is remarkable that the Raymond effect implies a transition of velocity profile which occurs over a non-negligible π distance, contrary to what usually happens in a no-sliding to sliding transition (see section 5.2). This is why the discussion of the path term is more complicated than for this previous example.
For particles in the upper part of the ice sheet, all their journey in the Raymond area is above this
crit
ical depth, and therefore > 0. This means that the path effect has the opposite outcome to the classical Raymond effect (see Fig. 4), rotating isochrones towards the surface. The zone of > 0 actually extends below the
crit
ical depth because of an inertial effect; it takes time for to return to zero.
For some particles deposited sufficiently close to the divide, the first part of their journey above the
crit
ical depth does not impact very much on the term (which is a horizontal integral) because their horizontal displacement is very small. Consequently, the net result of their journey is a negative term, rotating the isochrones toward the bed. The bump at the base of Figure 4 reflects the zone for which < 0, which is delimited by Ω ≈ 0.15 in the dome area in our example.
Finally, there is a transition between those two areas where > 0 and < 0, for which = 0. After leaving the Raymond area, a does not evolve (although the positive scale factor does) and the transition = 0 is passively advected along a trajectory.
The sum of the iso-Ω line slope and the path term slope leads to narrower bumps in the upper part of the ice sheet and broader bumps at the base. In our case, the decay to the flank velocity profile is sufficiently quick that the classical Raymond effect, defined by ∂Z Ω/∂X, does not obliterate the path effect. As a result, depressions appear around the Raymond bump in the middle of the ice sheet. Such depressions have been observed at Roosevelt Island, West Antarctica, (Reference Conway, Hall, Denton and GadesConway and others, 1999; Reference Martín and HindmarshMartín and others, 2006) though only on one flank of the Raymond bump and in the upper parts at Fletcher Promontory (Reference Vaughan, Corr and DoakeVaughan and others, 1999). A deeper example from Fletcher Promontory is displayed in Figure 5.
Example of a trough flanking a Raymond bump from Fletcher Promontory, Ellsworth Land, West Antarctica. Troughs are on right-hand side. Data collected on 5 December 2005 using a 4MHz nominal centre-frequency ground-based ice-penetrating radar.
If the transition from dome to flank velocity profile occurs over a sufficiently broad area, no depressions are visible around the Raymond bump. For example, isochrones for a parameter k defined
are illustrated in Figure 6. The question of which glacio-logical conditions are necessary for the transition to be sufficiently abrupt for the appearance of depressions is still open. For example, do these depressions appear with a Glen isotropic rheology? Or is it necessary to involve complex effects such as sliding just outside the dome area? Such questions could be answered in future studies using full Stokes ice-flow models.
Effect of reduced vertical velocity at a dome (Raymond effect), with a hyperbolic transition from dome to flank velocity profile. Particle trajectories are black thick dashed curves. Thick white curve is the crit ical line. (a) Isochrones in the (x, z) coordinates (note that there are in this case no depressions surrounding the Raymond bump); (b) iso-Ω lines in the (x, z) coordinates; (c) isochrones in the (x,Ω) coordinates; (d) isochrones in the (π, θ) coordinates. The abscissa axis is discontinuous for this last panel.
6. Conclusions and Perspectives
We have studied the link between the velocity field and the isochrone geometry along a steady flowline of an ice sheet using an analytical method. We determined the stream function of the flow and defined the normalized stream function Ω as its vertically normalized form. We further defined a logarithmic flux coordinate system (π, θ), in which ice particles have linear trajectories in the steady case. Using this new coordinate system, we showed that the slope of the isochrones is equal to the slope of the iso-Ω lines, plus a path term, which expresses the relevant past history of the ice layer along its path.
We illustrated this path term in the case of varying basal melting, varying basal sliding (Weertman effect) and varying velocity profile at the divide (Raymond effect). We show that in those examples the path term generally counteracts the slope of the iso-Ω lines. In the case of the Raymond effect, it can even lead to depressions on either side of the bumps if the transition from dome to flank velocity profile is sufficiently abrupt.
This analytical development opens new perspectives to reconstruct the velocity field in a steady-state ice sheet from measured isochrones. We may attempt, for example, to remove the path effect from the slope of the isochrones to derive the normalized stream function in the domain. The feasibility of such an inverse method and its application to radar profiles from Greenland and Antarctica should be studied in future works.
A further interesting question is how the ice flow can be constrained by isochronal layers in the full 3-D case. We do not expect such simple analytical solutions as in the present study, as there is no simple generalization of the stream function. The inversion of the flow from isochronal layers will in this general case be an ill-posed problem, and additional information from mechanical equations will be needed to reconstruct the velocity field.
Acknowledgements
This work is part of the French project MIDIGA (ANR) (Modeling and Identification for Drillings Interpretation in Greenland and America (Agence National de la Recherche)). The paper benefited from comments by two reviewers: A. Salamatin provided a very thorough and constructive review and encouraged us to simplify some of the manuscript while G.K.C. Clarke challenged us to deliver something accessible to a broader spectrum of readers. We also thank C. Martin for giving access to his modelling results for Roosevelt Island, E. King for filtering the radargram, and E. King and H. Corr for assistance in preparing the radar.
Appendix
Slope of Isochrones
We derive an analytical formula for the slope of the isochrones. The slope depends on the spatial coordinate system, and will be qualified according to, for example, the (π, θ) slope or the (x, z) slope. We first derive the (π, θ) slope and then convert this formula to the (x, Ω) and (x, z) coordinate systems.
In the following, P and S are two functions of (π, θ), defined as
and
A1. Slope in the (π, θ) coordinate system
Consider two particles initially at the upper surface of the ice sheet at arbitrary neighbouring positions x0, x0 + δx0 (which correspond to π 0, π 0 + δπ 0). After time Δt the particles will have positions x, x + δx (which correspond to π, π + δπ) with elevations z, z + δz (with corresponding θ, θ + δθ). Our goal is to derive an expression for the isochrone (x, z) slope by first deriving the (π, θ) slope. As indicated by the notation, we will assume that δπ, etc., are small quantities. The proof which follows is similar to that in Reference Parrenin and HindmarshParrenin and others (2006) (replacing T and ψ by S and P, respectively), with the exception that (1) P is now a function of π and θ and not just θ; and (2) the product SP is differentiated with respect to π and not with respect to θ.
As a consequence of Equation (27), the trajectories in the (π, θ) space of both particles are two parallel lines with a slope of – 1, separated by the horizontal distance δπ 0 (Fig. 7). This is the fundamental reason why we can easily derive an analytical expression for the slope of the isochrones. The (π, θ) slope of the isochrones can be written (assuming that δπ0 → 0)
Illustration of the derivation of the formula for isochrone slope.
(The notation | x is used throughout the paper to define constant age, i.e. an isochronal layer.)
To proceed further, we must relate δπ to δπ 0. Firstly, we need to evaluate the normalized time δτπ required for the second particle to travel from π + δπ 0 to π + δπ (see Fig. 7). When δπ, δπ 0 are small, these quantities can be derived from Equation (23) as:
Secondly, we need to evaluate the difference in duration between the two journeys of the particles, represented by δτ π0→ π , from π 0 to π for the first particle and from π 0 + δπ 0 to π + δπ 0 for the second particle. Using Equation (23), we can determine δτ π0 →π to be:
If we differentiate the SP product with respect to π, we obtain
By construction, both particles have the same age and therefore δτ π0 →π = δτπ . This gives
with
which can also be written in terms of (x, Ω) as
Finally, we replace δπ by (1 – α)δπ 0 in Equation (A3) to reach the first goal of deriving an expression for the (π, θ) isochrone slope:
The (π, θ) slope therefore depends on the integration along the particle trajectory of ∂(SP)/∂π. The initial position of the particle can be evaluated from Equation (29), and thus the slope of the isochrone can be estimated directly by integrating Equation (A8).
We can express α ina different way through an integration by parts of Equation (A8):
with S 0 = S(π 0) and P 0 = P(π 0, θ = 0). From Equation (A7), it follows that α > 1 corresponds to the fact that isochrone slope exceeded the vertical, which is not possible for a vertically increasing horizontal velocity. Indeed, such a velocity field implies ∂ 2Ω/∂z2 positive for all x, ∂P/∂θ negative and thus (1 – α) positive from Equation (A11). Hence, we will refer to a as the path-sign parameter as in general its sign determines the sign of the isochrone slope in the (π, θ) coordinate system.
A2. Slope of isochrones in (x, Ω) and (x, z)
We can now write the slope of isochrones in the (x, Ω) coordinate system as:
Now, dΩ/dx| θ = 0, and from Equations (20) and (21) we obtain
and
which gives
Returning to the physical (x, z) coordinate system, we can write
which gives the slope of the isochrones in the (x, z) coordinate system as
MS received 16 October 2006 and accepted in revised form 1 April 2007
