Effects of nonlinear rheology, temperature and anisotropy on the relationship between age and depth at ice divides

Ice flow in divide areas is strongly anisotropic. The evolution of ice fabric, from the onset of divide flow towards steady state with a fully developed fabric, has been shown to profoundly affect both the stratigraphy and surface topography of ice divides. Here, we investigate the effects of ice flow on the age-versus-depth relationship at ice divides by using a full Stokes thermomechanical model with a non-linear anisotropic constitutive relation between stress and strain rates. We compare our results with analytical approximations commonly employed in age–depth predictions, such as the Dansgaard and Lliboutry approximations. We show that these approximations systematically underestimate the age of ice at fully developed divides by as much as one order of magnitude. We also show that divides with fully developed fabric are ideal locations for ice-core extraction because ice under them can be up to one order of magnitude older than ice at the same depth at the flanks. In addition, these divides have a distinctive morphological structure that allows them to be clearly identified from satellite imagery or ground-penetrating radar data.


Introduction
Ice cores contain a record of Earth's climate, and are used, for example, to understand how recently observed changes in climate fit within a long history of natural climatic variability.Ideally, an ice-core timeline is obtained by counting annual accumulation cycles (annual layers) or seasonal variations in its chemical constituents, but because of cumulative effects of vertical compression and ice flow distortion these techniques can often not be used in the lowermost sections of an ice core.There ice flow modelling is sometimes the only op-tion available to determine the age of the ice (e.g.Dansgaard et al., 1982;Parrenin et al., 2007).Due to the complexities in modelling flow at ice divides, one-dimensional analytical approximations of ice flow, such as Dansgaard and Johnsen (1969) and Lliboutry (1979), are often used to estimate the age-depth profile at drilling sites.
The underlying difficulty in using flow models to determine the age of ice cores is the complex nature of ice flow.The flow of ice is commonly described by a non-linear rheology known as the Glen flow law (Glen, 1955) where the parameters are poorly constrained by data and are known to be dependent on many factors such as ice impurities, temperature and ice fabric (e.g.Paterson, 1994, ch. 5).For a non-Newtonian fluid such as ice, the effective viscosity is a function of the deviatoric stress state.For ice divides this implies that close to the bedrock, where deviatoric stresses tend to be small, the ice will be harder to deform than in the surrounding areas that are characterized by, in comparison, larger deviatoric stresses.Beneath ice divides this non-linear effect gives rise to the well-known Raymond bump (Raymond, 1983), i.e. an anticline of the internal isochrones.The ice flow in the divide is also known to depend on temperature (Hvidberg, 1996;Nereson and Waddington, 2002), basal sliding (Pettit et al., 2003;Martín et al., 2009b), and along-ridge flow (Martín et al., 2009b).Finally, polycrystalline ice under deformation can develop a preferred crystal orientation fabric and behave as a highly anisotropic material (e.g.Duval et al., 1983).Previous studies have shown that, independently of the model used, anisotropy strongly affects the age distribution around ice divides.The range of models used, includes stationary models with fixed geometry and fabric (e.g.Mangeney et al., 1996), models where fabric is not induced by flow (e.g.Pettit et al., 2007Pettit et al., , 2011)), models that assume stationary geometry (e.g.Seddik et al., 2011), and transient models with flow-induced anisotropy (e.g.Martín et al., 2009a).
In this paper we study ice flow in the proximity of an ice divide and the age-versus-depth relationship.We use a model that is an extension of the model described in Martín et al. (2009a).The improvements in the model, described in Sect.2, allow us to fully explore the effects of the model parameters on the fabric development, and the effects of temperature on ice flow and fabric.The model includes a number of processes not fully accounted for in models commonly used to study age-depth profiles in ice cores.In particular, the model used here includes the effect of flowinduced anisotropy on ice stratigraphy.As shown in Martín et al. (2009a), flow-induced anisotropy has the potential to profoundly affect the position and the shape of isochrones around ice-divides, but the effects of anisotropy on agedepth relationship have, so far, not been described.

