Main features and restrictions of the model

The model presented here is of the same nature as the VPSC model, but additional assumptions are made in order to achieve a simpler formulation which could be used for ice-flow modelling purposes.

The fabrics which are likely to have a significant influence on the flow of ice caps are strong fabrics with caxes clustered around the vertical, which form under compression and/or simple shear, and that with coplanar caxes, which form under extension. In the first case, the enhancement factor on the strain rates is about 10, whereas in the latter configuration, when the ice is loaded along the symmetry axis of the fabric, it is less than 0.1 (Pimienta and others, 1987). Taking into account these particular fabrics would be a major improvement of the existing ice-sheet models. Intermediate fabrics with two maxima, as observed by Bouchez and Duval (1982), can be modelled by Azuma (1994) and Castelnau and others (in press) but evolve towards a single-maximum fabric with increasing accumulated strain. Multi-maxima fabrics, found in the “migration-recrystallization” zone near the bedrock, cannot be modelled without including the recrystallization processes. However, laboratory experiments (Duval, 1981) and temperate-glacier flow simulations (Meyssonnier, 1989) have shown that an isotropic power law is quite convenient to model the mechanical behaviour of this ice.

To capture the essential features of the in-situ observed fabrics, we propose to model polar ice as a transversely isotropic medium. Rather than considering the ice polycrystal as an assembly of a discrete number of individual grains, we replace the discontinuous distribution of the orientations of the caxes with a continuous orientation distribution function (ODF), which gives the density of grains with a given orientation.

To allow analytical calculations, each grain of the polycrystal is modelled as a transversely isotropic medium, the symmetry axis of which is the caxis of the crystal, with a weak resistance to shear parallel to the basal plane. This reduction of the actual behaviour of the ice crystal may be justified as follows:

Further simplification is achieved, in this first version of the model, by assuming a Newtonian behaviour for the grains. This does not render the model inappropriate, as some indications of a power-law exponent less than 2 can be found in the literature (Doake and Wolff, 1985; Lliboutry and Duval, 1985; Pimienta, 1987). Extension to a more general (non-linear) behaviour is possible but this is left for future work.

The most severe limitation of the model, in its current form, is that it is restricted to loading conditions which respect the symmetry of the fabric.

General framework

The starling point is Eshelby’s relation for an elastic (homogeneous) ellipsoidal inclusion in an elastic matrix or HEM with a free boundary at infinity. If this inclusion were removed from the HEM, then unloaded, it would take the shape , leaving a hole of shape ⁰ in the HEM. The transformation of ⁰ to * corresponds to a strain є* which is called the “stress-free strain” of . According to Mura (1987), the displacement u inside the inclusion, corresponding to an imposed arbitrary stress-free strain є*, is

(1)

where is the stiffness tensor of the HEM, x are Cartesian coordinates in the reference frame defined by the (geometrical) principal axes of the inclusion, and the G_ijkl components are integrals, over the surface of the inclusion, which depend only on the elastic moduli of the HEM and on the inclusion-shape ratio.

Following Gilormini and Vernusse (1992), we use the results obtained by Mura (1987) for an ellipsoidal inclusion with a geometrical symmetry axis parallel to the symmetry axis of a transversely isotropic HEM. The major advantage of this approach is that the expressions for G_ijkl can be integrated analytically. However, this is a further restriction on the model, and we will consider only the simplest case of spherical grains which are assumed to remain spherical during their deformation.

In the following, all macroscopic variables related to the HEM are overlined.

The displacement gradient inside the inclusion is derived from Equation (1) in the framework of infinitesimal strain and rotation. By applying the superposition principle, it follows that when the HEM is submitted to a prescribed strain ⋷ at infinity, corresponding to a stress , the strain ϵ in the inclusion is

(2)

where S is the Eshelby’s tensor of components,

(3)

In order to satisfy static equilibrium, the stress σ in the inclusion must be

(4)

where I denotes the identity tensor.

When considering a heterogeneous inclusion (inhomogeneity), the stress-free strain ϵ* is defined by Equations (2), (4) and the behaviour law of the inhomogeneity

(5)

