Implementing an empirical scalar tertiary anisotropic rheology ( ESTAR ) into large-scale ice sheet models

The microstructural evolution that occurs in polycrystalline ice during deformation leads to the development of anisotropic rheological properties that are not adequately described by the most common, isotropic, ice flow relation used in large-scale ice sheet models – the Glen flow relation. We present a preliminary assessment of the implementation in the Ice Sheet System Model (ISSM) of a computationally-efficient, empirical, scalar, tertiary, anisotropic rheology (ESTAR). The effect of this anisotropic rheology on ice flow dynamics is investigated by comparing idealised simulations using ESTAR with 5 those using the isotropic Glen flow relation, where the latter includes an overall flow enhancement factor. For an idealised embayed ice shelf, the Glen flow relation overestimates velocities by up to 17% when using an enhancement factor equivalent to the maximum value prescribed by ESTAR. Importantly, no single Glen enhancement factor can accurately capture the spatial variations in flow over the ice shelf. For flow-line studies of idealised grounded flow over varying topography or variable basal friction – both scenarios dominated at depth by bed-parallel shear – the differences between simulated velocities using ESTAR 10 and the Glen flow relation vary according to the value of the enhancement factor used to calibrate the Glen flow relation. These results demonstrate the importance of describing the anisotropic rheology of ice in a physically realistic manner, and have implications for simulations of ice sheet evolution used to reconstruct paleo-ice sheet extent and predict future ice sheet contributions to sea level.

mechanical response (the strain-rate state) does not change by rotation of the sample as, by definition of ESTAR, the fabric instantaneously adapt to the stress state, and thus to any rotation of the stress state.The dependance of the flow law to the stress state, results from the mechanical anisotropy of the ice crystal and the anisotropy of the crystal orientations (the fabric) that develop during tertiary creep; but the mechanical behaviour predicted by ESTAR is isotropic.Note that this property is not due to the scalar relationship, the CAFFE model (e.g.Seddik et al., J. Glaciol. 2008) is anisotropic as its scalar enhancement factor depends on the direction of the stress state with respect to the ice fabric, and thus the mechanical response depends on the orientation of the stress state relative to the ice sample.
We are certainly presenting ESTAR as a physically-based alternative to the standard Glen flow relation.
The reviewer is incorrect in stating that ESTAR is isotropic.The conceptual origin of the description of anisotropy in ESTAR is based on considerable observations from ice cores, boreholes and laboratory experiments (e.g.Jacka and Maccagnan, 1984;Dahl-Jensen and Gundestrup, 1987;Gao and Jacka, 1987;Budd and Jacka, 1989;Etheridge, 1989;Li and Jacka, 1998;Morgan et al., 1998;Wang and Warner, 1999;Treverrow et al., 2012).Observations such as these demonstrate that if the stress configuration is changed within a fixed reference frame then the compatible fabric that develops in response to that stress is correspondingly altered, i.e. specific fabric patterns are associated with specific stress configurations (note that in the following we use the term 'fabric' to describe the distribution of crystallographic c-axis orientations).This is valid provided that the stress configuration remains stable for a period that exceeds the time (or strain) required for the microstructure to adapt to that state of stress.Since polycrystalline anisotropy is fundamentally linked to microstructure (including fabric, but not exclusively so) the anisotropic response of an aggregate can be directly related to the present stress configuration, if that is a suitable indication of a persistent recent strain history.
Identification of the local non-rotating shear plane (normal, n), the simple shear stress, τ , acting on that plane and its magnitude relative to the total effective stress, τ e , summarise the computational framework for specifying the anisotropic enhancement factor, E(λ S ), in ESTAR (Equations 4-6).The parameter λ S also characterises the fabric pattern, i.e., when λ S = 1 there will be a strong single maximum, while when λ S = 0, there will be c-axes concentrated around a conical surface centred on the compression axis -a small circle girdle.The anisotropy of E(λ S ) is derived from the spatial variability in the flow field and its dependence on the deformational component of the vorticity.
The concepts of rheology and its description in a flow relation can sometimes become confused.ESTAR is the representation of the deformational response of ice with a crystal fabric anisotropy compatible with prolonged deformation under a stable stress regime -this is the 'tertiary' aspect of ESTAR.When the compatibility requirement is considered, it is clear that there is no sense in contemplating arbitrary rotations between the stress configuration and anisotropic crystal fabricsthis relates to the deformability of anisotropic ice.Regarded in this light, it seems to us that the suggestions of both reviewers that one could meaningfully consider an arbitrary rotation of spatially anisotropic material properties and the applied stresses breaks the intrinsic point of ESTAR.The point also seems semantic -the rheology of polycrystalline anisotropic ice versus an anisotropic rheological description for polycrystalline ice.
This reviewer notes correctly that the enhancement factor in the CAFFE flow relation is anisotropic since it is defined by the local stress field and its orientation relationship with respect to a polycrystalline aggregate, and rotation of one relative to the other produces a different effect.In CAFFE the explicit inclusion of fabric in defining the deformability of polycrystalline ice allows the instantaneous flow response for any arbitrary combination of stress and fabric to be calculated.CAFFE and ESTAR are structurally similar in featuring a collinear relationship between stress and strain rate tensor components and hence each can be interpreted as an enhancement factor, defined as a function of some experimentally-based enhancement parameters and the local stress field.A key difference between the two relations is that for ESTAR fabric is not explicitly included in specification of the enhancement.By linking the enhancement directly to the relative magnitude of the components of σ via E(λ S ), ESTAR is restricted to a subset of the fabric/stress combinations that can be contemplated with microstructure based approaches, such as CAFFE.Importantly, these combinations where ESTAR is applicable are those where fabric has evolved to be compatible with the local stress field.Put another way, ESTAR only allows contemplation of tertiary creep (addressed on P2L5-13; P5L5-18; P7L1-4).This is not a problem since tertiary creep is relevant to the vast majority of the most dynamically active regions of an ice sheet.In the revised manuscript, we have noted examples where this might not be the case (P5L21-P6L12): "There are of course zones within an ice sheet where the assumption of compatible tertiary flow will not apply; however, we note that these zones will be restricted in their extent.We contend that ESTAR will apply to the vast majority of the dynamically active regions of an ice sheet, in particular those zones where creep deformation makes a significant contribution to the overall flow.Specific zones where the assumption of tertiary creep may be inappropriate can be summarised as those where fabric has not yet evolved compatibility with the flow, where there is a rapid transition in the flow configuration, or where creep deformation makes only a minor contribution to the overall dynamics.For example, in very cold ice in a low stress setting, such as the uppermost layers of the polar ice sheets, the time required to accumulate the strain necessary to develop a compatible fabric may lead to a near-surface zone in which the assumption of tertiary creep is not valid.Since the development of anisotropic fabrics provides an indication of the existence of tertiary flow, their observation at modest depths, (e.g., 100 − 200 m; Morgan et al., 1997;DiPrinzio et al., 2005;Treverrow et al., 2016) allows estimation of the maximum extent of the zone where tertiary creep is not occurring.The observation within polar firn of microdeformation processes that are necessary for the development of fabric throughout ice sheets (e.g.Kipfstuhl et al., 2009;Faria et al., 2014) suggests that it may be appropriate to even further restrict the extent of the nearsurface zone for which the assumption of tertiary creep is not valid.Additionally, the nonlinear nature of polycrystalline ice rheology leads to very high viscosities in low temperature and stress environments, so that incorrectly estimating deformation rates due to the assumption of tertiary flow in such regions may be of limited importance to simulations of ice sheet evolution.
Regions where rapid transitions in dynamic conditions can lead to abrupt changes in the pattern of applied stresses and a potential breakdown in tertiary flow compatibility include ice shelf grounding zones and other locations where basal traction is lost or abruptly changes, e.g., where ice flows over a subglacial lake, or with the onset of basal sliding in ice streams.The convergence zones where tributary glaciers or ice streams merge with a larger flow unit at a high angle may also lead to a transition in dynamic conditions that is problematic for the assumption of tertiary compatibility.Of course the more highly dynamic the evolving flow regime, the more rapidly a new compatible anisotropy will be established, so that the spatial interval where the flow relations are inapplicable may be limited.While there is little guidance on how to extend empirical flow relations to parametrise ice rheology in these transition regions, we note that similar difficulties exist for a Glen-type flow relation, which unlike ESTAR does not have the benefit of being able to correctly describe anisotropic enhancement throughout the remainder of the ice sheet."That tertiary creep is relevant to the majority of the ice sheet is crucial and is dealt with further in responding to Point 2.
Point 2: Tertiary creep R1: In the paper tertiary creep is presented as "the predominant mode of deformation in ice sheets", under the justification that it is usually reached after few percent of deformation in laboratory creep experiments.For a detail review on the ductile deformation of ice and textures in polar ice masses one can refer to the book "creep and fracture of ice" by Schulson and Duval.Tertiary creep in laboratory experiments is explained by softening processes associated with migration recrystallization.Textures resulting from migration recrystallization are indeed stress-configuration depend and the argument that the microstructure evolves more rapidly than the flow configuration sufficiently robust to assume a scalar isotropic relationship in the context of large scale ice flow modelling.In the central part of the ice-sheets, temperature and strain energy can be too low to initiate migration recrystallization and there is no evidence that it occurs except in the warm ice near the bed (the last 100-200 meters usually).See texture measurements in all deep ice cores in the central regions of Greenland and Antarctica.There is a general agreement that in the regions of normal grain growth and rotation recrystallization (most of the ice thickness in the central regions) basically develops as the results of lattice rotation due to intractrystalline slip; and thus reflect the deformation history.
The suggestion being made by the reviewer here seems to be that for laboratory experiments conducted at relatively high temperatures, tertiary creep is heavily influenced by migration recrystallisation, and that the results of such experiments are not representative of polar ice dynamics.The extension of this argument is that at lower temperatures boundary migration (a process which allows microstructural recovery during deformation) is less effective and does not contribute to fabric development.It seems this line of reasoning is then used to discount our expectation that appropriate fabrics, and perhaps tertiary creep conditions will be found in polar ice masses.
Following Schulson and Duval (2009), the reviewer states that fabric development at low temperatures is largely attributed to rotation recrystallisation -a deformation process, without the accompanying recovery afforded by boundary migration.Consequently for different polycrystalline aggregates, deformed under identical conditions of stress, it is suggested that the resultant fabric patterns will vary as a function of temperature due to the activity or otherwise of migration recrystallisation.As per Schulson and Duval (2009) the temperature at which migration recrystallisation is claimed to be sufficiently depressed, so that it no longer contributes significantly to fabric development, is ∼ −15 • C. Differences in the fabrics that are supposed to form in uniaxial compression at various temperatures are provided by Schulson and Duval (2009) as evidence to support this view.The small circle girdle (or cone-type) fabrics that form at higher temperatures are claimed to give way to single maximum (or cluster-type) fabrics at lower temperatures and stresses.The suggested implication of this for ESTAR is that the enhancement parameter derived for uniaxial compression from laboratory experiments is not valid for ice sheets.There are two problems with this argument and these appear to relate directly to the reviewer's concerns about how widespread tertiary creep is within an ice sheet.These issues are addressed in detail below, following a summary reiterating our view, before returning to the rest of the reviewer's comments.
More broadly, as described in Budd et al. (2013), tertiary creep describes a state of deformation where the microstructure within a polycrystalline aggregate has evolved to be compatible with the stress field, and this results in the development of polycrystalline anisotropy.A conspicuous feature of tertiary creep is the development of a fabric and to an extent one can remain agnostic about the actual microdeformation and recovery mechanisms that contribute to the development of this state.The anisotropic nature of the macroscopic response, and its spatial variability, is of primary importance to modelling large-scale ice sheet dynamics, not the specific mechanisms of microdeformation.
In short, we consider tertiary creep to be far more widespread than reviewer 1 claims it to be.Of course there are zones within the ice sheet where tertiary creep is not relevant and the manuscript has been altered to more clearly highlight these (e.g., P5L21-P6L12 -as quoted directly in our response to Point 1 above).