Model
In this section we present the numerical model and summarize model equations and boundary conditions.The model builds upon and extends the model used in Martín et al. (2009a).In contrast to Martín et al. (2009a), we include the effect of temperature on flow, we use a better optimised set of parameters to describe the ice fabric, and we use a more appropriate lateral boundary condition based on a zeroth-order anisotropic shallow ice approximation (SIA).

Field equations and boundary conditions
We solve flow in an xz-plane orthogonal to the axis of an ice divide.The z-axis is aligned vertically, the x-axis represent the horizontal direction of flow and in the y direction we assume plane-strain.The ice surface and bed are given by z = s(x, t) and z = b(x), respectively.
The Stokes system and its boundary conditions are Equation (1a) expresses the conservation of mass and Eq.(1b) the conservation of momentum.Here σ is the Cauchy stress tensor, ρ is the density of ice, g is the vertical component of the gravitational acceleration, and v = (u, 0, w) is the velocity vector.Equations ( 1c) and (1d) are the boundary conditions at the surface and bed, respectively.At the lateral boundaries, we assume for the velocity field the anisotropic shallow ice approximation.This approximation is discussed in Sect.2.4.
The kinematic boundary condition at the surface is where a is accumulation rate of ice, expressed as a volume rate per unit area.
The heat equation and boundary conditions are where θ s is the prescribed surface temperature, κ is the thermal diffusivity of ice, c is the specific heat capacity, Q D = trace(S • D) is the dissipation power, D and S are the strain rate and deviatoric stress tensors, K is the thermal conductivity, and Q G is the geothermal heat flux.We assume that the horizontal gradient of temperature is zero at the lateral margins.
Isochrones are lines connecting ice particles with equal age , where age is calculated from

Ice fabric: description and evolution
To describe the ice fabric we use the second and fourth-order orientation tensors, a (2) and a (4) , respectively (e.g.Gödert, 2003;Gillet-Chaulet et al., 2006).We assume that the ice fabric is primarily induced by deformation rather than recrystallisation.In that case, following Gödert (2003) and Gillet-Chaulet et al. (2006) the evolution of the second-order orientation tensor a (2) and boundary conditions can be written as where W is the spin tensor (the skew-symmetric part of the gradient tensor) and α is the interaction parameter.The interaction parameter α controls the relative weighting of the strain rate tensor (D) and the deviatoric stress tensor (S) in the fabric-evolution equation (Eq.5).
The evolution of the second-order orientation tensor a (2) given by Eq. ( 5) depends on a (4) , and for that reason the system of equations listed above is not closed.A common approach of obtaining a closed system is to express the components of the fourth-order orientation tensor as functions of those of the second-order orientation tensor (e.g.Advani and Tucker, 1990).We follow this approach and use the invariantbased closure approximation (IBOF).As shown by Chung and Kwon (2002), the general form of the IBOF closure approximation is (2) ml ) (7) (2) km a (2) ml ) + β 6 Sym(a (2) im a (2) mj a (2) kn a (2) nl ), where "Sym" denotes the symmetrical part of its argument and β i are six functions of the second and third invariants of a (2) .Following Chung and Kwon (2002), we assume that β i are polynomials of degree 5 in the second and third invariant of a (2) , and we use the coefficients computed by Gillet-Chaulet et al. (2006) so that a (4) given by Eq. ( 7) fits the fourth-order orientation tensor given by the orientation distribution function proposed by Gagliardini and Meyssonnier (1999).