After some algebra, the classical interaction formula is obtained as

(6)

The rotation tensor in the inhomogeneity, for a rotation ϖ prescribed at infinity on the HEM, is

(7)

where ω* is a consequence of the interaction between the inhomogeneity and the HEM. The components of ω* derive from Equation (1) as

(8)

where is the tensor of components

(9)

Using relations (7), (8) and (2), the rotation in the inclusion becomes

The aim of the following sections is to develop analytical forms for relations (6) and (10), then to apply these relations for ice by using the principle of correspondence.

Transverse isotropy

In the following, two Cartesian reference frames, and , are used (see Fig. 1):

Any (non-scalar) quantity, expressed as x in , is noted as x ⁱ when expressed in . The orientation of a grain is given by the two angles θ and ψ, defined in Figure 1.

The components a and a ^<i of the same vector with respect to and are related by a = R · a ⁱ , with

(11)

Using Voigt’s notation (the stress and strain tensors are expressed as vectors; σ = {σ₁, σ₂, σ₃, σ₂₃, σ₃₁, σ₁₂}, ϵ = {ϵ₁, ϵ₂, ϵ₃, 2ϵ₂₃, 2ϵ₃₁, 2ϵ₁₂}), the elastic compliance of the transversely isotropic HEM. expressed in , is noted

(12)

where . ᾱ is the ratio of the Young’s modulus Ē along the symmetry axis ( in Figure 1) to that in the plane of isotropy (, ӯ), ῡ and ῡ are the Poisson’s ratios in the planes parallel and perpendicular to the symmetry axis, respectively. is the ratio of the shear modulus in a plane which contains the symmetry axis to that in the plane of isotropy .

Fig. 1.

Macroscopic , and grain , Cartesian reference frames; angles λ and ψ define the c-axis orientation of a grain.

The elastic compliance of the inhomogeneity expressed in , takes a similar form in which Ḕ, ᾱ, , , ῡ₁, ῡ, are replaced with E, α, β, μ. ν₁, ν.

Incompressibility

Incompressibility is obtained by prescribing the values of ν, ῡ, ν₁ and ῡ as

so that

(14)

where p and denote the isotropic pressures in the inclusion and the HEM, then letting ῡ = ν approach 1/2. The analytical expressions for the Eshelby’s tensor components, at order zero in the small parameter 1 – 2ν, can be found in Gilormini and Vernusse (1992), who studied the case of an isotropic incompressible inhomogeneity in a transversely isotropic HEM. Taking the general symmetry relations into account, the non-zero components of S are

(15)

Additional relations which are useful to achieve results in closed form are

(16)

The expressions for the S_ijkl are given in the Appendix.

Interaction formulae for the ice polycrystal

The results obtained for elastic behaviours of the inhomogeneity and the HEM can be transposed to linearly viscous materials, of axial viscosities η and along the respective symmetry axes, by imposing incompressibility and replacing the strain ϵ with the strain rate ⋵, the rotation ω, with the rate of rotation , the stress σ with the deviatoric stress s, and the Young’s moduli E and Ē with 3η and 3, in relation (5)–(12). This is valid for infinitesimal strains and rotations from the current state and hence determines the instantaneous response from any deformed state.

Matrix given by relations (12) is now a viscosity matrix, and macroscopic parameter ᾱ is the ratio of the axial viscosity along the symmetry axis ( axis) of the ice polycrystal (HEM) to the axial viscosity in the plane of isotropy (i.e., the viscosity which would be measured in uniaxial compression along the or ӯ axis). Macroscopic parameter is the ratio of the shear viscosity in the or (ӯ−) plane to the shear viscosity in the plane of isotropy. This is summarized by

(17)

The same definitions apply to the grain (inclusion) parameters η, α and β.

The development of relation (6), using relations (12), expressed in frame , and relations (13)–(16), is lengthy but straightforward. Additional simplification was obtained by making α = 1 in the grain-behaviour form of relation 112); this does not conform to Pimienta’s (1987) results but preserves the essential feature of the grain behaviour (i.e., a weak resistance to shear parallel to the basal plane, obtained with β < 1). The corresponding form of relation (17) for the grain is then

