1. Introduction

As initially isotropic ice formed at the free surface of a large polar ice sheet moves downwards, it undergoes varying stress and deformation conditions, in response to which its internal polycrystalline structure evolves, giving rise to the development of directional properties of the material. Such a process, known as induced anisotropy, in which oriented fabrics in ice are created due to the crystal lattice rotation in the upper part of an ice sheet, and are subsequently slightly modified as a result of rotation recrystallization (polygonization) occurring in the middle part of an ice sheet, leads to the progressive strengthening ofthe anisotropic properties of the material. Evidence of this is provided by ice cores drilled in Greenland and Antarctica, showing, at large depths, strong fabrics with the majority of crystal c axes clustered around the vertical (Reference GowGow and others, 1997; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997). Near the base of an ice sheet, however, the above regular pattern changes dramatically, and ice fabrics with very coarse and interlocking grains develop. In these near-bottom fabrics, crystal c axes are usually broadly scattered in an irregular manner around the vertical (Reference GowGow and others, 1997), or multi-maxima fabrics resembling those observed in temperate glaciers are formed (Reference DuvalDuval,1981),with macroscopic properties of such fabrics close to isotropy (Reference LliboutryLliboutry, 1993). The transition from strong single-maximum fabrics to open fabrics may be quite abrupt, for it usually takes place over depths not exceeding the lowest 100 m of the ice sheet, compared to a typical ice-sheet thickness of 3–4 km. Such a rapid change in ice properties is due to the mechanism of dynamic recrystallization (also referred to as migration or discontinuous recrystallization) caused by migration of grain boundaries between dislocation-free and deformed ice crystals. As a result, under certain conditions, such as near-melting temperatures, high strains, strain rates and stresses, new grains are nucleated at the expense of old ones that ultimately disappear (Reference Lliboutry and DuvalLliboutry and Duval, 1985; Reference De La Chapelle, Castelnau, Lipenkov and DuvalDe La Chapelle and others, 1998).

Very few attempts have been made so far to model the process of dynamic recrystallization. Reference Van der Veen and WhillansVan der Veen and Whillans (1994), in their multi-grain model, have described the onset of this process by relating it to the total macroscopic strain in the polycrystalline aggregate. Recently, Reference MorlandMorland (2002) has proposed a model in which the onset of dynamic recrystallization mechanism is described by means of a temperature-dependent critical lattice distortion parameter, equivalent to a condition on a stored mechanical energy of dislocations. In another approach, Reference Faria, Kremer and HutterFaria and others (2003) construct a general theory of recrystallization processes in polycrystalline ice by applying the fundamental principles of thermodynamics. Inthispaper, we follow the approach formulated by Reference Staroszczyk and MorlandStaroszczyk and Morland (2001), based on the assumption that the onset of the migration recrystallization is controlled by strain rates and temperature. In this theory, fabric weakening at increasing strain rates is captured by introducing a scalar fabric-strength factor depending continuously on a single, temperature-dependent, effective strain-rate invariant, and it is supposed that the full strength of fabric anisotropy is retained below a lower critical value of this invariant, whereas above a higher critical value the fabric becomes isotropic. In Reference Staroszczyk and MorlandStaroszczyk and Morland (2001) the development of anisotropic fabric due to lattice rotation and its subsequent weakening due to migration recrystallization, and the effect of both processes on the viscous properties of ice, were illustrated for simple flow configurations encountered in laboratory conditions. Here, the theory is applied to investigate the influence of migration recrystallization on the large-scale behaviour of a polar ice sheet. Tothis end, a plane, gravity-driven, steady flow of an ice sheet on a horizontal bed with no-slip basal conditions, with fixed geometry and prescribed non-uniform temperature field, is analyzed. The results of numerical simulations, performed by a finite-element method, are compared with those for flows with no recrystallization involved (Reference De La Chapelle, Castelnau, Lipenkov and DuvalStaroszczyk and Morland, 2000) to illustrate the effect of dynamic recrystallization on depth profiles of velocities and deformation rates.

2. Problem Formulation

We consider a plane-strain, steady flow of an ice sheet taking place in rectangular Cartesian coordinates Oxz, with the x axis lying on the horizontal flat bed, and the z axis directed upwards and passing through the ice divide. For simplicity, the sheet profile is assumed symmetric about the plane x = 0.