Rheology
We assume that the monocrystal grain behaves as a viscous transversely isotropic medium and that there is a uniform stress distribution within the polycrystal (uniform-stress or static model, e.g.Gödert, 2003;Thorsteinsson, 2001;Gagliardini and Meyssonnier, 1999).Following Gödert (2003) and Gillet-Chaulet et al. (2005), we write the orthotropic rheology of the polycrystal as, where I is the identity matrix, the symbols • and : denote the contracted product and the double contracted product, respectively, and the three λ symbols are defined as The mechanical properties of the monocrystal can then be described by the basal shear viscosity η 0 , and the two relative viscosity ratio β, i.e. the ratio of viscosity of the grain for shear parallel to the basal plane to that in the basal plane, and γ , the ratio of the viscosity in compression or tension along the c-axis to that in the basal plane (e.g.Lliboutry, 1987;Meyssonnier and Philip, 1996).The viscosity ratio parameter β is known to be smaller than unity, and the parameter γ to be close to unity (Gillet-Chaulet et al., 2006).
In accordance to the Glen's flow law, we propose a nonlinear extension of the rheology described in Eq. ( 8) where where A(θ ) is the rate factor and it is temperature dependent, n the rheological index and "tr" denotes trace.We use the Dahl-Jensen (1989) relationship for the rate factor with n = 3 and θ c = θ − 273.16 (i.e. is given in • C).Ma et al. ( 2010) and Pettit et al. (2007) propose a similar non-linear extension to Glen's flow law (Eq.9) that depends only on the stress tensor instead of on the strain rate tensor.As observed by Ma et al. (2010), no theoretical or experimental results are available that allow us to discard either of these solutions.We adopt the non-linear extension described in Eq. ( 9) because it reproduces qualitative aspects of the layering in the stratigraphy of a divide evolving towards steady state (see discussion in Sect.4).Further numerical experiments (not reported here), showed that those observations cannot be reproduced using the model described above when using a non-linear extension depending on the stress tensor only.

Anisotropic shallow ice approximation and outflow
boundary conditions Mangeney and Califano (1998) proposed the extension of the SIA (Hutter, 1983) for anisotropic ice.Its use as lateral boundary conditions in full Stokes models is described by Gagliardini and Meyssonnier (2002).The zeroth-order shallow ice expansion approximate Eqs. ( 1a) and (1b) as (Mangeney and Califano, 1998;Gagliardini and Meyssonnier, 2002), where η xzxz is, from Eq. ( 8), the reduced component of the viscosity tensor linking strain rate and stress shear components.

Numerical details
The model presented in the previous sections is similar to that described in Martín et al. (2009a) aside from the coupling between heat equation and ice flow through the viscosity (Eq.9) and the dissipation power (Eq.3a).The numerical algorithms are, however, different.
The system represented by Eqs. ( 1)-( 5) with their respective boundary conditions, is solved with the open-source software Elmer (http:/www.csc.fi/elmer).The Stokes system, Table 1.Numerical values of the parameters used in the simulations.

Parameter Value
heat equation and free surface evolution are solved with the Elmer build-in finite element solvers, using linear elements and stabilised with free residual-free bubbles (see Elmer documentation for details).In contrast, fabric evolution and age equations are solved using a semi-Lagrangian method using a two-time-level scheme and linear interpolation (details in Martín et al., 2009a).The reason for not using a finite element solver in the equations of fabric and age, is that we find the semi-Lagrangian approach more stable close to the base of the divide where the ice is stagnant.

Ice flow analytical approximations
We include here, for reference, the two most commonly employed analytical approximations for ice flow at the ice divide.Dansgaard and Johnsen (1969) approximation for vertical velocity w can be expressed as where z = 0 represent the bedrock, H the ice thickness and z k is a free parameter.Similarly, Lliboutry (1979) approximation for vertical velocity w is where p is a parameter for the vertical profile of deformation that depends on the rheological index n and the gradient of temperature at the base of the divide (see Parrenin and Hindmarsh, 2007, for details).