(18)

By using the correspondence principle and separating the isotropic and deviatoric parts, the resulting interaction formulae, expressed in , give the strain rate ⋵ and the isotropic pressure p in a grain of c-axis orientation (θ, φ), as a function of the strain rate and of the isotropic pressure in the HEM. They are

(19)

where

(20)

and Γ₂₃. Γ31 are the expressions for the shear-strain rates and in the grain, expressed in by

(21)

with

(22)

Self-consistency

To describe the fabric, we use an QDF (Orientation Distribution Function) f(θ, φ), which gives the density of caxes with orientation (θ, φ), i.e. the relative number of caxes intersections with the unit sphere per unit area of the sphere. The relative number of grains whose orientations lie in the interval (θ, θ + dθ; φ, φ + dφ) is

(23)

Assuming that there is no correlation between the volume of a grain and its orientation, dN _r is also a volume fraction.

By definition,

(24)

For a transversely isotropic HEM, f does not depend on φ, i.e., f(θ, φ) = f(θ). The weighted average of a quantity x(θ, φ) is then defined as

(25)

For an isotropic HEM: f(θ, φ) = 1.

The self-consistency of the model is achieved by setting the averages of the isotropic pressure and of the strain rates equal to their corresponding macroscopic quantities, that is

(26)

Using expressions (19) for p and ⋵ _ij , with the complementary relations (15), (16) and (20)–(22), relations (25) and (26) result in the following non-linear set of three equations for the unknowns , ᾱ, :

(27)

The solution of Equation (27) can only be found numerically, except for the simple isotropic fabric. In this case ᾱ = = 1 and the following relations hold:

(28)

