1. Introduction

Generally (see, e.g., Reference HutterHutter, 1983), the creep of ice sheets and polar glaciers is modeled as the flow of an isotropic, heat-conducting, incompressible and structureless medium, for which there hold the following balance equations of mass, momenta, energy and entropy, respectively:

where the superposed dot represents the material time derivative, ν is the axial vector of the stress tensor and the superscribed T denotes the transpose of the respective tensor (further details of notation are explained at the end of this section). The usefulness of Equations (1–5) is actually twofold. On the one hand, Equations (1), (2) and (4) can be used to determine the basic thermodynamic fields of mass density ϱ(x, t) , translational velocity v(x, t) and specific internal energy e(x, t) (or alternatively the absolute temperature T(x, t), if a thermal equation of state relating both quantities is available) as functions of the position X and time instant t, provided that the external supplies of linear momentum g(x, t) and internal energy r(x, t) are prescribed, together with constitutive equations for the remaining fields, viz. the Cauchy stress t(x, t), and the heat flux density q(x, t). On the other hand, Equations (1–5) can be used instead to derive the explicit forms of such constitutive functions, which are constrained by the second law of thermodynamics through the entropic inequality which emerges from Equation (5), after elimination of the specific entropy production ζ (x, t). This inequality introduces two additional undetermined constitutive quantities, namely the specific entropy η (x, t) and its flux density ϕ(x, t), while s(x, t) is a prescribed external supply of entropy.

Hence, Equations (1–5) constitute a complete thermodynamic framework for the construction of constitutive theories, from which the celebrated laws of ice-sheet modeling, viz. Glen’s flow law and Fourier’s law of heat conduction, can be derived (Reference HutterHutter, 1983). Unfortunately, such an approach does not support adequate descriptions of evolving anisotropy and recrystallization phenomena, simply because these processes are essentially of microstructural nature. Notwithstanding, due to the increasing demand on improvements of Glen’s law, anisotropic enhancements have been proposed during the last 15 years, without much care about thermodynamical first principles. This has given rise to the critical question of whether most of the anisotropic ice-sheet models presented so far are indeed portraitures of physical reality or just ingenious casually coincident mimes of observed data (Reference Faria, Hutter, Kim and JungFaria and Hutter, 2001).

Owing to this, the aim of this work is to oppose the idiosyncrasy of ad hoc methods, by presenting a completely self-consistent thermodynamic approach for the construction of microstructurally based constitutive theories for polar ice. The texture (also called fabric¹) of the material is described by a continuous distribution of c axes, thus allowing the characterization of the most diverse types of anisotropy while recrystallization processes, like grain growth, polygonization and nucleation, are mainly driven by the density of dislocations. An essential feature of this theory, as in any continuum approach to ice-sheet dynamics, is the assumption that every ice-sheet “particle” encloses a huge number of grains within a so-called elementary volume, located at X (its center of mass). Despite the amount of grains contained in it, such elementary volume is regarded as so minute in comparison with the ice-sheet dimensions as to be treated as a “point” of the continuum.

Concerning the notation, scalar-, vector- and tensor-valued fields are represented by Greek and Latin letters. In particular, scalars are written in lightface, vectors in slanted or italic boldface, and tensors in upright boldface. Components of vectors and tensors are expressed in Cartesian indicial notation using italic Latin indices (i = 1, 2, 3), together with the usual summation convention of repeated indices. For instance, some frequently used definitions are (a, A , a, A, denote arbitrary vectors and tensors, respectively)

where e_ijk and ν_i represent the components of the permutation tensor and of the axial vector of stress.

2. Continuous Diversity Concept

The high anisotropy of the visco-plastic response of ice crystals is notorious: mechanical tests indicate that, within the range of temperatures found in glaciers and ice sheets, ice single crystals manifest transversal isotropy about their crystallographic axis of symmetry, the so-called c axis (McConnel, 1891; Reference SteinemannSteinemann, 1958; Reference GlenGlen, 1975; Reference Budd and JackaBudd and Jacka, 1989). Owing to this, it can be said that, in a polycrystal, ice grains with different c-axis orientations might be regarded, from a thermomechanical point of view, as “distinct materials”, even though they are made of the same chemical substance. Following this conjecture, we conclude that the determination of the mechanical response of a large polycrystalline ice mass via homogenization should resemble somewhat the derivation of the properties of a mixture from the respective attributes of its components. Nevertheless, unusually for a chemical mixture, there exists in polycrystalline ice no definite number of “constituents”, but rather something like a continuous diversity of possible grain types or, more precisely, lattice orientations.

