Hostname: page-component-586b7cd67f-t8hqh Total loading time: 0 Render date: 2024-11-30T15:23:57.225Z Has data issue: false hasContentIssue false

Influence of a non-uniform velocity field on isochrone geometry along a steady flowline of an ice sheet

Published online by Cambridge University Press:  08 September 2017

Frédéric Parrenin
Affiliation:
Laboratoire de Glaciologie et Géophysique de l’Environnement (CNRS–UJF),54 rue Molière, BP 96, 38402 Saint-Martin-d’Hères Cedex, France E-mail: [email protected]
Richard Hindmarsh
Affiliation:
British Antarctic Survey, Natural Environment Research Council, Madingley Road, Cambridge CB3 0ET, UK
Rights & Permissions [Opens in a new window]

Abstract

The relationship between velocity field and isochrone geometry along a steady flowline of an ice sheet is examined. The method is analytical and based upon the stream function and its vertically normalized form, the normalized stream function (NSF). We show that the slope of the isochrones is the slope of the iso-NSF lines, plus a path term which is the cumulative result of the past trajectory of the ice particles. We illustrate this path term in three different examples: varying basal melting, varying basal sliding (Weertman effect) and varying velocity profile around a divide (Raymond effect). The path term generally counteracts the slope of the iso-NSF lines. In the case of the Raymond effect, it can even lead to depressions surrounding the bumps if the transition from dome to flank velocity profile is sufficiently abrupt.

Type
Research Article
Copyright
Copyright © The Author(s) 2007 

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:

(1)

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

(2)

and the bedrock uplift flux as

(3)

See Table 1 for a list of notation used throughout the paper.

Table 1. 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

(4)
(5)

To start with, the horizontal velocity can easily be deduced from the vertical derivative of the partial flux:

(6)

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

(7)

Fig. 1. Illustration of the mass conservation equation on a vertical column truncated at vertical level z.

Taking the limit dx → 0, we obtain

(8)

We can now easily check that

(9)

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

(10)

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

(11)

This new variable can be related to the horizontal flux shape function ω through

(12)

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)

(13)

and with μ being the ratio of the basal flux to the total horizontal flux, i.e.

(14)

Let us now write the velocity field in terms of Ω. From Equation (4), the horizontal velocity ux is given by

(15)

From Equations (5) and (10), the vertical velocity uz can be derived as

(16)

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:

(17)

which simplifies to

(18)

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:

(19)

We now transform from (x, z) to a new system (π, θ) defined by

(20)
(21)

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π = /dt and uθ = /dt. The logarithmic flux horizontal velocity uπ can be derived from the horizontal velocity in the standard x coordinate:

(22)

Using Equations (10) and (15) leads to

(23)

Similarly, we can derive the vertical velocity uθ from u Ω :

(24)

Replacing u Ω by its expression in Equation (19) leads to:

(25)

The parameter κ, defined as either 1/uπ or –1/uθ , will play an important role in the following:

(26)

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.

(27)

which is equivalent to

(28)

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

(29)

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:

(30)

where α, the path-sign parameter, determines the sign of the (π, θ) slope and is given by a simple integration along the particle trajectory

(31)

The slope in the (x, Ω) coordinate is

(32)

Returning to the physical (x, z) coordinate system, the slope is

(33)

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:

(34)

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

(35)

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)

(36)

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

Fig. 2. 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

(37)

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

(38)

Fig. 3. 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:

(39)

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:

(40)

and k is defined

(41)

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.

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

(42)

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.

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

(43)

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.

Fig. 6. 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

(A1)

and

(A2)

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)

(A3)

Fig. 7. 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:

(A4)

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:

(A5)

If we differentiate the SP product with respect to π, we obtain

(A6)

By construction, both particles have the same age and therefore δτ π0 π = δτπ . This gives

(A7)

with

(A8)

which can also be written in terms of (x, Ω) as

(A9)

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

(A10)

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

(A11)

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:

(A12)

Now, dΩ/dx| θ = 0, and from Equations (20) and (21) we obtain

(A13)

and

(A14)

which gives

(A15)

Returning to the physical (x, z) coordinate system, we can write

(A16)

which gives the slope of the isochrones in the (x, z) coordinate system as

(A17)

MS received 16 October 2006 and accepted in revised form 1 April 2007

References

