1. Introduction

Ice-core data from ice sheets reveal the evolution of distinct fabrics with increasing depth. Consequently, according to its depth, ice exhibits a different viscous response. Therefore, following Reference LliboutryLliboutry (1993), a realistic simulation of ice-sheet flow first requires a model to predict the ice fabrics. Then the mean material properties can at least be estimated from those of the single crystal by homogenization. The homogenization procedure proposed by Reference LliboutryLliboutry (1993) was based on a continuous axi-symmetric distribution of the c axes, where the axis of revolution, as well as the distribution itself, were assumed to be known a priori, for instance from borehole data. Recently and more general, polycrystalline ice has been modelled by computer-based models, where the polycrystal is represented by a finite (commonly a large) number of single crystals, e.g. Reference Castelnau, Thorsteinsson, Kipfstuhl, Duval and CanovaCastelnau (1996a, Reference Castelnau, Duval, Lebensohn and Canovab), so that the fabric may freely evolve as a result of the applied loading conditions. Unless those models are found to be in agreement with laboratory observations, incorporation into numerical simulations of the motion of large ice masses is questionable. Therefore, the aim of this work is to propose a model which can describe the evolution of the mechanical properties of polar ice due to deformation-induced fabric evolution, where considerable reduction of the number of unknowns is achieved by adopting a continuum mechanical description to the fabric evolution. In so doing, we extend the classical continuum description, in the following referred to as the large- or ice-sheet-scale formulation, by the so-called small-scale description, which accounts for the c-axes orientations of the polycrystalline aggregate. To achieve the transition from small-scale to large-scale entities, we adopt the concept of a region of influence around each material point of the continuum. If its position, size and shape are neglected, the grain may be identified uniquely by its,-axis unit vector. If, moreover, a statistical description is applied, each large-scale material point may possess all,-axes orientations, which are thus describable by means of a continuous (and differentiable) orientation-distribution function (ODF), f(x,n,t), where × is the position of each point, n denotes a point on the unit sphere S² associated with the orientation of a particular crystal, and t is the time (see, e.g., Reference ClementClement, 1982). Obviously f(•,n,•) = f(•,–n,•), and to pass from an arbitrary small-scale quantity, Xn, to its large-scale counterpart, X, a volume-averaging procedure is applied as in Equation (1) below.

Let x and X be the positions of a material point of the macro-continuum in its present and reference configuration, respectively. Then, the motion of the macro-continuum, i.e. the mapping × = (X,t), is obtained as the homogenized small-scale motion of the crystals n(X,n,t), i.e.

where d²n = sin 0 d dΦ and S2 = 0 : [0,π] X Ψ : [0,2π]. Given the small-scale deformation gradient Fn, its multiplicative decomposition yields

in which R_n, n, and n₀ denote the inelastic (plastic, viscous) deformation, rotationFootnote ¹ , basal-sliding vector and the c-axis unit vector in the reference configuration, respectively. Material time differentiation then yields the small-scale actual velocity gradient, which may be decomposed into its symmetric, and skew-symmetric, parts respectively, e.g. Reference DafaliasDafalias (1984). As in Reference Godert and HutterGodert and Hutter (1998), describes the evolution of the individual c-axis orientation; if W n = W is assumed, we obtain

2. Balance Laws, Constitutive Equations and Motion

In addition to the classical mechanical balance laws for mass, linear and angular momentum, balance of crystal orientation must be considered. Depending on the underlying large-scale formulation, the local orientational balance takes the form

with respect to Eulerian, , and Lagrangian, f, coordinates, respectively. The divergence operators div_× and div_n are referred to the large- and the small-scales, respectively, and production (recrystallization) as well as small-scale conductive flux of orientation are neglected. The material under consideration is characterized by constitutive assumptions for the small-scale state variables, Z _n, which are generally given as a function, , of the large-scale velocity gradient, L, the Gauchy stress tensor, T, the structure tensor of order 2k, and the k-th order small-scale structure tensor ,

with

If no tensor order is explicitly given, the tensor is assumed to be of second order. In this work Zn may be identified with the stress Tn, the strain rate D_n, or it takes the role of the small-scale orientation rate n. However, only the simplest cases will be considered, where either velocity gradient (Voigt-Taylor, VT) or stress (Sachs-Reuss, SR) is assumed to be uniform within the region of influence.