Then system (27) reduces to the equation (3)(–2β = 0, which has only one positive root given by

(29)

Thus, for small values of β, the macroscopic viscosity of the HEM is controlled In the main viscosity along the c-axis direction.

Fabric evolution

According to Equation (11), the components c and c ⁱ of the c-axis vector of a grain are linked by

(30)

Under the velocity gradient L in the grain, c transforms into c + dc such that

(31)

Differentiating Equation (30) with respect to time leads to

(32)

It follows that

(33)

By introducing the rate of rotation tensor , defined by L = ⋵ + , and using the transformation formula ⋵ = R η ⁱ η R ⁻¹, expression (11) for R and c ⁱ = ͻ0, 0, 1ͽ, the change in the orientation of a grain follows from Equation (33) as

(34)

Assuming that the basal planes of a grain remain parallel to each other during the deformation, the component of the velocity in a grain along the z axis (caxis) is a function of z only, so that

(35)

The components of the rate-of-rotation tensor in the global reference frame are given by Equation (10) (written for rates). The Ω _ijkl components which appear in this equation, not given by Gilormini and Vernusse (1992), need to be calculated from relation (9) and Mura’s (1987) expressions for the G_ijkl . For a transversely isotropic HEM, the only non-zero components are

(36)

From Equation (10) and (36), the non-zero components of in are

(37)

Using Equation (35), (37) and the interaction formulae (19)–(22) in relations (34), completely determines and .

The ODF evolution derives from the continuity Equation (24), by requiring that the net flux of grains which enters the interval (θ,θ + dθ; φ, φ + dφ) during a unit time is equal to the increase in the number of grains in this time, i.e.

(38)

This equation, associated with Equation (34), describes the evolution of the fabric.

Calculation procedure

The equations to be solved are the self-consistent system (27), which gives the mechanical behaviour of the HEM , and Equation (38) for the evolution of the ODF, using Equation (34).

Under uniaxial compression or tension along the HEM symmetry axis, the velocity gradient prescribed at infinity, on the HEM boundary, is such that

(39)

Under these conditions, Equation (34) simplifies to

(40)

As Γ₃₁ does not depend on φ, these equations confirm that the symmetry is maintained when compression, or tension, is exerted along the symmetry axis of the HEM.

Even though the present model is only valid for small strain and rotation, in this uniaxial geometry with no macroscopic rotation, a finite uniaxial strain can be reached by summing the small strain increments which arise in small time increments. The results can only be achieved by numerical computation, as the self-consistent set of Equation (27) is non-linear and the ODF is not explicated analytically. Starting from an isotropic configuration (f(θ) = 1, = 1), system (27) was solved by Newton’s method which was found to be always convergent. At the end of each time step, the macroscopic parameters , ᾱ s, then the components of Eshelby’s tensor and the ODF were updated. Solution of Equation (38) was achieved by discretizing f(θ), then balancing the net flux of grains entering each interval (θ, θ + dθ) during the time step dt with the increase in f(θ) by using an upwind scheme which takes into account the sign of at the ends of each θ interval.

Results

figure 2a, and b shows the evolution of the ODF, under compression and tension, respectively, as a function of the equivalent strain . This equivalent strain has no physical meaning, as the model does not hold for finite strain, but it is a convenient reference measure for intercomparison between different models. As expected, the caxes tend to gather around the compression axis, or to fall in the plane perpendicular to the direction of tension, respectively, with increasing equivalent strain. As the relative number of grains in the interval (θ,θ + dθ) is f(θ) sin θ dθ (see Equation (23)), the usual representation of these fabrics by a Schmidt’s diagram would be two girdles, the radii of which decrease or increase respectively, with increasing equivalent strain. These curves were computed with β = 0.04 which is the solution of , with given by Equation (29). With such a value of β, according to Equation (18), a polycrystal with all its caxes parallel, undergoing simple shear parallel to the basal plane of the grains, would deform ten times faster than an isotropic polycrystal, in accordance with Pimienta and others’ (1987) estimate (all other directional viscosities being equal to that of a grain along its symmetry axis).

Fig. 2.

Orientation distribution function f(θ) computed for different values of , with grain parameter β = 0.04. under (a) uniaxial compression, (b) uniaxial tension.

Fig. 3.

Evolution of the macroscopic parameters , as a function of the equivalent strain ϵ_eq, for grain parameter β = 0.04, under (a) uniaxial compression, (b) uniaxial tension.

The evolution of the mechanical parameters , and , are shown on figure 3a and b. Under compression, in the extreme case of aligned caxes, the behaviour of the polycrystal should be that of the isolated grain (i.e., ). This final state seems difficult to reach:; it ϵ_eq = 1, the macroscopic parameters are found to be . Taking the viscosity of the isotropic polycrystal as a reference, the corresponding viscosities at ϵ_eq = 1 are, according to Equation (17): . Under tension, in the extreme case of coplanar caxes, the viscosity along the polycrystal symmetry axis tends towards the axial viscosity of the grain. At ϵ_dq = 1, the macroscopic parameters are , and the corresponding viscosities are: .

The influence of the grain parameter β on the evolution rate of the fabric, under uniaxial compression, is shown in figure 4a. The influence of β on the mechanical parameters of the HEM is shown in figure 4b; the stronger effect is on . Note, however, that at an equivalent strain of 1, the values of computed for β = 0.1 and β = 0.01 are only in a ratio of 3, so that this influence is not too great.

ODF Parameterization

All of these results were obtained by discretizing the ODF on 90 intervals between 0 and π/2. By doing so, there is little advantage in using an ODF formulation rather than discretizing the polycrystal in a finite number of grains, as has been done in Castelnau and others (in press) model. To remedy this situation, an attempt was made to parameterize the ODF with a small number of parameters. By analogy with the results of Giessen and Houtte (1992), who derived analytical formulations for the ODF in a model based on Taylor’s assumption, it was found that, both in compression and tension, the present ODF can be very accurately fitted with a function of the form

(41)

A theoretical demonstration of this, in the framework of the present model, is not straightforward, as the coefficients in Equation (40) for depend on ᾱ and which are implicit functions of f(θ).

The ODF’s obtained, with an evolving discretization of the ODF, or by using approximation (41) at the end of each time step to replace the actual ODF with f*, then discretizing f* at the beginning of the next time step (which mimics what should be done during a flow simulation), are practically identical, either under compression or tension, even at a large equivalent strain.

The generalization of such a parameterization, for loading conditions other than compression or tension, will retain the storage capacity required to solve a flow problem to a minimum.