Let σ be the Cauchy stress tensor, with components σ _xx, σ_zz and σ_xz in the plane of flow, and v be the ice velocity vector, with the horizontal and vertical components u and w, respectively. The deviatoric stress σ' is defined in terms of σ and the mean pressure p by

where I is the unit tensor and tr denotes the trace operator. In the absence of inertia forces, the momentum-balance equations are the equilibrium relations

where p is the ice density and g is the gravitational acceleration. The mass-balance equation is here the incompressibility condition div v = 0 given by

It is supposed that the ice temperature is significantly below the melting point at the bed so that there is no melt and no basal sliding at the interface, and at the ice divide we have the symmetry conditions. Hence, the boundary conditions at z = 0 and x = 0 are expressed, respectively, by

The stationary ice-sheet surface z = h(x) is assumed to be traction-free, and is subject to a net ice flux q(x) (positive for accumulation). The zero traction condition, and the kinematic condition prescribing the surface mass flux, are then given by

where n _s is a unit outward vector normal to the free surface, and u_s and w_s are the free-surface horizontal and vertical velocities, respectively.

In a manner typical for ice-sheet modelling, we eliminate physical dimensions from Equations (2.2) and (2.3) by adopting characteristic magnitudes: h*, a typical ice depth, used as a length scale; v*, a typical accumulation rate, used as a velocity unit; and τ* = pgh*, used as a stress unit. In terms of these scaling parameters, we introduce dimensionless variables defined by

Further, we apply coordinate stretching in the horizontal direction by adopting a small aspect ratio ϵ = h* /L* ≪ 1, where L* is the characteristic lateral dimension of the glacier, and introduce differentially scaled variables

As a result, both and coordinates, and both and velocity components, become order unity, and the maximum spatial derivatives in both directions have the same status. With the normalizations (2.4) and (2.5), the momentum and mass-balance equations (2.2) and (2.3) become

The deviatoric stresses in Equations (2.6) are prescribed in terms of the velocity components by the viscous constitutive relation described in section 3, while the pressure, which is not determined by the deformation due to the ice incompressibility, is a workless constraint stress. Therefore, the problem reduces to solving the system of three equations with three unknown functions

3. Constitutive Relation

The mechanism of induced anisotropy of ice, involving the strengthening of ice fabric due to lattice rotation, and weakening of fabric due to migration recrystallization, is described by an orthotropic constitutive theory formulated by Reference Staroszczyk and MorlandStaroszczyk and Morland (2001), which is now summarized.

The adopted viscous flow law expresses the deviatoric stress σ' in terms of the current strain rate D, the left Cauchy–Green deformation tensor B, and three structure tensors M ^(s) (s = 1, 2, 3) defined by the current (evolving) directions of principal stretches. Let the physical particle reference coordinates be OXi (i = 1, 3), with the equivalence X₁₌X and X₃ = Z, and apply to them the scaling analogous to that described by Equations (2.4) and (2.5) for spatial coordinates x and z. Then, with similar notations for spatial coordinates and velocity components, that is x ₁ = x,x₃=z,v₁= u and v₃ = w, the deformation gradient F, and the normalized velocity gradient L = L =D* and the strain rate with D* = v*/h*, have the following components expressed in the stretched coordinates:

In terms of Equation (3.1), the components of Cauchy– Green deformation tensor B = FF_T are

Following Reference Staroszczyk and MorlandStaroszczyk and Morland (2001), the constitutive law is expressed in an additive form consisting of two parts: the isotropic response of ice with no fabric, HI, and the additional anisotropic response depending on the evolving fabric, H_A, with the latter vanishing in an initial isotropic state B = I. Hence, in terms of the scaled variables defined by Equations (2.4) and (2.5), the normalized deviatoric stress is prescribed by

where the structure tensors M^(s) , describing orthotopic symmetries in the material, are defined by the tensor products

with e^(s) being the unit vectors in the directions of principal stretches (determined by the eigenvectors of the current deformation tensor B), and e⁽²⁾ = (0, 1, 0) in the reference frame adopted here. The isotropic and anisotropic terms in Equations (3.5) are given by

and involve two fabric response functions and each depending on a single deformation invariant argument: b_s (s = 1, 2, 3), being the squared principal stretches equal to the eigenvalues of B, and K = tr B = b₁ + b₂ + b₃. The vanishing of the anisotropic part H_A in the initial isotropic state B = I provides the normalization condition

