Anisotropic law and fabric evolution

In this section, we briefly present the anisotropic ice-flow model used in this study. For more details the reader can refer to previous publications related to the model itself (Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and CastelnauGillet-Chaulet and others, 2005, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006; Reference DurandDurand and others, 2007; Reference Martín, Hindmarsh and NavarroMartín and others, 2009). The fabric is described using the concept of orientation tensors. The second- and fourth-order orientation tensors are defined as

where c is the c-axis unit vector and () denotes the average over all the grains that make up the polycrystal. In this continuum description of the fabric, the polycrystal represents the local behaviour of a representative elementary ice volume. By assuming that the fourth-order orientation tensor, a⁽⁴⁾, is given as a tensorial function of a⁽²⁾ (Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and RuokolainenGillet-Chaulet and others, 2006), the fabric can be described in a very condensed way using only a⁽²⁾. By definition, tr a⁽²⁾ = 1, so that only the first two eigenvalues, a⁽²⁾ and a|,²⁾, and three Euler angles are needed to completely define the fabric. As a consequence, modelled fabrics are orthotropic, i.e. the c-axis distribution presents three orthogonal symmetry planes. Although the orthotropy is a simple form of the most general anisotropy, it is thought to be a good compromise between physical adequateness and simplicity. The second- order orientation tensor allows us to describe all the observed fabric patterns: for random c-axis distributions the diagonal entries of a⁽²⁾ are a⁽²⁾ ₁₁ = a⁽²⁾ ₂₂ = a⁽²⁾ ₃₃ = 1 /3, for a singlemaximum fabric with its maximum in the third direction, , and for a girdle-type fabric in the plane . When the material symmetry axes are those of the general reference frame, as for the three particular previous fabrics, the nondiagonal entries of a⁽²⁾ are nulls.

The behaviour of the polycrystal is described by the general orthotropic linear flow law (GOLF; Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and CastelnauGillet-Chaulet and others, 2005). In this law, ice is assumed to behave as a linearly viscous orthotropic material. In the present paper, following Reference Pettit, Thorsteinsson, Jacobson and WaddingtonPettit and others (2007) and Reference Martín, Hindmarsh and NavarroMartín and others (2009), the GOLF law is extended to a nonlinear form by adding an invariant in the anisotropic linear law. The simple choice is either to add the second invariant of the strain rate (Reference Martín, Hindmarsh and NavarroMartín and others, 2009) or the second invariant of the deviatoric stress (Reference Pettit, Thorsteinsson, Jacobson and WaddingtonPettit and others, 2007). No theoretical or experimental results are presently available that allow us to discard either of these two solutions, and other solutions based on anisotropic invariants of the deviatoric stress and/or the strain rate are also possible. Here the nonlinearity of the law is introduced through the second invariant of the deviatoric stress. With this definition, the anisotropy factors of the polycrystalline law for a given stress are identical in the linear and nonlinear cases. In other words, for a given fabric and a given state of stress, the corresponding strain rate relative to the isotropic response is the same for the linear and nonlinear cases. Using the strain-rate invariant, in the same way as Reference Martín, Hindmarsh and NavarroMartin and others (2009), leads to different anisotropy factors (as defined here) in the linear and nonlinear cases. Therefore, the proposed expression of the nonlinear GOLF law is

where A is the temperature-dependent Glen’s law parameter for isotropic ice. The six dimensionless anisotropy viscosities η _r (a⁽²⁾) and η_r+3(a ⁽²⁾) (r = 1,2, 3) are functions of eigenvalues of the second-order orientation tensor, a ⁽²⁾, which represent a measure of the anisotropy strength. The three structure tensors, Mr, the dyadic products of the three eigenvectors of a⁽²⁾, which then represent the material symmetry axes. In the method proposed by Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and RuokolainenGillet-Chaulet and others (2006), the six dimensionless viscosities, η _r (a⁽²⁾), are tabulated as a function of the fabric strength (i.e. the ai²⁾) using a micro/macro model. When the ice is isotropic, ηp_r = 0 and η_r+3 = 1 (r = 1,2, 3), and Equation (4) reduces to Glen’s flow law (Equation (1)).