Baldwin, D.J. Bamber, J.L. Payne, A.J. and R.L. Layberry. 2003. Using internal layers from the Greenland ice sheet, identified from radio-echo sounding data, with numerical models. Ann. Glaciol., 37, 325330.CrossRefGoogle Scholar
Conway, H. Hall, B.L. Denton, G.H. Gades, A.M. and E.D. Waddington. 1999. Past and future grounding-line retreat of the West Antarctic ice sheet. Science, 286(5438 280283.CrossRefGoogle ScholarPubMed
Fahnestock, M.A. Abdalati, W. Luo, S. and S. Gogineni. 2001. Internal layer tracing and age–depth–accumulation relationships for the northern Greenland ice sheet. J. Geophys. Res., 106(D24 33,789–33,797.CrossRefGoogle Scholar
Fujita, S. and 6 others. 1999. Nature of radio-echo layering in the Antarctic ice sheet detected by a two-frequency experiment. J. Geophys. Res., 104(B6 13,013–13,024.Google Scholar
Hempel, L. Thyssen, F. Gundestrup, N. Clausen, H.B. and H. Miller. 2000. A comparison of radio-echo sounding data and electrical conductivity of the GRIP ice core. J. Glaciol., 46(154 369374.CrossRefGoogle Scholar
Hindmarsh, R.C.A. 1999. On the numerical computation of temperature in an ice sheet. J. Glaciol., 45(151 568574.CrossRefGoogle Scholar
Hindmarsh, R.C.A. Leysinger Vieli, G.J.M. Raymond, M.J. and G.H. Gudmundsson. 2006. Draping or overriding: the effect of horizontal stress gradients on internal layer architecture in ice sheets. J. Geophys. Res., 111(F2 F02018. (10.1029/2005JF000309.)CrossRefGoogle Scholar
Leysinger Vieli, G.J.M. Hindmarsh, R.C.A. and Siegert., M.J. 2007. Three-dimensional flow influences on radar layer stratigraphy. Ann. Glaciol., 46, 22–28.Google Scholar
Lliboutry, L. 1979. A critical review of analytical approximate solutions for steady state velocities and temperatures in cold ice sheets. Z. Gletscherkd. Glazialgeol., 15(2 135148.Google Scholar
Martín, C. Hindmarsh, R.C.A. and F.J. Navarro. 2006. Dating ice flow change near the flow divide at Roosevelt Island, Antarctica, by using a thermomechanical model to predict radar stratigraphy. J. Geophys. Res., 111(F1 F01011. (10.1029/2005JF000326.)CrossRefGoogle Scholar
Massey, B. 1998. Mechanics of fluids. London, Spon Press.Google Scholar
Nereson, N.A. Raymond, C.F. Jacobel, R.W. and E.D. Waddington. 2000. The accumulation pattern across Siple Dome, West Antarctica, inferred from radar-detected internal layers. J. Glaciol., 46(152 7587.CrossRefGoogle Scholar
Parrenin, F. Hindmarsh, R.C.A. and F. Rémy. 2006. Analytical solutions for the effect of topography, accumulation rate and lateral flow divergence on isochrone layer geometry. J. Glaciol., 52(177 191202.CrossRefGoogle Scholar
Pinglot, J.F. Hagen, J.O. Melvold, K. Eiken, T. and C. Vincent. 2001. A mean net accumulation pattern derived from radioactive layers and radar soundings on Austfonna, Nordaustlandet, Svalbard. J. Glaciol., 47(159 555566.CrossRefGoogle Scholar
Raymond, C.F. 1983. Deformation in the vicinity of ice divides. J. Glaciol., 29(103 357373.CrossRefGoogle Scholar
Reeh, N. and Paterson, W.S.B.. 1988. Application of a flow model to the ice-divide region of Devon Island Ice Cap, Canada. J. Glaciol., 34(116 5563.CrossRefGoogle Scholar
Siegert, M.J. Hindmarsh, R. and G. Hamilton. 2003. Evidence for a large surface ablation zone in central East Antarctica during the last Ice Age. Quat. Res., 59(1 114121.CrossRefGoogle Scholar
Siegert, M.J. and 7 others. 2004. Subglacial Lake Ellsworth: a candidate for in situ exploration in West Antarctica. Geophys. Res. Lett., 31(23 L23403. (10.1029/2004GL021477.)Google Scholar
Vaughan, D.G. Corr, H.F.J. Doake, C.S.M. and E.D. Waddington. 1999. Distortion of isochronous layers in ice revealed by ground-penetrating radar. Nature, 398(6725 323326.CrossRefGoogle Scholar
Weertman, J. 1976. Sliding–no-sliding zone effect and age determination of ice cores. Quat. Res., 6(2 203207.CrossRefGoogle Scholar
Figure 0

Table 1. List of symbols

Figure 1

Fig. 1. Illustration of the mass conservation equation on a vertical column truncated at vertical level z.

Figure 2

Fig. 2. 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.

Figure 3

Fig. 3. 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.

Figure 4

Fig. 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 critical 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.

Figure 5

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

Figure 6

Fig. 6. 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 critical 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.

Figure 7

Fig. 7. Illustration of the derivation of the formula for isochrone slope.