In Equations (3.7), is a dimensionless viscosity of order unity, defined by where μ₀ is the physical viscosity of isotropic ice, commonly assumed to be a function of temperature T and the strain-rate invariant tr D₂. The temperature dependence of μ₀ is usually described by an Arrhenius-type relation. We, however, use the representation proposed by Reference Smith and MorlandSmith and Morland (1981), expressed by

with T μ a dimensionless temperature defined by 273.15 K, which seems to describe better the ice behaviour at high homologous temperatures. In this study we assume that is independent of tr D ₂ , that is suppose linearly viscous (Newtonian) behaviour of ice, a simplification which is justified for low (<0.1 MPa) deviatoric stresses (Reference Lliboutry and DuvalLliboutry and Duval, 1985).

To account for the dynamic recrystallization process leading to the fabric weakening with increasing strain rate, we extend the law given by Equations (3.5) and (3.7) by introducing a fabric-strength scaling function r and applying it to the anisotropic term H_A, which replaces Equation (3.5)1 by

The scaling factor r is supposed to depend only on a single strain-rate invariant tr D₂ and temperature through an effective strain-rate invariant (Reference Staroszczyk and MorlandStaroszczyk and Morland, 2001). We further introduce a critical strain-rate invariant at the centre of the range over which the fabric-strength reduction occurs. This range is defined by means of the lower and upper critical values and = Ī(1 + δ), where δ, a relative half-span of the critical range of is a free parameter of the model. The function is chosen to leave the viscous response unchanged when I I_{c l}, to restore the isotropic viscous response when and to vary continuously and monotonically between these two levels. Accordingly, is defined by

Now, in view of the relations (3.1), (3.3) and (3.4), the constitutive law defined by (3.10) with (3.7) gives the normalized stresses expressed in terms of the velocities by

where

We note that in the initial undeformed state, when B = I, that is b₁ = b₂ = b₃ = 1 and K = 3, in view of Equation (3.8), all the coefficients a _i(i = 1,. . .,4) become zero, regardless of the value of the factor r. And obviously, they are also zero for r = 0, that is for the effective strain-rate invariant I exceeding its upper critical level Ī_cu, after the migration recrystallization process has been completed. In either case, the formulae (3.12) reduce to the isotropic viscous fluid relations.

The model response functions and can be determined by correlation with the experimental data in the way described in detail in Reference Staroszczyk and MorlandStaroszczyk and Morland (2001) and Reference Morland and StaroszczykMorland and Staroszczyk (2003). Typical behaviour of ice under uniaxial compression and simple shear, with dynamic recrystallization involved, is illustrated in Reference Staroszczyk and MorlandStaroszczyk and Morland (2001). Here we use the constitutive theory to see how recrystallization affects the overall flow of a large polar ice sheet.

4. Numerical Results

The problem given by Equations (2.6) and (3.12), defined in terms of the velocities and and the pressure is solved by applying the finite-element method. Triangular elements with six nodes are used, with the velocity field approximated by quadratic interpolation functions and the pressure field approximated by linear interpolation functions. In total, 6414 elements, with 29 376 degrees of freedom, are used to discretize the ice-sheet domain. To find a steady solution, the flow problem is treated as time-dependent, starting from the isotropic state of ice throughout the glacier, that is assuming that F = I everywhere at the start of computation. Then an iterative procedure is followed, in which, successively, the velocities and pressures are first calculated from Equations (2.6), and these are used to calculate the changes in the deformation gradient components by solving the equation

where t denotes time, the superposed dot denotes the material derivative, and the summation convention for repeated indices applies. The changed components of F are then used to update the fabric described in terms of the current deformation tensor B = FF ^T , determining thus new values of and and the whole iteration process is repeated until steady flow is reached.