Comment to the reviewer on fabric development and the extent of tertiary creep in ice sheets
First we discuss the temperature range over which small circle girdle fabrics occur in response to uniaxial compression and then we address the issue of whether or not a single maximum fabrics can be associated with uniaxial compression.

Small circle girdle/cone-type fabrics
From laboratory experiments there is considerable evidence (e.g.Jacka and Maccagnan, 1984;Jacka, 1984;Gao and Jacka, 1987;Jacka and Li, 2000;Treverrow et al., 2012) supporting the formation of small circle girdle fabrics in response to uniaxial compression.That this occurs is not in question here, rather it is the issue of to what extent similar patterns have been observed in the polar ice sheets.
Law Dome in East Antarctica provides a good example of an approximately radial dome whose dynamics are mostly isolated from the remainder of the Antarctic Ice Sheet.Observations from the A001 core drilled at Law Dome summit, where the ice thickness is ∼ 1200 m show a distinct small circle girdle fabric at a depth of 318 m (see Fig. 3a, Budd and Jacka, 1989).At this depth uniaxial compression is expected to be dominant and the in situ temperature is ∼ −22 • C. We highlight this observation since it demonstrates that small circle girdle fabrics can form in response to uniaxial compression at temperatures below the ∼ −15 • C threshold for the activity for migration recrystallisation that is suggested by Schulson and Duval (2009).Small circle girdle fabrics have been recorded at many other locations, e.g.cores drilled at Siple Dome (DiPrinzio et al., 2005), Byrd (Gow and Williamson, 1976), Dye 3, Greenland (Herron et al., 1985) and the Amery (Budd, 1972) and Ross ice shelves (Gow, 1963).Importantly, some of these observations correspond to locations where the in situ temperatures were even lower than the ∼ −22 • C described for the Law Dome, A001 example.
Perhaps one reason why small circle girdle fabrics are considered by some to be the preserve of uniaxial compression laboratory experiments conducted at high temperatures and/or stresses is that these fabrics are not frequently identified in ice cores since ideal dome summits are uncommon, i.e., the vertical shortening is typically associated with different horizontal strain rate components in radial directions, particularly at depth.Even at sites that provide a good approximation of a dome surface, the inevitable depth-evolution of the flow configuration leads to correspondingly different fabrics.Indeed, at ice divides a nearly 2D plane-strain flow can be expected.
Under conditions of plane-strain a restricted form of the cone-type fabric pattern is encountered.These are 2-pole fabrics where the individual maxima are aligned in the direction of horizontal extension and separated by the small circle girdle cone angle (e.g.Gow, 1963;Wakahama, 1974).In laboratory ice deformation experiments, similar levels of flow enhancement are observed for plane strain and unconfined uniaxial compression (e.g.Wilson, 1982;Wilson and Peternell, 2012;Budd et al., 2013).
For ESTAR this means E C has a constant value for both plane strain and unconfined compression.Assuming this equality holds for any combination of normal deformations allows the Budd et al. (2013) flow relation (ESTAR) to be generalised via specification of the ζ-parameter.