Following Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and RuokolainenGillet-Chaulet and others (2006), the six dimensionless viscosities, η _r (a⁽²⁾), are tabulated using the viscoplastic self-consistent (VPSC) model (Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others, 1996, 1998). The two crystal parameters in the VPSC model used to tabulate the GOLF law were chosen to reproduce the experimentally observed polycrystal anisotropy. Following Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and RuokolainenGillet-Chaulet and others (2006), we use the shear strain rates ratio for a polycrystal with a single-maximum fabric and an isotropic polycrystal both experiencing the same shear stress. This anisotropy factor in shear is hereafter denoted k _s and, according to the experimental results of Reference Pimienta, Duval and LipenkovPimienta and others (1987), k _s ≈ 10. In other words, the VPSC parameters are chosen so that under simple shear a polycrystal with a single-maximum fabric is k _s times easier to deform than the corresponding isotropic polycrystal. The experimental results of Reference Pimienta, Duval and LipenkovPimienta and others (1987) indicate that an isotropic polycrystal is much easier to deform than a single-maximum fabric polycrystal experiencing the same uniaxial compressional stress. These experiments allow the definition of a second anisotropy factor for uniaxial compressional stress, denoted k _c. In what follows, the anisotropic polycrystal response deduced from the GOLF law (Equation (4)) is such that k _s = 10 and k _c = 0.4. The value k _c = 0.4 is in accordance with the experimental results of Reference Pimienta, Duval and LipenkovPimienta and others (1987) obtained from uniaxial compression tests on polycrystals with a single-maximum fabric. The anisotropy factors k _s and k _c are independent of Glen’s flow-law exponent, n, with the adopted nonlinear formulation.

Assuming that recrystallization processes do not occur and that the ice fabric is induced solely by deformation, the evolution of the second-order orientation tensor, a⁽²⁾, can be written as

where W is the spin tensor, defined as the antisymmetric part of the velocity gradient. The tensor C is defined as

The interaction parameter,α, controls the relative weighting of the strain rate and the deviatoric stress in the fabric- evolution equation (Equation (5)). When α = 0 the fabric evolution is solely controlled by the state of the strain rate, whereas when α = 1 the fabric evolves solely under the influence of the deviatoric stress. In between, as for the VPSC, both the strain rate and deviatoric stress contribute to the fabric evolution. In what follows, the interaction parameter is set to α = 0.06, in accordance with the crystal anisotropy and the VPSC model used to derive the polycrystal behaviour (Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and RuokolainenGillet-Chaulet and others, 2006). In Equation (5), the fourth-order orientation tensor is evaluated assuming a closure approximation giving a⁽⁴⁾ as a tensorial function of a⁽²⁾ (Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and RuokolainenGillet-Chaulet and others, 2006).

Flow equations

As shown in Figure 1, the domain geometry represents a two-dimensional flowline ice sheet and ice shelf. This geometry is the steady-state solution obtained with the same model assuming an isotropic rheology for step 12 of the MISMIP 3a experiment (Marine Ice Sheet Model Intercomparison Project; Reference Durand and HondohDurand and others, 2009b). Here we summarize the set-up of this experiment. The bedrock topography, in metres, is given as

where L _g is the length of the grounded part of the ice sheet. Using this bedrock and a uniform and constant surface accumulation of 0.3 ma^-1, the purpose of the MISMIP experiments was to compute various steady-state geometry configurations by varying Glen’s flow-law parameter, A. In this paper, only diagnostic simulations are performed and the upper and lower free surface, z = z _s(x) and z = z _b(x), respectively, are fixed and equal to the steady solutions obtained at step 12 (i.e. for A = 7.89 MPa^-3 a-1). The grounding-line position is also fixed and is arbitrarily taken to be x = 0. Diagnostic simulations are justified since the objective here is to infer the velocity enhancement due to ice anisotropy by comparing isotropic and anisotropic flows for the same geometry.

Fig. 1. Geometry of the ice sheet, notation and boundary conditions for the various experiments.

Ice-flow variables, velocity and isotropic pressure, are computed by solving the following field and boundary equations:

Equations (8) and (9) express mass and momentum conservation, respectively. Here u is the velocity vector, σ = τ −pI is the Cauchy stress tensor, p is the isotropic pressure, p _i and p _w are the density of ice and sea water, respectively, l _w is the sea level, n and t are the normal and tangential unit vectors to the free surfaces and g = (0, 0, -g) is the gravity vector. The upper surface is stress-free (Equation (10)) and for the lower surface the boundary conditions depend on whether the ice is grounded or floating. Where ice is grounded, the normal velocity is null and a Weertman-type friction law applies (Equation (11)), with C the friction parameter, whereas for the ice shelf bottom both normal velocity and tangential friction are null (Equation (12)). The left boundary condition (Equation (13)) expresses a symmetry axis in x = -L _g, and the right boundary condition (Equation (14)) is the end of the ice-shelf front, x = L _s. In the case of evolving fabric, the Stokes equations are coupled using the constitutive relation (Equation (4)) to the orientation tensor evolution equations (Equation (5)), as explained by Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and RuokolainenGillet-Chaulet and others (2006). The values of the various parameters in the equations are listed in Table 1. This set of equations is solved using the finite-element code Elmer/Ice; details on the numerics are given by Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and RuokolainenGillet-Chaulet and others (2006).

Definition of the enhancement factors

The objective of these simulations is to infer the value of the enhancement factors for the grounded part, E _SIA, and the ice-shelf part, E _SSA, by comparing anisotropic and isotropic results. It is easy to define, for each depth, values of the enhancement factors as a function of the local anisotropy and a given state of stress, but the objective here is to quantify the vertically integrated effect of anisotropy. Indeed, some models are vertically integrated, and the enhancement factors must then be determined from the global anisotropy effect from the surface to the bedrock.

For the grounded part, the dominant stress is the horizontal shear stress. From our full-Stokes anisotropic model, the enhancement in shear should be computed by comparing the shear strain rate, , for the anisotropic and isotropic cases. Because the shear strain rate varies with depth, the comparison has to be done on the resulting horizontal velocity, which is the horizontal surface velocity, u _s, minus the basal velocity, u _b (where basal sliding exists). Therefore, the enhancement factor in shear is determined as

The enhancement factor, E _SIA, is thus only evaluated from the part of the horizontal velocity induced by shear deformation. With our definition of the polycrystal anisotropy, we know that E _SIA ≤ k _s, equality being reached in the particular case of a homogeneous single-maximum fabric over all the domain.

For the ice shelf, the dominant stress is horizontal tension. The enhancement factor for the ice shelf must then be computed by comparing the horizontal strain rate, , for anisotropic and isotropic cases. In the case of perfect basal sliding, is almost constant with depth, so the enhancement factor can be evaluated at any elevation in the ice shelf. However, the SSA is also used to determine the ice- stream basal velocity. In this case, the tension enhancement factor should be inferred from the bottom value of the horizontal strain rate, . Therefore, we define the tension enhancement factor by

From the definition of the polycrystal anisotropy, the tension enhancement factor, E _SSA, should be greater than the anisotropy factor in compression, k _c (E _SSA ≥ k _c). Here again, the equality is reached for a homogeneous single-maximum fabric over the whole ice shelf.

Enhancement factor for grounded ice

The experiments are listed in Table 2. For all anisotropic experiments, the corresponding isotropic experiment (not listed) is used as the reference to compute the enhancement factor (Equations (15) and (16)). For the grounded part, two experiments were performed for various basal conditions, one assuming no sliding at the base (exp. 1) and the second with basal sliding (exp. 2). For these two experiments, the fabric is given and its evolution is not computed. Due to the fixed geometry and the strong changes in the flow condition at the grounding line, it was not possible to achieve a steady- state fabric field. Nevertheless, as shown below, the fabric was computed on the ice shelf where we did not have data input, using a simplified geometry model.

As shown in Figure 2 (dashed curve), for these two experiments the applied (fixed) fabric is a single-maximum fabric which becomes more and more concentrated from the surface to the bedrock as a function of the reduced depth. This is a typical fabric profile under a dome where ice is submitted to a compression stress, but it applies also along a flank where, due to shear, the fabric is known to evolve to a single-maximum pattern near the bedrock (Reference Seddik, Greve, Placidi, Hamann and GagliardiniSeddik and others, 2008). In what follows, this imposed fabric profile is referred to as the single-maximum (SMAX) fabric profile.

Fig. 2. Vertical fabric profiles given as the evolution of the diagonal components of the second-order orientation tensor, a⁽²⁾, for the imposed single-maximum (SMAX) fabric (dashed curve) and the computed steady-state solution in the ice shelf at x = 320 km (solid curve) and x = 800 km (dotted curve).