Clearly, in order to describe the physics of what could be called a“mixture of orientations”, some generalization of the classical theory of mixtures is required. In fact, such a generalized mixture theory has been recently addressed, though in a broader context, under the name mixtures with continuous diversity (Reference FariaFaria, 2001). Succinctly, a mixture with continuous diversity can be regarded as a continuum composed of an infinite number of interacting constituents whose distinctive properties vary smoothly from one to another. Mathematically, such a generalization of classical mixture theory can be readily understood by reckoning, for instance, the constituent field of mass density in an ordinary chemical mixture of N components, viz. ϱ^α (x , t) , with α = 1; 2; . . . ; N. To derive the respective field for a mixture with continuous diversity we simply allow the species label α to become a real variable (instead of an integer number) defined in an interval like [α_min, α_max] = A C IR. The extrema α_min and α_max are so prescribed as to enclose in the interval A all species which can possibly be found in the medium. In this theory, the species label α should have some physical meaning, attributed according to the relevant distinctive property of the species. Moreover, it is straightforward to demonstrate that multiple labels can be introduced (Reference FariaFaria, 2001), usually associated with the components of vector- and tensor-valued distinctive properties.

In particular, for the case of polycrystalline ice we find that at least two of those labels are necessary, since any c-axis orientation is uniquely defined in the three-dimensional Euclidean space by its corresponding polar and azimuthal angles, θ and φ, respectively. I n principle, we could simply set α₁ = θ and α₂ = φ, but calculations using this choice of labels are rather cumbersome. Instead of this, we introduce for convenience a (dimensionless) unit vector n , called orientation vector, whose spherical coordinates are just 0, φ and r = |n| = 1. It is easy to verify that n is nothing else than the radius vector of the closed, unit spherical surface S₂ embedded in the usual three-dimensional space (see Fig. 1). Owing to this, S₂ is referred to here as the orientation space.

Fig. 1. Relationship between the orientation η and its characteristic angles, θ and φ.