The results presented below have been obtained for a simple free-surface geometry defined by a parabola (h/H) + (x=L)2 = 1, where H is the ice-sheet thickness at the divide x = 0 and L is the lateral extent of the glacier. The simulations have been carried out for the aspect ratio ϵ = H/L = 10–₂. The temperature depth variation is assumed to be that determined for the GRIP ice core in central Greenland (Reference Gundestrup, Dahl-Jensen, Johnsen and RossiGundestrup and others, 1993), for which the isotropic viscosity ratio μo ( z = H)/μ0 (z = 0), given by Equation (3.9), is about 33.4. The magnitudes of the critical effective strain-rate invariant defining the onset of the recrystallization, have been chosen in such a way that the maximum thicknesses (occurring at x/L ~ 0.6 for the adopted glacier geometry) of the recrystallized ice are equal, respectively, to about 0.025, 0.05 and 0.075 of the sheet thickness h(x). In subsequent plots, the results corresponding to the latter values are labelled (2), (3) and (4), respectively, and label (1) indicates the results without recrystallization involved. The other free parameter of the model, δ, which determines the rate of the recrystallization process, has been assumed (arbitrarily due to the lack of any empirical indications) to be equal to 0.2.

Figure 1 illustrates the variation of the horizontal and vertical velocities, and at the glacier free surface for different levels of the critical strain-rate invariant I_c (solid lines). The results have been obtained for the ice enhancement factors: E_a = 1 for compression, and E_s = 5 for shear. For comparison, the results for E_a = 1/3 and E_s = 5 (dashed lines) and E_a = E_s = 1 (isotropic ice, dashed-dotted lines) are also plotted. It turns out that the global velocity field, for flows with and without dynamic recrystallization involved, is little sensitive to the value of E_a, and is very sensitive to the value of E_s, the shear enhancement factor. We note that the significant influence of recrystallization on the free-surface velocities is, essentially, confined to the part of the ice sheet where the thin bottom layer of recrystallized ice occurs. For instance, the presence of a layer of such ice of a relative maximum thickness 5% reduces the maximum free-surface horizontal velocity by about 18%; in the case of the vertical velocities, this effect is even more pronounced. In the near-divide or near-margin regions, where dynamic recrystallization is absent due to small shear strain rates there, the velocities are practically unaffected by the process.

Fig. 1. Free-surface horizontal and vertical velocities for the enhancement factors E, = 1 and E_s = 5 for different values of the critical strain-rate invariant I_c at which recrystallization occurs (solid lines). The dashed lines show the results for E, = 1/3 and E_s = 5, and the dashed-dotted lines show the results for isotropic ice, both with no recrystallization involved.

In Figure 2 we show, for different locations the depth profiles of both velocity components, and normalized by the respective free-surface values and and plotted against the normalized elevation z/h. Again, we note that the flow patterns are little changed by the recrystallization process both in the vicinity of the divide (plots for x/L = 0.2) and near the margin (plots for x/L = 0.8). In particular, this applies to the horizontal velocities whose profiles, even in the middle part of the ice sheet (x/L ~ 0.5), are found to be little sensitive to the recrystallization process (though the absolute velocity magnitudes differ from each other; recall Fig. 1). On the contrary, the normalized vertical velocity profiles in the middle part of the glacier are strongly affected by the presence of underlying isotropic ice near the base. Specifically, significant changes in the flow pattern are observed in the lower part of the ice sheet.

Fig. 2. Depth profiles of the normalized horizontal and vertical velocities at different locations x/L for different values of the critical strain-rate invariant I_c. The dotted lines show the results for isotropic ice with no recrystallization.

Finally, Figure 3 illustrates the variation of the shear strain rates for different flow regimes at two locations x/L in the middle part of the glacier, and shows as a function of the relative elevation z/h. We see that the magnitudes of the strain rates undergo abrupt changes in the thin transition layer separating the recrystallized “hard” ice from the overlying strongly anisotropic “soft” ice. We also note that with the increasing thickness of the bottom layer of “hard” ice, the shear strain rates in this layer decrease and approach the values for the isotropic ice flow with no recrystallization (dashed-dotted lines in the plots).

Fig. 3. Depth profiles of the normalized shear strain rate two locations x/L for different values of the critical strain-rate invariant Īc (solid lines). The dashed-dotted lines show the results for isotropic ice with no recrystallization.

5. Conclusions

The results of numerical simulations of a plane, steady flow of an ice sheet have shown that the process of dynamic recrystallization significantly affects the global flow pattern of the glacier. However, the results obtained should be treated with some at caution due to the restrictiveness of the constitutive model which, at this stage, describes the mechanism of dynamic recrystallization in terms of only temperature and strain rate. It seems that other factors are also important, of which the total deformation experienced by ice is likely to be the most crucial (personal communication from S. H. Faria, 2002). Therefore, future work should concentrate on appropriate extension of the constitutive theory, accompanied, if possible, by relevant experimental work.