Results for the grounded part are shown in Figure 3. As expected (Reference Mangeney and CalifanoMangeney and Califano, 1998; Reference Gagliardini, Meyssonnier, Hutter, Wang and BeerGagliardini and Meyssonnier, 1999; Reference Staroszczyk and MorlandStaroszczyk and Morland, 2000), the surface velocity, u _s, for the anisotropic fabric is greater than that of the isotropic fabric, which confirms that the anisotropy makes the ice softer in the grounded part, where shear stress dominates. The enhancement factor evolution along x is calculated using Equation (15) for both with basal sliding and no basal sliding. Even if the surface velocities in the case of basal sliding are much greater than those in the case of no sliding, for the same anisotropy, both give the same enhancement factor E _SIA = 4.5. The same values of E _SIA derived from two experiments with various basal conditions indicate that, here, basal sliding has almost no effect on the ice shear deformation. It should be noted that the enhancement factor is also uniform on all the grounded part. Reference Mangeney, Califano and CastelnauMangeney and others (1996) obtained an enhancement factor of 1.7 and Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (2002) found a value of 1.75. These small values are explained by the uniform-stress polycrystal model used in their approach, which had a limited anisotropy factor, k _s < 2.5. With k _s = 5, Reference Staroszczyk and MorlandStaroszczyk and Morland (2000) obtained an enhancement factor of 4. We recall that in our work Ks is set to its experimental value of 10.1.

Fig. 3. (a) Horizontal velocity at the surface of the grounded part of the ice sheet as a function of distance from the grounding line, for the isotropic fabric (dotted curves) and SMAX fabric profile (dashed curves), and with and without sliding (triangle and no symbol). curves, respectively). (b) Corresponding enhancement factor in shear evaluated from Equation (15) in the case of sliding (triangles) and no-sliding (no symbol).

Certainly because the simulations are performed using a fixed geometry, irregularities of the velocity field are observed in the vicinity of the grounding line (over a distance of a few ice thicknesses). Because the problem of the grounding-line dynamics is not the intended subject of this paper, results in the vicinity of the grounding line are not discussed. These irregularities of the velocity field also explain why it was not possible to compute the steady-state fabric field associated with this ice-sheet geometry. At the dome, because of the vanishing horizontal velocity, the enhancement factor, E _SIA, is not defined.

Enhancement factor for the ice shelf

As the ice shelf is dominated by longitudinal stretching and lateral shearing instead of horizontal shearing, application of the SMAX fabric profile is inappropriate. We performed two furher experiments, exp. 3 and exp. 4, to infer the influence of fabric type on the flow. For these experiments, the domain was restricted to the ice-shelf part in order to be able to compute the steady-state fabric field, as depicted in Figure 4. For this restricted geometry, all the flow equations are the same, except that the left-side boundary conditions (Equation (13)) are modified to impose an inflow flux given by u(0, z) = 50 ma^-1 and v (0, z) = 0. This imposed inflow velocity corresponds approximately to the mean horizontal velocity at the grounding line for exp. 1.

Fig. 4. Geometry and boundary conditions for the ice-shelf experiments.

In exp. 3 a constant girdle fabric (a⁽²⁾ = 0.01, a)²⁾ = a3²⁾ = 0.495), representing the expected fabric for a longitudinal uniaxial tension stress, is imposed on all of the ice-shelf domain. For exp. 4 the fabric field is computed as the steady- state solution of the coupling between the flow equations and the fabric-evolution equations, assuming a SMAX profile as a boundary condition for the inflow boundary, x = 0. We checked that, whatever the initial fabric in the ice shelf (isotropic or SMAX profile), the problem converges to the same steady solution. This test with various initial conditions allows us to check whether the steady state of the ice flow is achieved or not and validates the numerical method used to compute the fabric field. The steady-state fabric profiles for the ice shelf at x = 320 and 800 km are shown in Figure 2 (solid and dotted curves, respectively). The evolution with depth of the computed fabric is different from the SMAX one but, due to the plane strain hypothesis, the fabric does not evolve to a girdle fabric. Comparison of these fabric profiles taken at two different locations indicates that the vertical evolution of the fabric is fairly similar throughout the ice shelf.