It should be noticed that the new vector-valued species label, n , has in fact achieved the status of a new field variable, in addition to X and t, in such a way that the “constituent” mass density field in a mixture with continuous diversity reads ϱ*(x, t, n), where the “*” indicates that the respective field is orientation-dependent, i.e. it is defined, for any instant t, in a five-dimensional space whose coordinates are (x, n) . Of course, the same procedure can be extended to any other material quantity, thus enabling the definition of the orientation-dependent fields of translational velocity V*(x, t, n), Cauchy stress t*(x, t, n), specific internal energy e*(x, t, n), heat flux density q*(x, t, n), specific entropy η* (x, t , n and its flux density y ϕ*(x , t, n), etc. Furthermore, there is a remarkable characteristic off mixtures with continuous diversity which is not typically found in chemical mixtures; it is called familiarity. Generally, two species are said to be familiar if their distinctive properties (and consequently their physical behavior) are alike, though not identical. Clearly, the concept of familiarity is a direct counterpart of the notion of closeness in the Euclidean space and it is necessary in order to put the orientation n and the position X on the same footing, in a kind of “generalized spatial (Eulerian) description”. From this reasoning, it has been shown by Reference Faria and HutterFaria and Hutter (2002) that the following balance equations of mass, linear and angular momenta, energy and entropy derive from fundamental principles

where we have introduced the notation (cf. also section 1; the symbol a* represents any arbitrary tensor field)

It is worth noticing that, in spherical coordinates, the components of the orientational gradient operator (grad) are identical to those of the spatial gradient operator (grad) in two dimensions, i.e. subjected to the constraint r = 1.

From Equations (12) and (13) we immediately recognize that ice crystallites are treated here as microstructured media, with m* and σ* = Is* denoting respectively the fields of couple stress and intrinsic angular momentum per unit mass (i.e. specific c-axis spin), where S* denotes the angular velocity of spin (or c-axis spin velocity) and I is a material constant, called spin inertia. Moreover, in Equations (10–14) the fields g* , C*, r* and s* denote respectively the specific external supplies of momenta, energy and entropy, while and ς* are the specific production terms of mass, momenta, energy and entropy, due to inter-species interactions and dissipative processes. In particular, it is convenient to split the production of energy ζ* into diffusive (h*) and non-diffusive (ε*) parts, viz. ζ* = h* + ε*. Further, ξ* and φ are the conductive flux densities of momenta, energy and entropy in the orientation space, due to interactions between familiar species (low-angle interactions). Finally, 77* represents the specific internal energy production due to the work of stresses, viz.

while u*, called transition rate, measures the rate at which matter alters its distinctive properties. For the particular case of polycrystalline ice it has the connotation of a kind of orientational “velocity”, describing the rate at which c axes continuously change their orientations (e.g. by lattice rotation or polygonization). Consequently the transition rate can be related to the c-axis spin velocity through the relation u* = s* x n.

Equations (10–14) constitute the fundamental equations for polycrystalline ice modeled as a mixture with continuous diversity. Notwithstanding, in the present form they are neither sufficient nor adequate to describe fabric (texture) evolution and recrystallization in ice sheets. The justification for this and the derivation of the appropriate balance equations are addressed in the next section.

3. Balance Equations for Polar Ice

There exists a marked similarity between Equations (10–14) and the balance equations for chemical mixtures, which arises from the fact that grains in a polycrystal, just like the constituents of a chemically reacting mixture, are strongly interacting open systems. Nevertheless, the contrasting features of both systems are also evident. To understand this, we start recalling that, in the present description, a single ice-sheet particle located at X encompasses in fact a huge number of grains (see section 1). In this sense, the point X is macroscopic. On the other hand, we observe that there exists just one translational velocity assigned to a given n in this particular point, which implies that V* actually describes the (instantaneous) average translational motion of the crystallites whose c axes are parallel to n , at position X . Yet experience tells us that diffusive motion of grains is not a typical feature found in ice sheets. In other words, crystallites with a particular lattice orientation do not pass through the material, but rather they move with the bulk, even though they do not experience the same deformation as the aggregate. This allows us to vanish the translational diffusion by setting

We emphasize that Equation (19) has no relation at all with constraints imposed on the crystalline deformation (like the so-called Voigt–Taylor and Sachs–Reuss upper and lower bounds), since the scale considered here is indeed macroscopic. This means that strain incompatibilities on the grain level have already been smeared out by homogenization, in such a way that we can consider simply the motion of the ice-sheet particles, just as any other macroscopic continuum theory of ice sheets would do. The reason for this is that, in contrast to most models of anisotropic ice, the description of kinematics and microstructure evolution in the present formalism is not based on an averaging of the mechanics of individual crystallites, but rather on the general diffusion theory of structured continua. As a consequence, the results of any microscopically based homogenization method (e.g. upper bound or visco-plastic self-consistent approaches; see, e.g., Reference Castelnau, Thorsteinsson, Kipfstuhl, Duval and CanovaCastelnau and others, 1996) should be compatible, in a large scale, with the relations presented in this work. Whereas the assumption (19) could eventually be controversial for a small sample of ice in the laboratory, where a material particle contains merely a few grains, for large ice masses like glaciers and ice sheets this seems to be incontestably the most natural choice.

Akin conclusions can be established for the rotational motion of the lattice: since the spin inertia of the c axes is negligible (I -+ 0, i.e. crystallographic axes do not continue rotating after stresses are suspended) and the spin velocity 3* is limited, it follows that σ* = Is* → 0. In contrast to Equation (19), however, it must be emphasized that the spin velocity g*(x , t, n) generally depends on the orientation n . Finally, we assume for simplicity that g* = g, C* = 0, r* = r and s* = s, reflecting the fact that external body supplies of heat (through solar radiation) and momenta (through gravity)do not produce couples and are effectively independent of the distribution of c axes (Reference HobbsHobbs, 1974). Hence, from the above arguments we obtain the reduced balance equations for polar ice

where we made use of Equations (1), (9) and (19). Moreover,

Further simplifications can still be derived through a suitable constitutive theory (e.g. the vanishing of couple stresses, as will be discussed later). Nevertheless, before the constitutive analysis is approached, it must be recognized that there is still something lacking in the set of Equations (20–24). Indeed, while Equation (24) expresses thermodynamical irreversibility in a way to be explained later, Equations (20– 23) determine the basic fields ϱ*, V, S* and e* (or alternatively, temperature) associated with mass exchange, motion, c-axis rotation and thermal processes, respectively. Yet, although recrystallization phenomena are occurring in all these equations (through the mass production I * and other interaction terms), experience shows us that no combination of the basic fields mentioned above suffices to describe recrystallization. Hence, we need an additional balance equation for a new quantity, related to the recrystallization energy. We know that such an energy is stored from deformation. And we know that the deformation energy is mainly stored in linear lattice defects: the dislocations. Therefore, what we need is a balance equation for the dislocation density.

4. Dislocation Density Evolution

In order to derive the balance equation of dislocation from fundamental principles, we must temporarily move to its minute world, the micro crystalline level. A naive sketch of this is given in Figure 2a. Following Reference KrönerKröner (2001), we introduce the Burgers vector 6, characteristic of the dislocation, and a (dimensionless) unit vector t, tangential to the dislocation line. With these vectors we build the diadic t ⊗ b, called dislocation tensor. As discussed by Reference Kröner, Kanninen, Adler, Rosenfield and JafeeKröner (1970, Reference Kröner2001), this tensor completely describes the microscopic state of a dislocation and is closely related to the well-known (macroscopic) dislocation density tensor of classical continuum theories (see Reference Truesdell, Noll and S.Truesdell and Noll, 1965, and references therein). Therefore, we will use the dislocation tensor here to define, through homogenization, the fundamental macroscopic quantities which must appear in the dislocation balance equation.

Fig. 2. (a) Microscopic volume Vμ enclosing a basal dislocation and a grain boundary on its right. (b) Sketch of a volume of size between Vμ and V. The unshaded region near the boundary has a thickness of the order of the dislocation mean free path, λ.

Let V represent the spatial region occupied by an elementary volume of ice (i.e. an ice-sheet particle; see section 1) located at X, containing a huge number of dislocations. It is convenient to assume for a while a static situation in which c axes do not rotate and no mass flows across the boundaries of V. The Burgers vectors and the tangential unit vectors of the dislocations within V are respectively described by the distributions and where r* denotes a position within the elementary volume, i.e. on a microscopic scale. If L represents the set of those microscopic points which are instantaneously occupied by dislocations, then it can be easily recognized that

where, for simplicity, we have considered only one type of Burgers vector, whose magnitude b is constant. Hence, by dividing the volume V of V in a great number N of atomic cells of volume b3, we can introduce through homogenization the orientation-dependent Kröner autocorrelation function, viz. (Reference Kröner, Kanninen, Adler, Rosenfield and JafeeKröner, 2001)

where d3r denotes a volume element. In the relation above, we have used Equation (26) to convert the homogenization integral of Equation (27) into a sum over all atomic cells. The quantity n*(x, t, n) defines the instantaneous amount of such cells containing dislocation segments (for given elementary volume and lattice orientation). Now, if l ≈ b denotes the mean dislocation length within an atomic cell, we can derive, through multiplication of Equation (27) by b-4, the identities

which define the orientation-dependent dislocation density p*(x, t, n) as the total length of dislocations (within grains whose c axes are oriented to n) per unit volume. Hence, p* constitutes a new internal variable on the macroscopic scale, related to the microscopic distribution of dislocations through the homogenization integral (27).

From the master balance law of continuum theory (Reference Faria and HutterFaria and Hutter, 2002; see also Reference HutterHutter, 1983; Reference LiuLiu, 2002) we know that the time rate of any additive quantity must be balanced by external supplies, internal production, convective and conductive fluxes. For dislocations there are obviously no external supplies, and, due to the static situation temporarily assumed here, convective fluxes are absent. Dislocation production (by Frank–Read sources, grain boundary absorption and emission, etc.) is denoted by П* . To compute the conductive flux of dislocations, we adopt the same homogenization method employed above, but now in two steps.

First, we consider the case of a microscopic region V μ like that shown in Figure 2a, whose volume, Vμ = Nμb3, is some fraction of the mean grain-size (the subscript , 1 indicates that the value of the respective quantity depends on the size and geometry of the chosen region Vμ). In this case, the density of dislocations is denoted by and, owing to Equation (28), it can be easily derived from Equation (27) by simply substituting Vμ and Nμ for V and A , respectively. Further, let us suppose that the (orientation-dependent) average propagation velocity of such dislocations is When the dimensions of the region Vμ are comparable with the dislocation mean free path A, most dislocations can move through Vμ without serious obstructions and the flux density of dislocations is given simply by This is a conductive flux, because no mass transport is involved with it.

Yet, when the volume of the microscopic region Vμ is increased up to the size V of a macroscopic elementary volume V, we conclude from Equation (27) and Figure 2b that the conductive flux of dislocations is drastically reduced by a factor proportional to λ/V₁/₃. This is because dislocations separated from the boundary of V by distances larger than A are usually annihilated or immobilized before having the chance to reach that boundary. Hence, for sufficiently large elementary volumes, as in the case of a glacier or ice sheet, only the dislocations situated in an infinitesimally thin layer of V (represented by the unshaded region in Fig. 2b) contribute to the conductive dislocation flux density, which consequently becomes negligible. Therefore, we conclude that in the absence of matter transport and c-axis rotation, the time rate of the dislocation density is simply balanced by its internal production, i.e.

Finally, in a general dynamic situation in which the material in V moves and rotates arbitrarily, additional convective fluxes of dislocations must be accounted for, due to mass transport, p* V* , and c-axis rotation, p*U* . Consequently, from the above arguments and the master balance law already mentioned, there follows the generic form of the dislocation balance equation

which can be rewritten with the aid of Equations (1), (9) and (16) as

Equation (30) is the wanted balance equation of dislocation for polar ice. Admittedly, this equation is based upon a model (which can be improved, for example by distinguishing between different kinds of dislocations), but in spite of this, it is rigorous within this context.

5. ON Constitutive Relations

What we have presented so far is a continuum approach for the thermodynamics of polycrystalline polar ice involving five orientation-dependent basic fields of mass density ϱ* , translational velocity V , c-axis spin velocity S* , internal energy e* (or alternatively, temperature) and dislocation density p*, which can be determined by the balance equations (20–23) and (30) provided that constitutive representations for the quantities and are prescribed. Notwithstanding, such constitutive relations are not arbitrary, since they must satisfy the second law of thermodynamics, which expresses the inherent irreversibility of natural processes through the orientation-dependent entropy inequality, ς* ≥ δ* (Reference FariaFaria, 2001). Insertion of it into Equation (24) yields an entropic inequality with four additional constitutive quantities, namely η* , ϕ*, φ* and the so-called entropy production deviation, δ*(x, t, n). This last represents the part of the entropy production due to interspecies interactions, i.e. interactions between crystallites with different c-axis orientations, and can be written as δ*=ς* −ς* ^_P . Evidently, when the material at hand is a single crystal, no deviation exists. Conversely, the vanishing of 6* in a polycrystal indicates that grains with different lattice orientations behave independently, as isolated systems. In this case, ς* = ς* ^_P , which is the specific entropy production of an isolated, hypothetical ice single crystal under similar conditions.

Thermodynamically consistent constitutive equations for the quantities already mentioned have been derived through exploitation of the orientation-dependent entropy inequality by using standard methods, such as Liu’s Lagrange multipliers (Reference LiuLiu, 2002). Details of such an exploitation are still under analysis and will soon appear elsewhere, but a few of the relevant issues can already be advanced. For instance, it can be shown that the time rate of the dislocation density is an important constitutive variable, which should not be discarded because, among other reasons, it reflects the influence of recovery on the material response. Further, we have found that each of the production terms ε *, δ* and can be decomposed in two parts, one of them directly proportional to the mass production There fore, recrystallization processes involving intercrystalline mass exchanges (e.g. nucleation or grain growth) have inevitable effects upon the energy, entropy and dislocation density of the material. Finally, the most interesting outcome of the constitutive theory at this stage is that the couple stress, m*, must vanish for ice sheets. Intuitively, this can be readily justified by the fact that couple stresses are generally proportional to gradients of the spin velocity S* (cf. Equation (25)). Yet experience shows us that the spin velocities of c axes in polar ice are roughly of the same order as the strain rate (Reference Kamb, Heard, Borg, Carter and RaleighKamb, 1972; Reference Jacka, Jun and T.Jacka and Li , 2000). Consequently, spin velocity gradients have the magnitude of gradients of strain rate, which are usually neglected in constitutive theories for ice sheets.

The vanishing of couple stresses leads to a curious conclusion concerning the mechanics of polycrystalline ice: whereas the stresses on individual grains are generally non-symmetric (cf. Equation (22)), the polycrystalline stress, averaged over all c-axis orientations, is indeed symmetric. To demonstrate this, we recall the homogenization formulas for all fields occurring in Equations (20–23) and (30), which can be obtained as a special case of the general expressions presented by Reference Faria and HutterFaria and Hutter (2002), viz.

where d₂n denotes an area element on the spherical surface of S ^₂ . If we make use of Equation (31) and the following theorem (for a proof see Reference FariaFaria, 2001)

valid when the arbitrary field a* is continuous in S₂, we can derive, after integration of Equations (20–23) and (30) over the whole orientation space, the ordinary balance equations (1–4), together with the habitual dislocation balance equation

6. Conclusion

The view of polycrystalline polar ice as a mixture with continuous diversity has yielded the construction of a thermodynamic continuum approach to ice-sheet dynamics, which constitutes the first step towards a thermodynamically consistent constitutive theory for large polycrystalline ice masses. In contrast to the fundamentals of other macroscopic models (Reference Morland and StaroszczykMorland and Staroszczyk, 1998; Reference Staroszczyk and MorlandStaroszczyk and Morland, 2001), the present theory naturally incorporates recrystallization processes into its balance equations and blends the virtues of an orientational description of the c-axes distribution with the convenience of a macroscopic description.

Compatibility with microscopic models (e.g. Reference Azuma and Goto-AzumaAzuma and Goto-Azuma, 1996; Reference Castelnau, Thorsteinsson, Kipfstuhl, Duval and CanovaCastelnau and others, 1996; Reference Meyssonnier and PhilipMeyssonnier and Philip, 2000) is ensured by the absence of restrictive constraints upon stresses and strains on the crystalline level, while theories based on the concept of an orientation distribution function (Reference Svendsen and HutterSvendsen and Hutter, 1996; Reference Gödert and HutterGödert and Hutter, 1998; Reference Gagliardini, Meyssonnier, Hutter, Wang and BeerGagliardini and Meyssonnier, 1999) are also encompassed, in fact with advantages. For instance, it is easy to grasp that, in the present theory, the counterpart of the orientation distribution function (ODF) is the mass fraction f* = ϱ*/ϱ. Due to Equation (31), f* is always a normalized function, even when recrystallization is active, in contrast to the situation found in usual ODF-based theories (cf Reference McconnelGödert and Hutter, 1998; Reference Gagliardini, Meyssonnier, Hutter, Wang and BeerGagliardini and Meyssonnier, 1999).

Concerning recrystallization effects, explicit expressions for the dislocation production 77 have recently been proposed (Reference GromaGroma, 1997; Reference Acharya and BeaudoinAcharya and Beaudoin, 2000; for the particular case of ice, see Reference Montagnat and DuvalMontagnat and Duval, 2000), which can likely be extended to its orientation-dependent correlate, 77*. The consequences of these expressions for the interaction and production terms occurring in Equations (20–24) are under current analysis, but it seems certain that some concepts from the thermodynamics of continuous reactions (Reference Aris and GavalasAris and Gavalas, 1966) should play a fundamental role. For instance, notions akin to that of chemical potential and activity were formerly helpful in clarifying the role of interaction terms like and ξ* (Reference Condiff and BrennerCondiff and Brenner, 1969; Reference FariaFaria, 2001), while the latent heat of mass transformation (or, rather, of recrystallization) has recently been recognized as an essential concept for linking the productions of energy, ε* , and entropy, 6* , to that of mass by recrystallization, (cf. section 5).

Finally, we remark that mean grain-size may appear in the actual description just as a recrystallization parameter, reflecting the assumption that the thermomechanical response of ice should not depend directly on it (Reference Budd and JackaBudd and Jacka, 1989). Notwithstanding, an extension of the theory, including grain-size as an additional distinctive grain property, is straightforward and will be considered in the near future.