It is known (see, e.g., Reference Castelnau, Duval, Lebensohn and CanovaCastelnau, 1996b), that the contributions of the pyramidal- and the prismatic-slip planes are quite small compared to sliding along the basal plane, which represents the most active glide plane. This prompts us to reduce the single-crystal deformation to basal slip only. The exact constitutive relation for the stress deviator according to Glen’s flow law is hardly amenable to an analytical treatment of the small-scale flow.

In view of many linear models already introduced into the literature (cf Reference Meyssonnier and PhilipMeyssonnier and Philip, 1996; Reference GagliardiniGagliardini and Meyssonnier, 1999), reasonable insight can be anticipated if an assumption of Newtonian fluid is applied. Furthermore, the deviatoric-stress level the fabric develops may be assumed quite low, so we finally reduce the constitutive relation between the resolved stress and the basal gliding to a linear relationship, . With that, the evolution equation (3) can be written more specifically for the Taylor or the Sachs assumption, if in (3), D n is replaced by D or , respectively; μ represents the basal fluidity.Footnote ² Rewriting the resolved quantities for the basal plane by , where the double dot denotes twofold contraction, one may deduce the fine-scale constitutive equation for the VT and the SR assumptions

respectively, where

denotes the projection of a second-order (tensor) quantity onto the basal plane. From this, the large-scale constitutive relations are simply obtained if the small-scale orientation tensors are replaced by their associated large-scale structure tensors, which yields

Because of the underlying small-scale incompressibihty, C is semi-definite, in the sense that it maps only the deviatoric part of a tensor uniquely. In principle, Equations (8) and (3) suffice to compute the fabric evolution due to large-scale loading. However, no influence of the grain-to-grain interaction has yet been taken into account.

One possibility to account for grain-to-grain interactions is to extend Equation (8) by an isotropic part, controlled through a parameter, /3, which represents a large-scale quantity, used to accommodate the small-scale flow, as was done by Reference GagliardiniGagliardini and Meyssonnier (1999). On the other hand, one may extend the small-scale flow directly, so that instead of the pure SR or VTcondition, a linear combination of the SR and VT is applied, weighted by a parameter v, .

This may be derived if we assume that the grain deformation, Dn, which is so far given solely by basal slip, also contains a certain amount of the large-scale deformation D,

where the grain-interaction coefficients v1, v2 are assumed to depend on the actual c-axes distribution. This assumption may be justified by deformation processes which are not describable by means of basal slip alone (like polygonization and the activation of non-basal-slip systems). In the following, the strength of fabric of the polycrystal is measured by

where M_max denotes the maximum eigenvalue of the second order structure tensor M and d {2, 3} gives the dimension of the underlying large-scale space. Hence, perfect alignment along a single direction m may be identified with M_a = 1, whereas for a random distribution of the c axes, M_a = 0. Homogenizing Equation (9), the large-scale constitutive equation takes the form

On the other hand, if M, 1, the structure tensor M → m ⊗ m and Equation (11) should reduce to the single-crystal constitutive equation, C σ C_m, with respect to m. Assuring further the positiveness of the dissipation rate, M _a → l implies

However, as a consequence of increasing alignment, the prerequisite assumptions, justifying the homogenization procedure, are violated, i.e. a small scale no longer exists and all crystals undergo nearly the same deformation D_n→D, which, in terms of the interaction coefficients, yields Vl→ 1 v2→0. In reality, perfect alignment was not observed (cf Reference AlleyAlley, 1992; Reference Castelnau, Thorsteinsson, Kipfstuhl, Duval and CanovaGastelnau and others, 1996a), so that v1 is restricted by a positive number a ≠ 1. Moreover, supposing that the limit case holds even if the alignment is imperfect, v2 is given by

From Equations (9) and (13), the modified small-scale constitutive equation (3) takes the form

where, for convenience, the subscript of v1 is omitted.

So far, the small-scale quantities are all parameterized by means of large-scale quantities and (see Equations (14), (11) and (8)). Assuming them to be stepwise constant, one may integrate n, which yields a functional relation for the motion of the small-scale material point in the form

in which the orientation at the present time, t, is related to its value, n₀, at the initial time, t₀, and to position x₀. Assuming that the number of c axes within a material region of S2 is preserved, the existence of the inverse mapping γ⁻¹ (relative to n and n₀), yields

for the present ODF, which is then constructable accordingly from its initial values if the small-scale motion is known, as we will see for planar flow in the following.

3. Planar Flow