Results for the ice-shelf part are presented in Figure 5. In the absence of horizontal shearing, the velocity is constant with depth and its value increases almost linearly with distance to the grounding line. The surface velocities for the isotropic case and the three anisotropic cases (SMAX, girdle and steady-state fabrics) are shown in Figure 5a. The surface velocity, u _s, for the isotropic case is greater than any one of the anisotropic cases, indicating that anisotropic ice is stiffer under the stress conditions prevailing in ice shelves. For each experiment, the corresponding enhancement- factor evolution along the flow direction, inferred from Equation (16), is shown in Figure 5b. As for the grounded part, for a given experiment, the value of the enhancement factor, E _SSA, is almost uniform over the whole ice shelf. Its value is 0.71 for the SMAX fabric and 0.58 for the girdle fabric. For the computed steady-state fabric, the mean value of Essa is ˜0.57.

Fig. 5. (a) Horizontal velocity at the surface of the ice shelf as a function of the distance from the grounding line, for the isotropic fabric (dotted curve), SMAX fabric profile (dashed curve), girdle fabric (dot–dashed curve) and computed steady-state fabric (solid curve). (b) Corresponding enhancement factor in tension evaluated from Equation (16).

We also tested imposing an isotropic fabric as a boundary condition for the inflow boundary x = 0 (instead of a SMAX profile) and varying the value of the inflow velocity (from u(0, z) = 50 up to 300ma^-1). All the results (not shown) indicate that the inferred enhancement factors, E _SSA, forthese cases range between 0.52 and 0.71, greater values being obtained for the SMAX fabric profile. For the flow regime prevailing in ice shelves, whatever the fabrics, the flow enhancements induced by ice anisotropy are very similar.

Effect of non-isothermal conditions

In these applications, enhancement factors were inferred using isothermal conditions, whereas ice-sheet models use a non-constant temperature field. For the grounded part, increasing temperature with depth leads to a greater value of the enhancement factor, E _SIA, because it increases the shear deformation near the bedrock where the fabric is the most concentrated. We performed simulations assuming a fixed temperature field function of the reduced depth as observed at the Greenland Icecore Project (GRIP) site (Fig. 2a in Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 2000). These results indicate that for a variable temperature the enhancement factor, Esia, is equal to 5.6. Finally, isothermal and non-isothermal simulations give a range for E _SIA of 4.5-5.6, the greater value obtained for a non-constant temperature certainly being the more realistic.

For the ice shelf, the steady-state fabric obtained with isothermal conditions was used to test the influence of variable temperature conditions, assuming the same temperature profile as for the grounded part. The value of the enhancement factor in tension for a variable temperature is E _ssa = 0.60, to be compared with the value of 0.57 obtained with the isothermal conditions. These results indicate that an evolving temperature with depth has almost no influence on the resulting enhancement factor. In tension, contrary to what happens in shear, the flow is controlled by the layer with the lesser fluidity, i.e. the upper layers where the temperature is low.

Which enhancement factor for ice streams?

In the previous subsections, we have shown that the calculated enhancement factors in shear, E _SIA, and in tension, E _ssa , are almost uniform over the grounded part and the ice- shelf part, respectively. Let us now consider the transition zone between these two parts, where both the longitudinal and shear stress have the same magnitude. For this purpose, exp. 2 (SMAX fabric profile and sliding) was repeated for different friction coefficients, ranging from 0.3C to 1.2C, where C is the value used in the previous simulations and is given in Table 1. As shown in Figure 6, the value of the enhancement factor in tension, E _SSA, decreases from 1, where sliding is negligible, to its ice-shelf value (E_SSA = 0.71) just downstream of the grounding line. The smaller the friction parameter is, the longer the stream, and the decrease of E _SSA follows the stream length. This variation of E _SSA characterizes the fact that the longitudinal stress becomes progressively greater than the shear stress as the ice moves down in the ice stream.

Fig. 6. Enhancement factor in tension, E _SSA, as a function of distance to grounding line, evaluated from Equation (16) for different friction parameters from faster flow, 0.3C (circles) and 0.5C (lozenges), to slower flow, C (triangles) and 1.2C (squares). The value of C is given in Table 1 and corresponds to that used in the previous experiments.

Therefore, the value of the enhancement factor used in the SSA, if it is being applied to an ice stream, should not be constant and equal to that of the ice shelf, but should decrease from 1 at the onset of the ice stream to its ice-shelf value at the grounding line.

Again, because of the fixed domain geometry, the diagnostic velocity field in the vicinity of the grounding line is not smooth, inducing large oscillations of the enhancement factor in tension.