Single maximum/cluster-type fabrics
As noted above, strongly clustered vertical single maximum fabrics are routinely observed at depth in polar ice sheets.Contrary to what we interpret as the opinion of reviewer 1, we do not suppose that these fabrics form in response to uniaxial compression at temperatures where the deformation regime is dominated by rotation recrystallisation and recovery by migration recrystallisation is negligible.Where strong single maximum fabrics have been observed in polar ice sheets these can be explained by the contribution of simple shear to the local dynamics -this includes ice-coring sites in central regions of Antarctica and Greenland (e.g.Gow and Williamson, 1976;Russell-Head and Budd, 1979;Azuma and Higashi, 1984;Tison et al., 1994;Gow et al., 1997;Morgan et al., 1997Morgan et al., , 1998;;Thorsteinsson et al., 1997;Azuma et al., 1999;DiPrinzio et al., 2005;Durand et al., 2007;Montagnat et al., 2012Montagnat et al., , 2014;;Weikusat et al., 2016).
A primary concern when selecting a site to drill an ice core for a palaeoclimate record is that as far as possible it should be dynamically quiescent.Despite this aim it is clear that considerable shear can exist well above bedrock at the dome summit or ridge divide locations chosen for palaeoclimate studies.At drilling sites where there are measurements of borehole inclination to accompany the corresponding fabric data it is apparent that the development of distinct single maxima fabrics is associated with increasing levels of a simple shear deformation in the overall flow regime at depth (e.g.Dahl-Jensen and Gundestrup, 1987;Etheridge, 1989;Morgan et al., 1998;Jansen et al., 2017).No summit/divide is ideal, nor are they necessarily stationary over glacial cycles.Combined, these factors increase the likelihood of encountering single maximum fabrics (and simple shear) at depth at a palaeoclimate coring site.Even at small distances (1-2 ice thicknesses) from an ideal summit, simple shear is present at modest depths and makes a significant contribution to the overall deformation.
In Greenland the GRIP-NGRIP-NEEM ice cores were drilled along a ∼ 700 km section of flow line.The coring sites are located nearly equidistant along the flow line, with GRIP at the highest elevation and NEEM furthest along the flow line.For the GRIP, NGRIP and NEEM sites Montagnat et al. (2014) have estimated the depth at which the vertical and simple shear strain rates are equal.Their estimates are 79%, 73% and 65% of the total ice thickness, respectively.As noted by Montagnat et al. (2014), these estimates do not make any allowance for the effect of mechanical anisotropy on enhancing strain rates.It is reasonable to expect that doing so would lead to a proportionally greater increase in the shear strain rates relative to the vertical strain rates.This would shift estimates of the depth where the simple shear and vertical strain rates are equal to lower values, i.e. closer to the ice sheet surface.Results from Law Dome provide some interesting context for these estimates since the simple shear strain rates have been quantified by repeat measurements of borehole inclination at several locations and the vertical strain rate can be estimated by similar methods to those employed by Montagnat et al. (2014).At the Law Dome, Dome Summit South (DSS) site an ice core was drilled 4.7 km (< 4 ice thicknesses) SSW of the dome summit.Here the simple shear and vertical strain rates are equivalent at ∼ 1 /2 the total ice thickness (Morgan et al., 1998;Treverrow et al., 2015), demonstrating that considerable simple shear can develop throughout the depth profile very close to a dome summit.In accord with the experimental observations of Budd et al. (2013), observations from ice cores reveal there is a progression towards increasingly strong single-maximum fabrics and relatively higher strain rates as the proportion of simple shear within the deformation regime increases (even when the effect of temperature on flow rates is considered, e.g.Etheridge, 1989;Wang and Warner, 1999;Morgan et al., 1997Morgan et al., , 1998)).These observations are part of the fundamental basis of ESTAR.
Laboratory experiments also demonstrate the formation single maximum fabrics in simple shear (e.g.Bouchez and Duval, 1982;Li and Jacka, 1998;Li et al., 2000;Wilson and Sim, 2002;Budd et al., 2013;Wilson and Peternell, 2012).As described in Schulson and Duval (2009) a two-maxima fabric forms initially during simple shear experiments.The dominant maxima is normal to a plane referred to as the permanent (Bouchez and Duval, 1982) or non-rotating (Budd et al., 2013) shear plane.With increasing strain the second maxima rotates towards the dominant maxima and vanishes (e.g.Bouchez and Duval, 1982).
The only experiments we are aware of which demonstrate that a single maximum fabric can form in response to uniaxial compression have been conducted at very high deviatoric stresses in conjunction with a high confining pressure.Prior et al. (2015) present a diffuse single maximum fabric for a sample deformed in uniaxial compression to a strain of 0.37 with a deviatoric stress of 15 MPa at −43 • C where the confining pressure was 50 MPa.Similar experiments by Qi et al. (2017) at −10 • C with 10 MPa confining pressure show that as the deviatoric stress increases from 1.30 MPa to 4.31 MPa the resultant fabric undergoes a transition from a small circle girdle to a diffuse single maximum.This behaviour is attributed to the increased activity of rotation recrystallisation at higher stresses.Importantly, these experiments do not support for the concept of a single maximum forming in response to compression within an ice sheet, since they clearly demonstrate that a small circle girdle is expected for deviatoric normal stresses ≤ 1.30 MPa.
To be clear, we are not suggesting that temperature and stress have no influence on the relative activity of microdeformation and recovery processes occurring within polar ice sheets and the formation of fabric.A detailed comparison of the fabrics encountered in Antarctic and Greenland ice cores and laboratory experiments is beyond the scope of the present response.However, it is worth noting that factors other than stress (both magnitude and configuration) can influence the development of fabric.For example, the presence of soluble and insoluble impurities is related to fluctuations in fabric and grain size over small spatial scales that are observed in ice cores (e.g.Morgan et al., 1997;DiPrinzio et al., 2005).
The series of uniaxial compression experiments of Jacka and Li (2000) provide some insight into the effects of temperature and stress magnitude on the rate at which fabric and strain rates evolve as a function of accumulated strain.Their results suggest that for temperatures below −15 • C fabrics do not evolve into the distinct small circle girdles that are observed at higher temperatures.For accumulated strains of ∼ 10% (or less) the observed fabrics are weakly anisotropic, displaying a minimal degree evolution from the initially isotropic distribution of c-axis orientations.Unfortunately these fabrics do not provide strong support for the development of either a small circle girdle or single maximum fabric.Jacka and Li (2000) conducted a single experiment at −45 • C, that unlike other experiments in the series, was conducted in two stages to allow higher strains to be accumulated (e.g.> 20%, Treverrow et al., 2012).This experiment provides some support for reduced rates of fabric and strain rate evolution at low temperatures and/or stresses.This experiment was incomplete when Jacka and Li (2000) was published, but nevertheless was showing signs of a transition from secondary to tertiary creep (E>1).At the end of the first stage the sample was machined back to a cylindrical shape prior to recommencing deformation.The experiment was terminated some considerable time after the initial results were published.Analysis of the sample (T.H. Jacka pers.comm) demonstrated that while a weak small circle girdle fabric did eventually form with continued deformation, the strain required for this at −45 • C was higher than the 8-10% observed in the experiments at higher temperatures and stresses.
That lower rates of fabric evolution may be expected at lower temperatures (or stresses) is consistent with the detailed assessment of dynamic recrystallisation within polar ice presented by Faria et al. (2014).They describe how multiple microdeformation and recovery processes are active during microstructural evolution, even in the very cold ice (T < −20 • C) present in the upper regions of the polar ice sheets.Importantly, Faria et al. (2014) note that the relative contribution of these processes to microstructural evolution can vary according to temperature and strain rate.
Since small circle girdle fabrics can form due to uniaxial compression in ice sheets -even if they are not widespread -then allowances must be made for the activity of boundary migration processes at lower temperatures than the suggested threshold values of ∼ −10 • C (Duval and Castelnau, 1995) or ∼ −15 • C (Schulson and Duval, 2009).Observations from polar firn and ice cores, e.g., Kipfstuhl et al. (2009); Weikusat et al. (2009b,a) reveal microstructural features which clearly demonstrate that internal strain energy is sufficient to allow dynamic recrystallisation processes, including straininduced boundary migration, to contribute to microstructural evolution at temperatures below −20 • C. Modelling of the microstructural evolution of firn (Steinbach et al., 2016) demonstrates the crucial role of air bubbles in the development of the strain localisation and strain energies that are sufficient to drive dynamic recrystallisation at low temperatures in firn.
The point of the preceding discussion is to demonstrate that tertiary creep and the development of polycrystalline anisotropy is associated with the formation of a statistically steady state microstructure (including fabric) that persists under stable conditions of temperature and stress.Kipfstuhl et al. (2009) and Faria et al. (2014) demonstrate that tertiary creep cannot be simply defined according to activity (or otherwise) boundary migration.Consequently it is incorrect to invoke this single specific microdeformation process as the determinant of fabric variations observed with depth (or temperature) in ice sheets.Kipfstuhl et al. (2009) demonstrate that strain-induced boundary migration is active below the threshold temperatures proposed by Schulson and Duval (2009).
Regarding the issue of just how widespread tertiary creep is within polar ice sheets, Faria et al. ( 2014) describe (with specific reference to the EDML site) how the densification of firn, driven by a reduction in pore space, must be accommodated by ice deformation.Total accumulated strains within firn are high enough to support the conclusion that tertiary creep is occurring there.Faria et al. (2014) describe how the occurrence of tertiary creep, without the development of a corresponding pattern of preferred c-axis orientations is due to the complex geometry of the pore space in the two-phase ice-air system which produces localised heterogeneity in the stresses driving deformation.Thus, any localised development of preferred c-axis orientations would be masked by the spatial variability in local stress and strain.
There are some counter-intuitive aspects to the argument for single maximum fabrics forming in response to compression.Firstly, for deformation largely accommodated by basal slip, the effective viscosity will continually increase if a single maximum evolves in response to compression.Consequently high levels of deformation would have to be accommodated by non-basal processes.Secondly, if a single maximum forms in compression due to the inactivity of migration recrystallisation, then it is hard to see why single maximum fabrics become stronger with increasing depth in an ice sheet where the ice is is also warmer.If a single maximum did form in an ice sheet due to compression, as the ice becomes warmer we should see a transition to a cone-type fabric, since the activity of migration recrystallisation would be increased.The problem is that we don't see such a progression in fabrics.
The multiple maxima fabrics observed in zones closest to bed in an ice sheet don't provide the answer here, since they are not the same as the cone-type small circle girdle fabrics associated with uniaxial compression.Multiple maxima fabrics form due to the high activity of migration recrystallisation in temperate ice.In polar ice sheets this often corresponds to warm ice near the bed that has undergone a stress relaxation due to the effect of underlying bedrock topography on disrupting the large-scale flow.The high levels of strain energy in relatively warm ice can drive migration recrystallisation.The existence of a topographic disturbance to the flow can be deduced from horizontal velocity profiles obtained from borehole inclination measurements.The inflection in these profiles close to the bed is consistent with boundary layer flows where there is an adverse pressure gradient, i.e. the reduction in velocity close to the stationary surface (the bed) exceeds expectations associated with either a neutral or favourable pressure gradient.
We now return to the latter part of the reviewer's comment about textures and tertiary creep, which we repeat for convenience: R1: Textures resulting from migration recrystallization are indeed stress-configuration depend and the argument that the microstructure evolves more rapidly than the flow configuration sufficiently robust to assume a scalar isotropic relationship in the context of large scale ice flow modelling.In the central part of the ice-sheets, temperature and strain energy can be too low to initiate migration recrystallization and there is no evidence that it occurs except in the warm ice near the bed (the last 100-200 meters usually).See texture measurements in all deep ice cores in the central regions of Greenland and Antarctica.There is a general agreement that in the regions of normal grain growth and rotation recrystallization (most of the ice thickness in the central regions) basically develops as the results of lattice rotation due to intractrystalline slip; and thus reflect the deformation history.
We agree with the reviewer that all development of anisotropic fabric reflects the deformation history, more particularly, the recent deformation history.It is certainly true that within the lowest 100 − 200 m warm ice, high levels of strain energy and stress relaxation due to the effect of underlying bedrock topography on disrupting the large-scale flow can lead to significant migration recrystallisation and formation of multiple-maxima fabrics.The topographic disturbance to the flow can be deduced from horizontal velocity profiles obtained from borehole inclination measurements.Such profiles are consistent with flows where an adverse pressure gradient exists.However, as discussed above, it is incorrect to suggest that there is 'general agreement' on a depth-progression of specific zones of normal grain growth and rotation recrystallisation etc throughout an ice sheet.A more up-to-date review on microdeformation in ice sheets is provided by Faria et al. (2014).The observations of Kipfstuhl et al. (2009); Weikusat et al. (2009a,b) demonstrate that strain energy sufficient to drive strain-induced boundary migration occurs throughout ice sheets, including firn layers.This view is also supported by the numerical simulations of Steinbach et al. (2016).
R1: The combination of anisotropic textures and the high mechanical anisotropy of the ice monocrystal makes the ice polycrystal highly anisotropic.This has been confirmed by laboratory experiments that have shown that the strain-rates for secondary creep of polycrystalline ice depends on the fabric and its orientation with respect to the test configuration; while the tertiary creep is independent of the initial fabric and depends only on the test configuration.
Yes, this final point is precisely the basic tenet of our theory, but only because the stress configuration has been applied for long enough to achieve the tertiary state.Tertiary creep of polycrystalline ice is dependent upon the stress configuration and whether or not it is occurring is not simply a matter of the activity of migration recrystallisation.
R1: All the studies presented in Sec.2.1 as 'microstructure approaches', where motivated by representing the mechanical anisotropy of polycrystalline ice that develops as a result of the development of strain-induced preferred orientations of the ice crystals.They can not be opposed or compared to the flow law presented here, as they don't represent the same physical processes.
We disagree with the reviewer.What all flow relations aim to do is represent the bulk deformation of a polycrystalline aggregate.Microstructure based flow relations, where fabric is an input, can be applied equally to predictions of secondary or tertiary creep, or by extension to any arbitrary combination of stress and fabric.ESTAR only claims to represent tertiary creep -a subset of the cases that can be compared with microstructure based flow relations.Consequently comparisons with microstructure based approaches are entirely valid provided that the appropriate combinations of fabric and stress/deformations are considered.Laboratory experiments and ice cores for which corresponding strain rate data are available from repeat measurements of the borehole inclination are valuable sources of data that can be applied in such comparisons.This is particularly true for laboratory measurements of tertiary creep since the applied stresses, strain rates and microstructure data are all available.
R1: Laboratory experiments for tertiary creep are certainly relevant for the regions where migration recrystallization occurs; i.e. in temperate glaciers, in the bottom part of the central regions of the icesheets and certainly in the margins where temperature and stresses are sufficiently high.These areas are potentially important for the flow dynamics of polar ice sheets.The flow law presented in this paper must be presented in this perspective, and not as an empirical approach to represent the strain-induced mechanical anisotropy that results from the development of ice fabrics as seen in all deep ice-cores.This implies to largely rewrite large parts of the abstract, introduction and second section, to discuss the validity of the assumptions of the flow law and their domain of applicability for large-scale ice flow modelling.
Based on the preceding discussion, we do contend that ESTAR is a legitimate '...empirical approach to represent the strain-induced mechanical anisotropy that results from the development of ice fabric'.We have slightly modified the manuscript to more clearly highlight the 'domain of applicability' of ESTAR (P4L31-P5L18): "A second approach comprising experimental and observational approaches (Li et al., 1996;Wang et al., 2002a,b), modelling (Wang andWarner, 1998, 1999;Hulbe et al., 2003;Wang et al., 2003Wang et al., , 2004;;Breuer et al., 2006;Wang et al., 2012), and theoretical studies (Warner et al., 1999), has focussed on the development and assessment of an anisotropic flow relation for polycrystalline ice in which the nature of the crystal fabric and the magnitude of strain rate enhancement, E, are both regarded as determined by the stress regime.This assumption is supported by experimental observations for pure polycrystalline ice, which demonstrate that an accumulated strain of ∼ 10% is required for the microstructure to evolve to a state that is compatible with the flow configuration, irrespective of its initial condition (Jacka and Maccagnan, 1984;Gao and Jacka, 1987;Li and Jacka, 1998;Treverrow et al., 2012).Specifically, this approach regards the fabric and the enhancement in tertiary flow as determined by the relative proportions of the simple shear and normal deviatoric stresses.For such flow relations, it is typically assumed that the spatial variation in dynamic conditions (e.g., flow configuration and temperature) only occur gradually in an ice sheet, so that the microstructure evolves to maintain compatibility with these conditions.Through most of an ice sheet we expect that the rate of microstructural evolution generally exceeds the rate at which the flow configuration varies, and that the distances travelled by a parcel of ice during the time taken to develop a compatible fabric are typically small compared to the relevant ice sheet spatial scales.
The anisotropic flow relation proposed by Budd et al. (2013) represents a continuation of this strand.They found that a scalar anisotropic flow relation, i.e., one maintaining the collinear relationship between the components of ε and σ (τ e is a scalar function of the second invariant of σ ) provides a good fit to laboratory data from combined compression and shear experiments.Such a scalar anisotropic rheology also simplifies the requirements for implementation within ice sheet models that are already compatible with the (scalar) Glen rheological description.Budd et al. (2013) proposed what we term ESTAR as a suitable candidate scalar anisotropic rheology generalised to arbitrary stress configurations (i.e., not restricted in its application to the limited set of experimental stress configurations described in Li et al. (1996) and Budd et al. (2013))."Rewriting 'large parts of the abstract, introduction and second section' is not required.It is important to bear in mind that this is a paper about the implementation of a flow relation, not its formulation (P3L12-17).
Point 3: Application to the embayed ice shelf R1: I think it would have been more interesting to test the model with larger transverse dimensions, i.e. to represent non-embayed ice-shelves, as I assume that simple shear along the margin will be less important and the difference with Glen and a uniform enhancement factor for shear should increase?
As mentioned in the manuscript (P11L17-22) we have carried out the embayed ice shelf experiment for the cases where the transverse dimensions are L ∈ [20, 60, 100] km.Figures 1-3 below show the ESTAR to Glen velocity ratio, ESTAR to Glen thickness ratio, and λ S for each of the three cases.Note that the figures are plotted to approximately preserve their aspect ratio.It is clear from these figures that as the real aspect ratio increases, the ice shelves become flatter, so the proportion of the ice shelf that is shear dominated does not change markedly.Specifically, approximately 60% of the L = 20 km ice shelf is shear dominated (i.e., λ S > 0.5), 63% of the L = 60 km ice shelf is shear dominated, and 67% of the L = 100 km ice shelf is shear dominated.This is an important result, demonstrating that, contrary to the reviewer's expectations, changing the aspect ratio does not result in an embayed ice shelf that approaches the unembayed situation (i.e., with free-slip at the side walls).
Furthermore, regardless of aspect ratio the pattern and magnitude of the velocity ratios is similar in all cases, the maximum difference being approximately 15% near the ice-ocean front where normal stresses dominate.
Hence, given the similarities in the influence of the rheology on the dynamics between the aspect ratios, we prefer to focus our discussion in the manuscript on only one transverse dimension, L = 20 km.We have updated the manuscript to be more explicit about the consequences of increasing the aspect ratios, as per P11L17-22: "The experiment was carried out for model domains with transverse spans x ∈ [0, L], for L = 20, 60, and 100 km and along-flow dimension y ∈ [0, 100] km.The initial ice thickness decreases uniformly from 1000 m at the grounded zone to 300, 600, and 850 m at the ice front for the L = 20, 60, and 100 km cases, respectively.The main features of the anisotropic effects are similar regardless of aspect ratio.This is principally because wider embayed ice shelves are flatter so that the influence of simple shear stresses on the dynamics is not particularly sensitive to aspect ratio.Accordingly, we focus our discussion on one transverse length scale: L = 20 km." Point 4: ISMIP-HOM experiments R1: When discussing the results for Glen in the main text it is never clear for which value of the enhancement factor.The fact that the velocity scales with the enhancement factor is not a result, but as written in the manuscript follows from the definition of the experiment.So it would be more clear for the reader, to discuss the results for only one value of the enhancement factor or maybe better, find the value of the Glen enhancement factor that minimize the velocity difference with the tertiary creep flow law, and discuss this value and the results for this value.
The ISMIP-HOM experiments are scenarios in which the bed-parallel shear is the main driver of the flow.Hence, it emerges as appropriate that we choose a value for the Glen enhancement factor equal to the shear enhancement factor, i.e., E G = E S = 8.Accordingly, we have adopted the reviewer's first suggestion and have deleted the paragraph from the manuscript that compares the results for different values of E G .Since in response to Reviewer 2 we have introduced additional experiments (see details later), and as an interested observer would recognise the overall scaling properties, we have not considered it useful to take up the alternative suggestion of fine-tuning E G for each experiment.On P13L23-25 of the manuscript for ISMIPB, we clarify that our discussion of the results concentrates on the value E G = 8, which is the most appropriate overall choice for this experiment: "In what follows, we consider the case when the Glen enhancement factor is equal to the ESTAR shear enhancement factor, i.e., E G = E S = 8.This is the most relevant case for the ISMIPB experiment as the dynamics here are driven by bed-parallel shear, as discussed below." For ISMIPD, we include the following lines, echoing the above lines regarding ISMIPB that here we compare only the results from the case when E G = E S = 8 (P15L17-19): "Consistent with ISMIPB, the control of the final deformation flow in the ISMIPD experiment is bed-parallel shear, so we consider the case when the Glen enhancement factor is equal to the ESTAR shear enhancement factor, i.e., E G = E S = 8." R1: The discussion in page 12 about the value of the shear fraction along the transition curve is not very clear.
We have amended the discussion about the shear fraction along the transition curve on P14L23-P15L5 to be clearer, as follows: "In order to examine the dynamics giving rise to the high shear-dominance peaks in Fig. 6d and Fig. 7d, we consider the following exact form of λ 2 S (for these two-dimensional flow fields) expressible using the cartesian frame strain rate components for some spatially varying coefficients α, β, and γ.Since there is no surface accumulation, velocities and hence local non-rotating shear planes at the ice sheet surface are parallel to the surface.The traction free surface boundary condition implies that the numerator ( ε 2 ) in Eq. 23, and accordingly λ S , vanishes at the surface, except that if εe also vanishes, λ S is technically undefined.Our implementation sets λ S = 0 for vanishing ε in such situations.It is apparent from Eq. 23 that along the transition curves, i.e., where εxx = εzz = 0, λ 2 S = β, independent of (non-zero) εxz strain rate.One can show that β → (1 − S 2 x ) 2 towards the surface (i.e., for surface slope in the x-direction S x ) along the transition curve, in order to satisfy the surface boundary condition.This indicates that λ S would be finite along the transition curves all the way to the surface, except that we enforced its vanishing there.For these locations, the Glen and ESTAR viscosities corresponding to Eqs. 2-3 would tend to infinity as εe vanished approaching the surface, but are limited to a maximum value in the ISSM implementation.
Note that away from the transition curves λ S goes to zero as we approach the surface, associated with vanishing shear on the non-rotating shear plane and the corresponding dominance of normal deformations.We return to these near-surface spikes in λ S in the discussion."