In order to keep the theory considerably simple and to obtain analytical results as far as possible, we restrict considerations to planar deformation, v ε IR2, as one would expect in the vicinity of an ice divide, but fully three-dimensional small-scale flow, n ∊ 7 s2 (where Ts2 ∊ IR3 denotes the set of tangential spaces to S²). If we identify e₃ with the anti-plane direction, the orientation can be written as n = sin ⊖ e_R + cos ⊖e₃, where e_R := cosΦe1 + sinΦe₂, with latitude angle 6 and longitude angle $. Following the calculations outlined in the Appendix, one may derive the following general expression for the ODF,

which depends on the in-plane orientational distribution fs1 and one additional parameter with so that the integration along 6 may be carried out analytically, fs1 describes the distribution of the c axes along the circle S¹ = S² ∩ IR² and is given as the mean value of fs2 with respect to 6, (Equation (18)). Therefore, to obtain the general ODF for this case, it is sufficient to know the in-plane ODF which may be derived if one applies the procedure outlined above to the in-plane motion of the c axes,

The in-plane evolution equation for the c axes is given by (cf. Appendix)

and η = (e₂•W_e1)/r, r + 0. Supposing the parameters r, η and a are constant, the solution of Equation (19) distinguishes four non-trivial cases depending on the value of η,

for which the solutions may be found in integral tables (e.g. Reference Gradshteyn and RyzhikGradshteyn and Ryzhik, 1965). It can be shown that all of them are structurally covered by

where K, δ(u), δ(1) β_(u), /3₍₁₎ are parameters (the indices u, I simply stand for "upper" and "lower") representing the influence of the loading conditions on the solution, and φ = Ψ + a. If one starts from an initially isotropic solution, differentiation of φo = Φo + a with respect to φ (cf. Appendix), leads to

Applying the definition given in Equation (21), the Relation (22) takes the simple form

irrespective of the applied loading conditions. Hence, we expect that Equation (23) represents the general form of the solution, where the parameters a^T = (a, b, c) are functions of the actual large-scale state variables, L, T and M, acting during the time increment At, i.e. the structure of Equation (23) has to be preserved even if one starts from an arbitrarily chosen fabric. For this case the actual ODF is given by where we generalized considerations to an arbitrary time-step At = tk - t _k. To obtain in terms of the actual orientation angle Ф(t_k ) – Фt_k−1 has to be expressed through the inverse of the in-plane motion, . Straightforward calculation then yields the initial in-plane ODF in terms of the actual orientation φ ^k,

For brevity, the following abbreviations are introduced:

with V_j=u or V_j = l, respectively.

Now, combining Equations (24) and (22), one may extract ^T = (l,sin(2$),cos(2$)), so that the actual ODF takes the form

There is no other dependency of (26) on Ф than the explicit dependence through Ф, so that (23) represents indeed the exactFootnote ³ solution of the planar flow problem and provides the rule to update the coefficient vector a successively,

with

where the constant updating matrices depend on the large-scale stress and deformation conditions acting during the underlying time-step, At, as well as on the present fabric. A further reduction of the set of three parameters (a, b, c) to just two independent parameters is possible if the probabilistic nature of the ODF, the constraint condition is taken into account.

To demonstrate how this approach works, we consider the pure shearing or biaxial deformation, d₁₂ = 0, where no large-scale spin occurs, so η = 0 for VT and SR conditions, respectively. If we further assume that positive loading is acting along the 1–direction, 2α = π/2 and r = –dn11 ≤ 0, so that the resulting evolution equation (19) is given by

After the integration, the in-plane motion takes the form

from which one obtains

with δβ = δuβ1 - δ1 βu. According to Equation (27), combination of the above quantities yields the transition matrix for pure shear,

where we used

Starting from an isotropic distribution, for which the coefficients are given by a ^T = (1, 0, 0), the in-plane ODF takes the form