Evolution of ice fabric and stratigraphy towards the steady state
In order to understand the ice fabric development and its effect on the age-depth relationship, we simulate the evolution of an ice divide from the onset of divide flow towards a steady state.Initial conditions are a flat surface over the whole model domain; ice fabric varying linearly from isotropic at the surface to vertically-aligned single-maximum fabric at the base; temperature following the Robin (1955) analytical approximation; and, for the age of ice, we use the steady state solution for anisotropic SIA.The values of the parameters used in the simulation are listed in Table 1.
A constant surface accumulation rate is prescribed.Ice is allowed to flow out from the left-and right-hand side margins of the domain at a rate that equals the surface influx of ice so that the total ice volume in the domain is constant throughout time.The results are presented in non-dimensional units where both the initial thickness (H ) and the accumulation rate (a) are equal to unity and time is proportional to the characteristic time of the divide (t D = H /a).
As discussed in Martín et al. (2009a), the initial values of fabric and age do not affect the evolution of the divide towards the steady state but they affect the initial stages of divide development (t t D ).Experiments not reported here show that the same applies to temperature.The initial conditions in this study approximate the steady state solutions for flow, temperature and fabric in an area at the flanks of a divide, where ice flow is dominated by shear.Consequently, we The Cryosphere, 6, 1221-1229, 2012 www.the-cryosphere.net/6/1221/2012/are simulating the development of a divide from the onset of divide flow in an area that has been previously dominated by shear.Different initial conditions will affect the solution in the initial stages of divide development (t t D ) but not the evolution towards the steady state.
Figure 1 shows the ice stratigraphy at different stages of the divide development (t = {1/10, 4, 5, 10}t D ). Results are presented for two limiting values of the interaction parameter α.As described above, α = 0 implies a fabric evolution driven only by strain rate, and α = 1 a fabric evolution driven only by stress (see Eq. 6).To reproduce the approximate dimensions of the observed ice stratigraphy (e.g.Martín et al., 2009a), we have, for each given value of α, chosen an optimised value of the relative viscosity ratio β (see Table 1).
Figure 1 shows how, once the ice starts to flow, the isochrones in the lower sections of the ice column start to take the form of anticlines.The amplitudes of the anticlines vary with depth, with the maximum located at around twothirds of the total depth.This feature is known as Raymond bump and is prominent during the first stages of the divide evolution (0 ≤ t D 4t D ).Subsequently, i.e. after about t = 4t D , a syncline starts to develop in the apices of the anticlines giving rise to double-peaked Raymond bumps.At t = 10t D an approximate steady state is reached.
Due to limitations of the model used in Martín et al. (2009a), ice stratigraphy could only be calculated for values of α 1.These limitations have now been addressed and Fig. 1 shows results for both α = 0 and α = 1.As can be seen from inspection comparison of the upper and the lower panels of that figure, despite some quantitative differences, the qualitative features of the ice stratigraphy, i.e. number and general location of synclines and anticlines, are identical.It can be concluded that the results presented in Martín et al. (2009a) for α 1 are in fact generally valid for any value of α.
Further details of the ice fabric and the difference in resulting ice stratigraphy for α = 0 and α = 1 are shown in Fig. 2. In both cases, the fabric varies gradually from being isotropic at the surface, to a vertical girdle fabric in the middle section, and then to a single maximum fabric at the bottom.For a fabric evolution controlled only by the stress tensor (α = 1), the angle between the reference and the orthotropic frame δ is larger and the girdle fabric more pronounced than for a fabric evolution controlled by strain rate tensor (α = 0).The difference between girdle and single maximum is illustrated with the Woodcock (1977) K-value (K < 1 indicate girdle and K > 1 single maximum).

Effects of divide evolution on age-depth and comparison with Dansgaard and Lliboutry analytical approximations
We present here detailed results of the simulation described in Sect.3.1 (Table 1) with parameter α equal to one.Figure 3 compares the vertical velocity and the age variation with depth at the ice divide for different stages of divide development with the space of possible solutions of Dansgaard (Eq.12) and Lliboutry (Eq.13) approximations.The space of solutions include all the possible values of the parameter z k for the Dansgaard approximation (0 ≤ z k ≤ H ) and a wide range of the parameter p (0 ≤ p ≤ 10) for the Lliboutry approximation.
After the onset of flow, in less than one-tenth of t D , the vertical velocity distribution at the divide is just outside the area that represents the space of possible solutions of the analytical approximations.From then onwards there is a gradual transition towards steady state, that it is reached for the vertical velocity after about 3 times t D .The difference between the numerical solutions and the analytical approximations increases with time as the divide evolves towards the steady state.In steady state numerically and analytically calculated vertical velocities differ by as much as 20 % (see Fig. 3).
Regarding the age-depth distribution, there are two stages.First, a gradual increase of the age of ice (t 3t D ) and then a sharp transition towards the steady state.The first stage is characterized by the genesis of a single Raymond bumps in ice stratigraphy and the second by double-peaked Raymond bumps (see Fig. 1).
Figure 4 shows the vertical velocity for the same simulation at a horizontal distance of five ice thicknesses from the divide.Unlike under the divide, the full system solution is within the space of possible solution of the analytical approximations.At this flanking position and for this particular simulation with the parameters described in Table 1, we find that the analytical approximation fits the full system vertical velocity with under 1 % error in velocity for the Dansgaard approximation (best fit z k varying with time from 0.15 H to 0.25 H ) and under 0.1 % for the Lliboutry approximation (best fit p varying with time from 8 to 4, being 4.21 the expected value for a steady state isotropic model).
Figure 4 also shows the transient age-depth relationship at a horizontal distance of five ice thicknesses from the divide.At this position, the transient effects on the age-depth are small (about 2 % in elevation) and within the space of solutions of the analytical approximations.Because the initial age-depth variation is set to that corresponding to the steady state solution of the SIA approximation, we can conclude from this small difference that, at the flanks of the divide area, analytical and SIA approximations predict a similar age-depth to that predicted by the full system.