REVIEWER 2
Point 1: Anisotropy R2: I have had a chance to look at the other reviewer comments and I have to say that I fully agree with him/her in the main two points of his/her review.First, the authors always refer to the method as "anisotropic", it is the "A" in the acronym.That is not simply misleading, it is wrong.In an isotropic rheology the relation between stress and strain-rate is a scalar number, twice the viscosity, in a more general anisotropic rheology it is a tensor.The method presented in the paper is simply not anisotropic.
The question of the anisotropic nature of ESTAR has already been addressed in our response to Reviewer 1. ESTAR is an anisotropic flow relation.The fact that ESTAR describes the rheology of anisotropic polycrystalline ice via a scalar enhancement factor does not preclude it from being an anisotropic flow relation.Reviewer 2 is incorrect in claiming that a scalar (collinear) relationship between the stress and strain rate tensor components is a necessary feature of an isotropic flow relation and that a tensor relationship between the stress and strain rate tensors is a necessary condition for an anisotropic rheological description.This issue has been addressed previously in the literature (Placidi et al., 2010).
The Continuum-mechanical, Anisotropic Flow model, based on an anisotropic Flow Enhancement factor (CAFFE) (Placidi et al., 2010) is an anisotropic flow relation that features an enhancement factor that is scalar function of a deformability parameter that is defined in terms the fabric and stress field.Faria (2008) provides a derivation of the anisotropy of CAFFE and Section 4 of (Placidi et al., 2010) provides further verification of its anisotropy, despite its collinear nature.Conceptually the same arguments apply to ESTAR, since it also features collinearity of the stress and strain rate tensors and a scalar enhancement that is an anisotropic function of the stress configuration.
As discussed in the response to Reviewer 1, the required ('tertiary') compatibility of the anisotropy with the stresses prevents any idea that stress can be rotated relative to the anisotropic material as the reviewer implies -this is the essence of anisotropy.
R2: Also in this point, the authors refer to Glen flow relation vs ESTAR throughout the manuscript.ESTAR as far as I can see is a Glen flow relation with a pre-exponential factor that varies spatially to account for ice fabric.A in Equation 1 is known to depend on ice fabric and a myriad of other things (not only on temperature as stated in the manuscript), the whole point of ESTAR is that it is giving a method to estimate the effect of ice fabric on A.
One could regard the ESTAR flow relation as it is presented, as a fabric-dependent enhancement factor for the Glen flow relation in tertiary flow.We regard this as essentially a semantic point since the task of implementing this factor is identical regardless of the viewpoint, and the enhancement function would still need a name that would encompass, empirical, tertiary and anisotropic.In terms of tertiary flow, it is not entirely certain (e.g.Treverrow et al., 2012) that tertiary flow does only differ from the Glen relation by this pre-factor.Regarding the "myriad of other things" known to influence A in addition to temperature, we contend that these are usually persistent material properties, like dissolved or particulate impurities, rather than dynamically evolving features such as crystal anisotropy.We now make this distinction more clearly in the manuscript, as per P2L18-20: "A(T ) is a flow parameter (Pa −n s −1 ), dependent on homologous temperature T and persistent material properties, for which various parameterisations exist based on laboratory tests and field measurements (e.g., Budd and Jacka, 1989;Cuffey and Paterson, 2010)." Point 2: Tertiary creep R2: And second, the empirical relation to extract the enhancement factor is based on laboratory experiments of tertiary flow, and it may well be that it can not reproduce the widespread observations of strain-induced anisotropy in polar ice.There is no discussion in the paper about how the proposed method explain observations or even expectations of the effect of ice fabric in polar flow.Do the vertical or horizontal variations of enhancement factor make sense according to observations of fabric or strain-rate fields?
The aim of this paper is to describe the numerical implementation of a flow relation or ice rheology within an ice sheet model (P3L12-17).While ESTAR is not directly concerned with determining straininduced anisotropy, there is a connection through the empirical observation that the flow-induced anisotropy is connected to the stress regime, so incorporation of ESTAR into an ice dynamics model is clearly an essential step in testing whether this can "reproduce the widespread observations of straininduced anisotropy in polar ice".
This flow relation, ESTAR, has been previously described (Budd et al., 2013) and its performance evaluated using ice core, borehole and laboratory data (Treverrow et al., 2015).It is outside the scope of this paper to present a comprehensive review of how ESTAR "explain(s) observations or even expectations of the effect of ice fabric in polar flow", or to provide an even more general review on the effects of polycrystalline anisotropy on polar ice sheet dynamics.The background to ESTAR is described in Sects.2.2-2.3.It is important to realise that implementation of ESTAR (or any other rheological description) in an ice sheet model provides an additional quantitative means by which that flow relation can be further tested and subsequently applied in simulating ice sheet dynamics.Regarding the validity of the assumption that tertiary creep is the predominant mechanism of deformation in the dynamically active region of polar ice sheets, we refer the reviewer to preceding comments in this response document.We also slightly modify the sections from the previous manuscript detailing conditions under which we expect ESTAR to be valid (P4L31-P5L18; quoted directly in the response to Reviewer 1 above), and provide new paragraphs discussing regions where it might not hold (P5L21-P6L12; again, quoted directly in the response to Reviewer 1 above).
In response to the last item above: R2: Do the vertical or horizontal variations of enhancement factor make sense according to observations of fabric or strain-rate fields?
We refer the reviewer to Treverrow et al. (2012) and Budd et al. (2013) and references therein for a detailed discussion on the variation in strain rate enhancement according to the stress configuration.Wang and Warner (1999), Wang et al. (2002a), andTreverrow et al. (2015) each discuss the vertical variations in enhancement that are required to explain observed strain rate profiles within an ice sheet and also present fabric profile information.
Point 3: ESTAR vs tuning hardness parameter R2: I also miss an important point in the paper, what is the point of using ESTAR in a large-scale ice sheet model?
We have interpreted this question as 'what is the point of using any anisotropic flow relation in any large-scale ice sheet model?'.
There has long been recognition that a flow relation that realistically incorporates the anisotropic rheology of polycrystalline ice is a necessary component of any ice sheet model.While a diverse range of anisotropic flow relations have been proposed, implementing such a relation within an ice sheet model in computationally efficient manner has proven challenging.Additionally, there has been a long standing awareness that using a Glen-type flow relation with a spatially-invariant enhancement factor does not provide an adequate description of polycrystalline ice rheology.What we are presenting is our attempt at incorporating a description of anisotropic rheology into an ice sheet model (P3L12-17).
R2: I may be missing something here but large-scale models are capable to initialize the local depth-averaged hardness parameter in the Glen flow law using satellite observations.What is ESTAR adding?
While it is possible to "tune" a "local depth-averaged hardness parameter" as well as a basal friction parameter using inverse methods to match modelled and observed surface velocities, such an approach provides no physical insight into the local controls on deformation rates and no indication how the "depth-averaged hardness" should change as the ice sheet evolves (note, we also do not want to invert for both ice rheology parameters and basal friction, as there would be an infinite number of solutions -they have "similar" effects on the surface velocities.This is why we only invert for B on floating ice, and friction on grounded ice, but never both at the same place).ESTAR is adding the ability to incorporate the spatial variability of large-scale anisotropy within ice sheets into simulations of ice sheet evolution, with a physical basis, in a fashion that will remain valid as the system evolves.It would be possible to explore using inverse methods to find appropriate global values for E C and E S , since ESTAR describes the spatial variations, encoded in λ S , e.g.(P19L32-34): "Indeed, with the implementation of ESTAR in ISSM, it might be possible to use inverse methods to search for values of E S and E C that improved the match between modelled and observed surface velocities."However, there is clearly too little data in a 2D surface velocity field to determine a 3D varying "hardness parameter" throughout an ice sheet.
R2: Could ESTAR inform the initialization of a large-scale model somehow?It could be applied to simulations were there is no satellite observations but then, as I suggested in my previous point, shouldn't we check how well does it do in a case where we have observations?Also, are there any significant vertical variations on enhancement factors that are captured by ESTAR but can't be inverted by a large-scale model?In any case, the paper needs be clearer about what is ESTAR and what can it do.
The purpose of implementing and testing ESTAR in a large-scale ice sheet model is to explore what advantages it offers.We have shown that it is a computationally-efficient, physically based description of anisotropic flow (Sects.5-7).We are reluctant to further expand Section 2.3 about the theory underlying ESTAR, as this paper is about the implementation of a flow relation that has already been discussed at length by Budd et al. (2013).In this manuscript, we have explained the implementation and explored what ESTAR can do in some preliminary test examples; its implementation is a necessary precursor to testing it in other more realistic situations, including the reviewer's suggestion to "check how well it does".Indeed, we are currently applying ESTAR to an Antarctic glacier system where observations are available to validate our simulations and these are forthcoming.However, it is outside the scope of this current manuscript to include full simulations of real glaciers.
As already discussed, large-scale models are capable of initialising the "local depth-averaged hardness parameter" in the Glen flow relation (relying, of course, on also tuning of the basal traction parameter over grounded ice) using satellite observations.However, this sort of initialisation, without physical insights into the processes controlling the tuned parameter, is problematic.As the system evolves, the various physical processes that a "black box" parameter represents may change, rendering the "initialisation" values inappropriate.ESTAR could potentially inform the initialisation of an ice sheet model, through determining appropriate global values of enhancement factors E S and E C , using the physics embodied in ESTAR through λ S to include spatial variations in 3D (P19L32-34).Indeed, in our ice shelf example, an essentially 2D flow problem in the absence of temperature effects, one could attempt to invert for the values of E S and E C to see if this improved agreement between modelled and observed velocities.However, the enhancement factor in ESTAR is not physically or dynamically equivalent to the overall flow parameter in the Glen flow relation.Furthermore, in this study, we have used the full-Stokes 3D implementation of ESTAR, which relies on more than just depth-averaged parameters.Indeed, there is insufficient data in satellite-derived surface velocities to derive a local, 3D ice flow parameter.
R2: Also, are there any significant vertical variations on enhancement factors that are captured by ESTAR but can't be inverted by a large-scale model?In any case, the paper needs be clearer about what is ESTAR and what can it do.Wang et al. (2002a) demonstrate the need for a significant vertical variation in enhancement factor -simply derived from observations -together with a convincing connection to both an ESTAR-type prescription that yields that depth variation, and to the anisotropic crystal orientation fabrics observed.Treverrow et al. (2015) continues this study by exploring how ESTAR compares with other candidate flow relations.Both these studies advocate the need for a depth dependent enhancement factor, though in the general context of an assumed driving stress.This further drives the motivation for implementing ESTAR in an ice dynamics model that can solve for the full stress field to extend the testing of ESTAR against observations.Point 4: Aspect ratios R2: Finally, I find really intriguing that the authors state that results are very similar with different aspect ratios in the experiments presented in Sections 5.1 (P9-L31) and 5.2 (P11-L25).What do they mean?In the embayed ice-shelf, I would expect the ice fabric induced by extension and shear at the margins to be different.The narrower the ice-shelf is I would expect that the overall effect of lateral shearing should be more important.In the ice flow over a bumpy bed, L is the wavelength of the bedrock undulations.I would expect the increase in basal roughness to have a strong control over the orientation of the fabric.What does it mean that aspect ratios don't affect ESTAR results?I would like to see some discussion about the results presented.
For the embayed ice shelf, as discussed above in response to reviewer 1, we have focussed our discussion in the manuscript on the case when the transverse dimension is L = 20 km.Figures 1-3 in this document show the Glen and ESTAR steady-state surface velocity ratio, the Glen to ESTAR thickness ratio, and λ S , for the cases when the transverse dimensions are L ∈ [20, 60, 100] km, respectively.The key point highlighted in these figures is that as the aspect ratio increases, the ice shelves become flatter, and the overall proportion of the ice shelf that is shear dominated does not change markedly (so the embayed ice shelf does not approach the unembayed situation).Hence, and in the interests of constraining the length of the manuscript, we have omitted the results of the cases when L = 60 and 100 km.We have amended the manuscript to be clearer on the impact of increasing aspect ratio, as highlighted on P11L17-22: "The experiment was carried out for model domains with transverse spans x ∈ [0, L], for L = 20, 60, and 100 km and along-flow dimension y ∈ [0, 100] km.The initial ice thickness decreases uniformly from 1000 m at the grounded zone to 300, 600, and 850 m at the ice front for the L = 20, 60, and 100 km cases, respectively.The main features of the anisotropic effects are similar regardless of aspect ratio.This is principally because wider embayed ice shelves are flatter so that the influence of simple shear stresses on the dynamics is not particularly sensitive to aspect ratio.Accordingly, we focus our discussion on one transverse length scale: L = 20 km." Regarding the ISMIP-HOM experiments (Sects.5.2 and 5.3), we wish to thank the reviewer for encouraging us to explore more deeply the cases of bed variations with shorter wavelengths.We did observe greater variations in the anisotropic case, and we have also developed a more extensive discussion of the results (Sect.6).The analysis also provides some indications of regions where the assumptions of the ESTAR rheology may not hold and this is now addressed in the discussion (Sect.6), though the consequences of that type of lapse in realistic ice sheets requires further study.
Given that ESTAR does not directly calculate crystal orientations our interpretations do not directly address the reviewer's interest in fabrics, though we do provide some interpretations in the discussion (P17L31-P19L22): "Our idealised test cases also provide some insights into the validity of the tertiary flow assumption underlying ESTAR, and the development of anisotropic crystal fabrics compatible with the current deformation regime.In the embayed ice shelf test the most significant change in the deformation regime is clearly the transition to extensive flow on approach to the ice shelf front.The contours of λ S here are relatively well aligned with the ice flow so that flowing ice experiences gradual changes in stress regime, and the magnitudes of strain rates (e.g., Fig. 4b) and velocities (Fig. 3) indicate that in the region near the ice front, where the Glen and ESTAR results vary appreciably, a progression of essentially compatible fabrics would be maintained.Indeed, under the prevailing deformation and flow conditions these would even develop from random fabrics over a few km.
The ISMIP-HOM experiments reveal potential violations of the tertiary flow assumption, although the significance for the flow field of these apparent short-comings needs to be assessed with regard to the somewhat artificial nature of the tests.Indeed as we saw, the difference between the ESTAR predictions and the Glen flow relations (provided E G = E S ) was small, except for ISMIPB with L = 5 km, although of course ESTAR makes no claim to correctly describe the transient rheology of ice with an evolving anisotropy.
The ISMIP experiments have a spatial periodicity, which could allow one portion of the repetitive basal conditions to dominate the overall flow.Also, there is no surface mass budget in these experiments so that, as remarked earlier, the ice surface is a flow line, whereas in a system with surface accumulation fresh snow is always being added and advected down into the ice sheet where it makes the transition to solid ice.Accordingly, in the flow regime of these prognostic experiments even the surface layers would be regarded as having developed some anisotropy just as the lower layers would, since they have in principle been deforming over an arbitrarily long time.
The main issue about the establishment of tertiary flow conditions in the periodic environment of our ISMIP experiments concerns the possible cycling of the flowing ice through a variety of stress regimes.This leads to transition regions where the stress regime and presumably the crystal anisotropy would be evolving, and the compatibility assumptions behind ESTAR could locally be violated.Clearly the spatial extent of transitional flow and the delay in attaining any new tertiary state depends on the magnitudes of the strain rates and the velocity of the ice.By combining these with a threshold for accumulated strain as the criterion for development of a compatible (tertiary) fabric under a persistent flow regime, the extent of a transition zone can be estimated.This scale can then be compared to the horizontal variation of the stress regime.Selecting the 10% strain required to develop a compatible anisotropy from initially randomly oriented ice should provide a conservative yardstick, when applied to gradual changes in stress regime.
The patterns of stress regimes revealed by the distributions of λ S (Figs. 6d, 7d, 9a, and 10f) indicate where along-flow variations in stress regime might be too rapid to sustain the assumption that a compatible crystallographic anisotropy had evolved.For the ISMIPB experiments this concern is essentially focussed to the near surface spikes in λ S around the two locations where longitudinal deformations vanish, since the anisotropy of deeper ice will be compatible with deformation dominated by simple shear.There may be some complications with a slow cycling of the upper levels between extensive and compressive flow.Very near the surface, the λ S peak intervals are narrow and the shear strain rates there are very small (corresponding to transition scales of several kilometres) so that there will be no appreciable development of a shear compatible fabric.Either side of the peaks, the main extensive and compressive flow domains for L = 20 km (see Fig. 6e) are ≈ 5 km long and have transition scales of < 1 km which suggests that the strongly normal stress dominated upper layer will be mainly in tertiary state.Turning to the transient shear intrusions into this layer: at 100 m depth over the bump the shear transition scale is 3 km while at 200 m depth over the depression the shear transition scale is ≈ 5 km, suggesting that that the λ S peaks do represent a local failure of ESTAR's tertiary assumption.Throughout the domain the lower half of the ice column has transition scales of ≤ 300 m which, given the gradual variations in λ S and the direction of ice flow, indicates that region is in the tertiary state.
For ISMIPB with L = 5 km, which displays a generally deeper band of normal stress dominated regime (Fig. 7d), the transition scales for the compressive and extensive regions are ∼ 500 m for regions of ≈ 1 km in extent, while the shear transition scales at 100 m depth above the bump and 200 m depth above the depression are now ∼ 10 km and ∼ 1 km, respectively.For most of the domain the transition rates in the lower half of the ice column are ≤ 100 m, although this rises to nearly 1000 m above the bedrock bump.
In the ISMIPD case, for L = 20 km, the pattern of λ S (Fig. 9a) shows there are also transitions between simple shear dominated and extension dominated deformation associated with the slippery region, with varying abruptness at different depths, and some of the contours of λ S in this instance almost orthogonal to the ice flow.A complication is that there is very little longitudinal deformation (Fig. 9e) occurring over the slippery region because the overall flow is controlled by the periodic sticky spot.Accordingly, there would not be any significant fabric evolution across this ∼ 4 km region (estimated transition scales there are > 40 km) so that the tertiary assumption and using E C (since λ S = 1) would be inappropriate.Once again, the low strain rates here (Figs.9e-f) translate into very stiff ice and might make the influence of ESTAR enhancement factors relatively unimportant.A factor of 100 in εe changes viscosity by a factor of 21.5, whereas the maximum viscosity contrast from E(λ S ) is 1.39.The shear strain rates are also very low, with transition scales > 1 km except very close to the bed over the sticky spot.Accordingly, while a compatible fabric could be expected where the large λ S values are shown in Fig. 6a, its presence would be due to the periodic flow, and the inability of the λ S ∼ 0 region to modify it.
For the last test, ISMIPD with L = 5 km, strain-rates are once again very low, and there is no simple structure to the picture of the stress regime portrayed by λ S in Fig. 10f.Below mid-depth there is a periodically continuous band of shear that might favour the development of crystal anisotropy, but clearly the tertiary flow assumption of ESTAR would not be particularly useful here." For the ISMIP-HOM experiments with larger horizontal extents (i.e., L > 20 km), the extent to which the longitudinal stresses impact flow and the differences between ESTAR and Glen results both decrease as L increases, such that the dominant sensitivity of the flow to the bed aspect ratio is similar for both rheologies.This is generally consistent with the results concerning the longitudinal stresses presented by Pattyn et al. (2008), allowing for the fact that our simulations have evolved to steadystate whereas those tests were to find flows satisfying initial stress balance.By contrast, for smaller horizontal extents (e.g., L = 5 km), the increasing importance of longitudinal stresses leads to greater differences between Glen and ESTAR experiments.We overview the expanded ISMIPB and ISMIPD results in turn.
In ISMIPB the more rapid bed variation shows greater reductions in ESTAR velocities (stiffer ice) and higher spatial contrast in the velocity ratios.This is evident when comparing the L = 5 km results in new Fig. 4 below (Fig. 7 from the updated manuscript) with the case when L = 20 km in Fig. 6 from the updated manuscript.The smaller aspect ratio (L = 5 km) leads to a difference of up to 25% between the Glen and ESTAR velocities in the topographic depression, and a difference of almost 20% in the surface layers over the topographic bump.We have amended the manuscript in section 5.2 for ISMIPB to include discussion on this smaller aspect ratio, as follows (P14L14-22): "In addition to the case where L = 20 km, we also investigated the impact of reducing the horizontal extent to L = 5 km.In this steeper bed scenario (Fig. 7), the ESTAR surface velocities are at least 11% slower than the Glen velocities in the surface layers across the whole domain, as much as 20% slower around the topographic bump, and up to 25% slower in the topographic depression (Fig. 7c).The much greater reductions in the magnitude of the ESTAR velocities for L = 5 km are a consequence of the increasing importance of longitudinal stresses in the stress balance equations for the smaller aspect ratio (Fig. 7e), and also in some areas the lower strain rates, which lead to correspondingly stiffer ice.Indeed we see a clear decline in the shear strain rate in the lower part of the bed depression in Fig. 7f, in contrast to Fig. 6f.The qualitative pattern of the longitudinal strain-rates in Fig. 7e is similar to the L = 20 km case, although the horizontal gradients are naturally accentuated, and the "transition curves" are displaced." Note that while the contrast due to the ESTAR incorporation of anisotropy is greater for the L = 5 km case, the structure and patterns of the results are broadly similar, e.g.compare Figs. 6 and 7 from the updated manuscript.We have extended the discussion section (as suggested) to cover these additional results (P17L20-P17L30): "For more rapidly varying bed topography in ISMIPB, with L = 5 km, the differences in velocity for the two flow relations reached 25%, with surface variations of 11%...These results suggest if major bed topography only varied on scales much longer than the ice thickness, close agreement between ESTAR and the Glen flow relation might be achieved more generally by choosing the tertiary shear enhancement factor as the Glen enhancement factor (E G = E S ).This might provide a physical rationale to replace the ad hoc enhancement factors typically used in largescale grounded ice sheet modelling with the value appropriate to flow dominated by simple shear.However, larger differences between velocities and vertical shear profiles emerged for the more rapid bedrock variation, where the importance of including longitudinal stresses in the momentum balance is already recognised (Pattyn et al., 2008), suggesting that adopting ESTAR would be preferable." We do not observe similar dramatic changes in flow for ISMIPD when the horizontal extent is reduced (Fig. 5 below).In the ISMIPD experiments, sliding dominates over deformation to an extreme extent for both length scales.For L = 20 km (new Fig. 8 of the manuscript) the only difference between the velocities for the Glen and ESTAR cases are observed over the sticky spot.The differences between the horizontal velocities in ISMIPD tests are less than 1% in both cases, which already had prompted us to examine the ratio of the deformational contributions to the velocities as well.While these are clearly very small contributors to the flow, we observe that there is both a marked spatial contrast in the ratio, and a considerable change in that contrast between the two spatial scales.
For L = 20 km the aspect ratio appears to favour a simple correlation between the variation in bed friction coefficient and the deformational regime.Over the slippery spot the stiffer ice of the ESTAR case correlates with stiffer ice a reduced horizontal deformational velocity, but for L = 5 km a different, more complex picture emerges.The structure of the longitudinal strain rates is similar for both aspect ratios, but the pattern of shear strain rates is more complex for L = 5 km.We note the resulting interesting pattern in λ S (Fig. 5 below) compared with λ S for the L = 20 km case (Fig. 9 from the manuscript).This is commented on in the explanation of the surface spike structures (see below).Again, and as remarked in the updated manuscript (P16L4-13), the only region for which λ S is relevant and plays a role in the dynamics is the region over the sticky spot, where deformation is the dominant control of flow near the bed.In this region, the proportion of shear-dominated deformational flow does not increase markedly over the L = 20 km case.We have included discussion on the impact of the smaller aspect ratio (P16L14-29): "The results of a FS prognostic run to steady-state for the ISMIPD experiment for L = 5 km are presented in Fig. 10.The v x velocity ratio (Fig. 10c) shows that very little difference is seen between results for the two flow relations, unlike the ISMIPB experiments in the previous section, where differences up to 25% between ESTAR and Glen cases emerged for the shorter bedrock periodicity.The tiny differences in overall velocities are enhanced but the patterns in Fig. 9c and Fig. 10c are similar.However, there is a significantly different picture in the ratio of deformation velocities seen in Fig. 10d, compared to Fig. 9e.The largest differences are now limited to the lower portion of the ice column and for much of the region over the slippery bed the ratio is almost unity.The pattern of deformation regimes mapped by λ S (Fig. 10f) is also more complex than that seen in the preceding L = 20 km case (see Fig. 9d).The general structure of the normal strain rates is similar to previous experiments, but here the persistence of a band of shear (Fig. 10h) above the slippery spot at intermediate depths prevents the establishment of a vertical block of flow dominated by normal stresses.The shear profile above the sticky spot is also much weaker in the upper layers.Accordingly, λ S reveals a band of unevenly shear-dominated deformation which is continuous across the periodic domain.Once again, shear dominated spikes extend towards the surface in association with the vanishing of the normal strain rates.
The spatial variations in the viscosity ratio (Fig. 10e) depart significantly from those of λ S , reflecting more strikingly than for L = 20 km (Fig. 9e) the combined influence of the pattern of enhancement (controlled by E(λ S )) and the effect of different strain rates, with values both above and below the range (1.0-1.39)directly controlled by E(λ S )/E G ."

Other modifications
As indicated at the start of this response document, we have considerably expanded the discussion section as suggested (Sect.6).The additional observations led to some alterations in ordering in that section.We also introduced additional diagnostics in the results section (Sect.5) -including ratios of total and deformational velocities and of viscosities to assist our interpretation of the influence of ESTAR versus Glen rheology.We changed the presentation of shear strain rates to use a logarithmic colour scale to better capture the contrast and assist in the interpretation.The conclusion (Sect.7) has also been expanded to cover the additional insights and to respond to reviewers' concerns.