which is exactly the same as if we had integrated Equation (28) directly. As one would expect, for d_n11 > 0 and 4> ε [0, 4 the maximum of fs1 for pure-shear is located at 4>_max = π/2; alternatively, if d_n11 ≤ 0, it will be found at Ψ _ma× = 0.

As a final demonstration, the trivial case, when r → 0, will be considered. The fine-scale differential equation is then given by with ω₁₂ = e₁. The resulting inverse (rigid-body) motion takes the form ΦK–1=ΦK+ω₁₂Δt, from which the update operator

may be derived by expanding the trigonometric function of Equation (23). As one would expect, turns out to be a pure rotation matrix operating on the ODF-parameter-space, where the a–parameter Oaxis serves as the axis of rotation, since a represents the isotropic part of the ODF, and hence remains unchanged under any rotation.

In summary, the actual ODF is given if the actual coefficient vector a is known. Therefore, the evolution of the ODF is given by the evolution of a.

4. Numerics

The finite-element approximation, to which the results below are referred, is based on the elastic-viscoelastic analogy using rectangular four-node-quasi-incompressible elements (cf Reference HughesHughes, 1987). Usually the ice flow is described by the Stokes equations referred to Eulerian coordinates. Hence, in addition to the (local) material (Equation (27)), the large-scale fabric evolution must be considered, i.e. reformulation of / S1 in terms of large-scale Eulerian coordinates (indicated by a tilde) becomes necessary if flux across the element boundaries takes place. Note, that although / s1 was originally defined with respect to Lagrangian coordinates, the Eulerian description of the ODF must reflect the same structure as Equation (23).

This will be done by adopting an explicit time-integration scheme based on a finite-volume approximation (see, e.g., Reference Ferziger and PericFerziger, 1996) using Equation (4), where the finite-volume mesh is chosen to coincide with the finite-element mesh. This enables us to make subsequent use of the velocity field, which was obtained from the finite-element analysis for the fabric computations.

In view of an explicit time-integration scheme, based on a finite-volume approximation, the orientational balances (Equation (4)) will be rewritten in their conservative forms with respect to a certain volume element Vij (ij are counting indices for the central position of the rectangular elements) which yields

where the superposed dot denotes the material and ∂_t() the local time derivative, respectively. Furthermore, ∂Vij represents the element boundary and u its outward unit-normal vector. Recalling that the initial "values" of and f are the same at each time-step, the application of the mean-value-theorem yields

with K ε [0,1], where the ODF approximation is supposed to be C° continuous. Explicit integration is obtained for Equation (33) if K = 1

Note that is still a function of*, which is uniquely defined by three parameters, for each element. So, these coefficients may be calculated in general from three independent values of the actual ODF. For example: if we choose (Ф₁, Ф₂, Ф₃) = (0, π/6, π/3) the coefficient matrix for each element is determined by

Thus, the calculation of the fabric reduces to the determination of three additional equations for each element. At the boundary between two adjacent elements, the ODF is defined as the linear interpolation of the two adjacent ODFs.

Ignoring recrystallization, orientation is conserved,

and the set of additional equations experiences a further reduction by 1 if Equation (35) is adopted and is then in agreement with the classical theory of planar orthotropy, which is uniquely determined by two parameters, one direction and one additional material parameter.

The original domain is to be surrounded by additional finite-volume elements, so-called boundary-volumes, to apply boundary conditions for the fabric by prescribing the coefficient vector a for each of them. In principle, we distinguish three different types of boundary volumes (see Fig. 1): free boundary, across which material as well as fabric transport occurs; stationary boundary, across which the ODF is assumed to be continuous, however material transport may occur; fixed boundary, across which neither material nor fabric may flow. In our problem, the free surface of a glacier partly represents a free boundary, where because of the randomness of the orientation of the ice (snow) crystals, a = (1,0,0) will be prescribed, whereas if the streamlines are directed outward of the free surface, the fabric within the domain will not be influenced; therefore this part should be modelled as a stationary boundary. Due to symmetry, the velocity field at the ice divide is restricted to the vertical direction, so no horizontal flux has to be considered. Therefore, similar to the bottom line, it possesses all the properties of a fixed boundary.

Fig. 1. Schematic ice-sheet discretization by a coupled finite-element volume approach.

5. Application

In this section some results for stationary planar and axisymmetric flow will be presented. These results are used to identify the grain interaction parameter v (see Equation (14)), by accommodation of the computational results to measurement data of the GRIP core as well as to those obtained by Reference Castelnau, Thorsteinsson, Kipfstuhl, Duval and CanovaGastelnau and others (1996a) using a so-called viscoplastic-self-consistent approach (VPSG).

To this end, the ODF is needed for axi-symmetric loading conditions, which is also covered by Equation (17) if one considers that If the stresses T = σ e₃ ⊗ e₃ are introduced into Equation (6)2, the large-scale strain rate, d₃₃ = e₃ (cf Equation (8)), takes the form

where M₃₃ and M₃₃₃₃ are the components associated to the loading axis of the second- and fourth-order structure tensors, respectively. Recalling, that the accumulated strain ε_c is given hereby

it can be written completely in terms of the applied loading if the large-scale deformation gradient F = γe ₁ ⊗ e ₁ + γe ₂ is taken into account;

In Equation (38) as well as in the following, uniform effective loading is assumed, i.e. σ_iμ_iΔt_i = σμΔt. In addition, the evolution of the large-scale structural quantities in the course of the time-step is neglected. On the other hand, after treating the first-order small-scale evolution equation (14) and using (36), one obtains

which represents the anisotropy of the polycrystal. In contrast to A, the anisotropy of the GRIP core data, displayed in Figure 2, is related to where the integral is to be taken over the hemisphere H² and || • || denotes the Euclidean norm. Accomplishing the integration, one obtains

Fig. 2. Evolution of degree of orientation R_o for pure SR (α = 0.0) and grain interaction ((3=0.5, α = 1.2) in comparion to GRIP data and VPSC model as a function of the accumulated linear strain, ∊_c, under uniaxial compression.

which relates both parameters. If the x-axes distribution reflects rotational symmetry, Ro ₀ → 1 if a single maximum fabric develops, whereas, if the c axes are distributed randomly, R₀ = 0. Hence, R₀ gives the intensity of the actual fabric. We already mentioned that we prefer M a as an alignment measure, however we will give results with respect to Ro to make them comparable to Reference Castelnau, Thorsteinsson, Kipfstuhl, Duval and CanovaGastelnau and others (1996a). Accordingly, the grain interaction coefficient is supposed to take the general form

where α and β are constant values, determined by comparison of the computational results and the GRIP data (cf Fig. 2). From this, best agreement was achieved with /3 = 0.5, which is assumed throughout the following considerations. That is, the identification of the effective grain interaction, v, is reduced to the determination of a. Note, that the data presented in Figure 2 are based on an assumed velocity field along the GRIP core, which certainly does influence the fabric evolution.

Despite of this, it is obvious from Figure 2 that the curve obtained with an interaction parameter α = 1.2 yields the best fit to the data, except when R₀ ≤ 0.4, for which the pure SR as well as the VPSG model give better approximation. Results (not presented in Fig. 2) very close to VPSG are obtained if α = 1, where if all crystals are aligned, the small-scale evolution equation is completely based on the VTassumption. On the other hand, if a > 1, the SR based contribution to the small-scale flow changes its direction when so that the fabric will not approach perfect alignment, i.e. R₀ ≤ 1. It is known from the work of Reference AlleyAlley and others (1992), that this behaviour used to be explained by polygonization (or rotation recrystallization), where the strain-induced small-scale motion appears to be balanced through the development of sub-grain boundaries, which under uniaxial loading is expected to result in a so-called girdle fabric. Considering the periodicity of the ODF, (Equation (17)), a girdle fabric cannot be taken into account exactly. Hence, beside controlling the small-scale directly by the actual fabric, the interaction parameters, a, β, also summarize effects due to deviations of the fabric from a single-maximum fabric.

Next, the fully planar case, i.e. , is considered, where the material structure is completely determined by only one element of the fourth-order structure tensor, .

Figure 3 is concerned with the material response under simple shear deformation K, with . Owing to the underlying linear stress-strain-rate relation, the shear viscosity, which after an initial growth decreases until it reaches half of its initial value, reveals softening behaviour as is qualitatively expected (Reference DuvalDuval, 1981). Moreover Figure 3 shows that simple shear deformation requires the development of normal stresses σ = (σ₁₁ – σ₁₂) similar to the shear stress response (cf. Reference Li, Jacka and BuddLi and others, 1999). It should be mentioned that the peaks of the stresses are strongly influenced by the choice of the interaction parameters a, β, e.g. if no grain interaction is considered at all, one obtains strong hardening behaviour.

Fig. 3. Normalized shear stress rand normal stress σ for simple shear deformation, where α = 12, β=0.5.

This can be deduced from the results presented in Figure 4, where the maximum of the ODF, f_m, is plotted against its orientation, Φm driven by simple shearing for different parameters a. The computation starts from a randomly distributed c-axes configuration. In spite of a, the maximum value f_m initially evolves at . The density of the dotted fines may be associated with the convergence of the material behaviour, i.e. a low dot density reflects fast changes on the small scale. Whereas for α = 1.2 and α = 1.1, f_m posseses the limit values f_m ≪ 16 and f_m ≪ 30, respectively, infinite growth takes place if a ∈ [0,1.0], where the pure SR model is represented by α = 0.0. Moreover, it is known that, under the applied loading, the so-called easy-glide configuration of a polycrystal is achieved if the mean, -axes orientation, which may be identified with Ψ _m, is directed close to the vertical axis, 4>_m = π/2. Therefore, the softest material behaviour under simple shear is obtained if α = 1.2; hardening (locking) was observed for a ≤ 1. On the other hand, approaching the vertical axis (Ф_m → π/2) is accompanied with a decrease of f_m and hence, a decrease of the effective anisotropy.

Fig. 4. Maximum of the ODF, f_m, vs its orientation, 4>_max in radians for pure SR (α = 0.0) and different grain interactions (α =1.0, α = 1.1, α = 1.2).

Finally, a rectangular domain, 10 times longer than high, was considered as a simple model for stationary plane flow in the vicinity of an ice divide, driven by gravity. The vertical direction was discretized by 10 equal size elements with aspect ratio 0.25. Perfect sliding was assumed at the bottom in conjunction with a vanishing slope. Symmetry implies pure vertical flow at the divide. At the free surface the ODF was expected to reflect isotropic behaviour, whereas at all other boundaries stationarity was assumed. Note that these large-scale boundary conditions correspond on the material level to uniaxial compression used to accommodate the interaction parameters α = 1.2, (3 = 0.5 (cf Equation (41)). Therefore, we may anticipate that the GRIP data will to some extent be reproduced by the simulation. Comparison of the evolution of the order parameter M₁₁₂₂ at the divide with the quadratic Schmidt factor ( obtained from the GRIP core, is given in Figure 5. As M₁₁₂₂, the Schmidt factor is a measure of the softness (fluidity); a significant effect of the orientation distribution on the mean material response is expected only below «1250 m (Thorsteinsson, 1997, p. 113). Accordingly, this depth was chosen to correspond to the height (H = 1) of the ice-sheet model, where M₁₁₂₂ = ⅛. The quadratic Schmidt factor was accommodated by linear transformation to the isotropic case. The graph in Figure 5 was selected from the series of solutions to represent a good fit to the GRIP data. Hence, it is in principle possible to obtain reasonable results by applying this approach to modelling ice-sheet flow.

Fig. 5. Evolution of alignment along the GRIP core. Finite-element-finite-volume computation (solid line) and measured data (dotted line).

6. Conclusion

A decoupled, two-scale (small/large) description of the strain-induced orthotopic anisotropy of ice polycrystals has been proposed. The evolution of the ice-crystal c axis, originally derived from pure kinematic calculations, was supplemented by small-scale phenomenological considerations to account for grain-to-grain interaction. Under stepwise plane-flow conditions, it could be shown that the fabric, represented by an orientation-distribution function, is completely determined by three parameters. A general transition matrix was derived from the fabric-evolution equation establishing stepwise computation of the fabric parameters.

It could be shown that this approach is capable of representing the fabric evolution of polycrystalline ice, at least qualitatively. A strong influence of the material response on the grain interaction was observed, but could not be considered in detail here. By comparison with measurement data, the importance of a grain-interaction contribution to the fabric evolution became obvious. This was finally done by combining the uniform-stress and the uniform-strain assumptions through a phenomenological interaction coefficient depending on the actual fabric alone. The specific accommodation of the interaction model to field data revealed its ability to model the material response due to effects normally explained by polygonization/rotation recrystallization and the activation of non-basal-slip systems.

To obtain the actual fabric in a large-scale ice-sheet flow, we generalized the update procedure to Eulerian coordinates. The flux of orientation was considered by applying a decoupled finite-element-finite-volume approach, where the actual structure is represented by only three additional material parameters at each integration point to account for the evolving orthotropy. To concentrate on the proper material response, not influenced by a varying free-surface as well as by a wavy bedrock, a rectangular domain was chosen with an initially overall random distribution of the c axes as the most simple case of nearby ice-divide flow. Even if a coarse discretization is applied, the proposed strategy turned out to be appropriate for producing reasonable results.

Finally, it should be mentioned that there is further need for more elaborate investigations, especially concerning alternative initial and boundary conditions as well as different interaction models.