Discussion
The main unconstrained parameters of the anisotropic rheology model used here (see Sects.2.2 and 2.3), are the stress exponent n, the relative viscosity ratios β and γ , and the interaction parameter α.We assume that γ equals unity (e.g.Gillet-Chaulet et al., 2006).We also fix the stress exponent to 3 and use the temperature-dependent rate factor proposed by Dahl-Jensen (1989) for that value (Eq.10).The effect of the interaction parameter α on modelled stratigraphy is described in Sect.3.1.
The main difficulty in selecting the model parameters n, α, β and γ , is that they are not independent.In Sect.3, we explore two cases that represent extreme values of parameter α: α = 0 and β = 0.01, and α = 1 and β = 0.1.The former corresponds with the model presented in Martín et al. (2009a), the later has been chosen so that the model rheology (Sect.2.3) reproduce the anisotropic properties of ice measured by Pimienta et al. (1987).Their experimental data suggests that single maximum fabric is about ten times softer to shear than isotropic ice.For the rheology discussed in Sect.2.3, that implies For the value n = 3, the above equation gives β = 0.1.In a previous study Martín et al. (2009a) used a value of α equal to zero and a value of β = 10 −2 .Similarly small values of α 1 have been also used for linear rheology by (e.g.Gillet-Chaulet et al., 2006;Durand et al., 2007).As discussed in Drews et al. (2012), the fit between modelled and observed stratigraphy of Halvfarryggen Ridge suggests a value of α close to unity and β ≈ 0.1.Within the context of our model, a value of α so close to unity implies a fabric evolution primarily driven by stress rather than strain rate.It can be argued that using a α value close to unity makes our model approach more consistent, since α = 1 implies that the stress acting on the microscopic crystals and the polycrystal are identical.This is indeed one of the assumptions made in the development of the rheology model we employ (i.e. the uniform stress approximation, see Eq. 8).
The rate factor A is highly sensitive to temperature and for that reason we have used a fully-coupled thermomechanical model in our study.A number of different modelling runs were performed to explore the effect of temperature on the general features of the flow stratiagraphy.We will not report in detail the results of these runs.The main finding is that although temperature can clearly have a strong quantitative effect on the age-depth distribution, the qualitative aspects of the flow (number of synclines and anticlines) do not change.
In particular, we find that the introduction of the temperature variable does not alter the results presented in Martín et al. (2009a) obtained using an isothermal model.In agreement with the results obtained by Hvidberg (1996) for isotropic ice, we find that the coupling between temperature and viscosity produces a zone of relatively warm ice at the base of the divide.The effect of this "warm spot" is to soften the ice and to increase ice flow in this area, with the effect that the ice, for a given depth, is younger than when calculated using an isothermal model.
There are some additional effects that are known to affect fabric and age of ice such as divide migration in a transient stage, recrystallisation, along-ridge flow and basal sliding.Their detailed effect is unclear but generally they tend to reduce the fabric development and the amplitude of the Raymond bump (Martín et al., 2009a).Because of this we expect our age predictions for a given depth to represent an upper limit.
As shown above, once the divide is in steady state, commonly used analytical approximations are unable to reproduce our numerical results (see Sect. 3.2).These analytical approximation can, however, reproduce modelled age-depth distribution in the early stages of divide development.Also, at a horizontal distance of more than 5 times the ice thickness from the divide, the age-depth profiles calculated with our numerical model and the analytical Dansgaard and Lliboutry approximations are similar, differing at all stages of the divide evolution generally by less than 2 %.
Giving our modelling results discussed above the question arises as to what represents the ideal location for an ice-core extraction in ice-divide areas.Particular conditions around the divide, such as non-uniform accumulation (Drews et al., 2012) or thermo-dynamical conditions at the ice-bedrock interface (Seddik et al., 2011), may add elements to the discussion or define the location, but we will focus here on the effects of non-linear flow and flow-induced anisotropy under the divide described in this paper.Ice cores are typically extracted outside but close to an ice dome divide (less than 3-5 times the ice thickness).At this short distance away from the divide, the ice flow regime is arguably simpler to simulate as the effects of anisotropic flow are less pronounced than directly at the divide.From the point of view of ice core interpretation this can be considered to be an advantage.However, ice under the divide could be up to one order of magnitude older than ice at the same depth at the flanks.If the primary aim of an ice-core project is to extract as old an ice as possible, the site should therefore clearly be at the divide and not some distance away from it.
In this context we point out that, as discussed in Martín et al. (2009a), given sufficient time anisotropic flow at ice divides gives rise to distinct morphological surface features.These are commonly clearly visible in satellite imagery as linear features running parallel to the ridge.These linear surface features results from subtle changes in slope in the vicinity of the dome.A further indication of long-term stability, as well as a product of anisotropic flow, is the formation of double-peaked Raymond bumps.Hence, stable ice divides can easily be identified.
These features in satellite imagery and ground-penetrating radar data are widespread in coastal areas of West Antarctica (H. Pritchard, personal communication, 2012).The reason is that the particular conditions in these areas of high snow accumulation and modest ice thickness can produce fully developed divides after as little as a few hundred years of divide position stability (Martín et al., 2009a).

Summary
In this paper, we have improved on the model presented in Martín et al. (2009a) exploring further the influence of fabric, temperature and rheology in ice stratigraphy and age-depth prediction.We have been able to fully explore the effects of the interaction parameters α on the fabric development, and have included the effect of temperature on ice flow and ice fabric.We find that although temperature and the value of the interaction parameter α can clearly have a strong quantitative effect on the age-depth distribution, they do not alter qualitative aspects of the flow (number of synclines and anticlines).
Ice-core dating often rely on ice flow modelling in areas where counting annual layers or seasonal variations of chemical constituents is not possible.We show that the commonly used Dansgaard and Lliboutry approximations systematically underestimate the age of ice at the divide, but are flexible enough to produce results similar to the full system solution when used at a certain distance from the ice divide.
Our study suggests that divides that show double-peaked features in their radar stratigraphy provide ideal locations fore ice-core drilling.The existence of such double-peaked features is indicative of a fully developed ice fabric and suggests stable flow conditions over a period of time comparable to or longer than the characteristic time (ice thickness divided by accumulation).Finally, we have shown that the ice in the summit area can be up to one order of magnitude older than ice at the same depth at the flanks.

Fig. 1 .
Fig. 1.Modelled ice stratigraphy (isochrones) at different stages of the divide development assuming that fabric evolution is driven only by stress (upper panels: α = 1) or strain rate (lower panels: α = 0).Colour represents the age of ice distribution as a multiple of the characteristic time of the divide (t D = H /a). Sub-figures show the solution within an horizontal distance from the divide of five times the initial ice thickness.

Fig. 2 .
Fig.2.Contours of the maximum eigenvalue of the orientation tensor a 33 ( a 33 = 1/3 represents isotropic ice and a 33 = 1 single maximum fabric), the angle in the divide plane between the reference and the orthotropic frame δ and K-Woodcock-value (K < 1 indicate girdle and at K > 1 single maximum) at quasi steady state (t = 10t D ) assuming that fabric evolution is driven by stress (α = 1; top) or strain rate (α = 0; bottom).Sub-figures only show the solution within an horizontal distance from the divide of five times the initial ice thickness.