Introduction
Antarctic ice shelves are largely composed of ‘meteoric’ ice, which has primarily formed from compacted snow on the surface of the ice sheet that has subsequently flowed onto the ocean. However, many ice shelves also contain a layer of ‘marine’ ice, consolidated from accreted platelets which have formed in a mixture of salty and fresh water frozen in situ beneath the ice shelf (Lewis and Perkin, Reference Lewis and Perkin1986; Jenkins and Doake, Reference Jenkins and Doake1991; Jenkins and Bombosch, Reference Jenkins and Bombosch1995; Bombosch and Jenkins, Reference Bombosch and Jenkins1995). This marine ice layer commonly comprises between a sixth and a third of the total ice shelf thickness (Oerter and others, Reference Oerter1992; Eicken and others, Reference Eicken, Oerter, Miller, Graf and Kipfstuhl1994; Treverrow and others, Reference Treverrow, Warner, Budd and Craven2010), with its thickness varying spatially across the base of the ice shelf (Fricker and others, Reference Fricker, Popov, Allison and Young2001). Due to differences in its structure, chemistry and in situ temperature, marine ice is likely to have different rheological properties to meteoric ice (Moore and others, Reference Moore, Reid and Kipfstuhl1994; Treverrow and others, Reference Treverrow, Warner, Budd and Craven2010; Dierckx and others, Reference Dierckx, Peternell, Schroeder and Tison2014).
Compared to an ice shelf of the same geometry without a marine ice layer, marine ice can accelerate ice shelf expansion and thinning (e.g. Hulbe and others, Reference Hulbe, Johnston, Joughin and Scambos2005), cause stresses to localise in the upper part of the shelf (e.g. Lange and MacAyeal, Reference Lange and MacAyeal1986) and affect rifting behaviour (e.g. Kulessa and others, Reference Kulessa, Jansen, Luckman, King and Sammonds2014). In turn, changes in the ice shelf properties can have far-reaching consequences for buttressed ice upstream (Fürst and others, Reference Fürst, Goelzer and Huybrechts2015; Pegler, Reference Pegler2018; Reese and others, Reference Reese, Albrecht, Mengel, Asay-Davis and Winkelmann2018). However, the sensitivity of ice shelf models to the presence of marine ice, through its impact on the thermal profile and ice physical properties, has not been systematically studied. Here, we address this knowledge gap using a model of an idealised ice shelf.
Rheology in ice-sheet models
The relationship between applied stress and resultant strain rate in ice shelf models is defined using a material constitutive relation, also known as a ‘flow relation’. The most commonly used flow relation for ice shelf modelling applications is given by the following equation (Glen, Reference Glen1952, Reference Glen1955, Reference Glen1958):
where $\dot {\epsilon }$
is the strain rate tensor (s−1), A(T′) is an empirically derived parameter dependent on temperature relative to the pressure melting point, τe is the effective stress (Pa), σ′ is the deviatoric stress tensor and n is the stress exponent, which is generally assumed to be equal to 3, although this has been challenged by some (e.g. Treverrow and others, Reference Treverrow, Budd, Jacka and Warner2012; Qi and others, Reference Qi, Goldsby and Prior2017; Bons and others, Reference Bons2018; Millstein and others, Reference Millstein, Minchew, Pegler and Sciences2022). In this form, the law assumes that ice is mechanically isotropic. More recent studies (e.g. Treverrow and others, Reference Treverrow, Budd, Jacka and Warner2012; Budd and others, Reference Budd, Warner, Jacka, Li and Treverrow2013) have considered the addition of an enhancement factor, E, to Eqn (2), as follows:
where E is defined as the ratio of tertiary to secondary creep rates
This is a convenient way to parameterise the effect of microstructure or chemistry on enhancing strain rates during tertiary creep.
However, a constant enhancement factor is limited in its ability to accurately represent steady-state strain rates in a mass of ice where the nature of the applied stresses, particularly the proportion of simple shear to compression, varies spatially. Budd and others (Reference Budd, Warner, Jacka, Li and Treverrow2013) more recently proposed the Empirical Scalar Tertiary Anisotropy Regime (ESTAR) flow relation, as a constitutive flow relation to account for the impact of different stress configurations on the overall flow regime. This flow relation captures the differing effects of simple shear and compression via an enhancement factor E(λS) and is given by Graham and others (Reference Graham, Morlighem, Warner and Treverrow2018):
where E(λS) is defined as:
E C and E S are experimentally derived dimensionless enhancement factors obtained from the creep deformation of initially isotropic ice, in either compression or simple shear respectively. Specifically they are obtained from the ratio of strain rates measured at the tertiary and secondary (minimum) stages of creep in each case. Treverrow and others (Reference Treverrow, Budd, Jacka and Warner2012) found that the ratio E S/E C is generally equal to 8/3, irrespective of stress configuration.
λS is the shear fraction – the ratio of the magnitude of shear on the non-rotating shear plane ||τ′|| to the effective stress – defined as:
The shear fraction λS takes values in [0, 1]: under compression alone scenarios, λS = 0 and E(λS) = E C; under simple shear alone scenarios, λS = 1 and E(λS) = E S. λS and the physical interpretation of the non-rotating shear plane are discussed in depth by McCormack and others (Reference McCormack2022).
It should be noted that the ESTAR flow relation assumes that the flow of most ice sheets and glaciers is steady-state (tertiary) creep, meaning the ice sheet has undergone sufficient deformation that the microstructure is compatible with the underlying stress configuration; i.e. the crystallographic preferred orientation (CPO) and grain size of the ice are what is expected to occur in meteoric ice in the tertiary (steady-state) creep stage. A detailed discussion of the underlying assumptions of ESTAR and its limitations can be found in Graham and others (Reference Graham, Morlighem, Warner and Treverrow2018) and McCormack and others (Reference McCormack2022).
Characteristics of marine ice
Our knowledge of the rheological behaviour of ice is largely informed by ice deformation experiments. These are typically performed on either laboratory-made ice with no pre-existing CPO or chemical impurities (e.g. Kamb, Reference Kamb1972; Durham and others, Reference Durham, Stern and Kirby2001; Goldsby and Kohlstedt, Reference Goldsby and Kohlstedt2001; Qi and others, Reference Qi, Goldsby and Prior2017), or on meteoric ice (e.g. Gao and Jacka, Reference Gao and Jacka1987; Dahl-Jensen and others, Reference Dahl-Jensen, Thorsteinsson, Alley and Shoji1997; Treverrow and others, Reference Treverrow, Budd, Jacka and Warner2012; Craw and others, Reference Craw, Qi, Prior, Goldsby and Kim2018). Few deformation experiments have been performed on marine ice, and so our knowledge of its rheological properties, and how they differ from those of meteoric or laboratory-made ice, is limited.
Given the scarcity of deformation experiments performed on marine ice, in order to understand how its rheology might differ from the overlying meteoric ice, we must consider variations in its temperature, chemistry and microstructure.
Temperature
Marine ice forms at the in situ freezing point of the ocean water in an ice shelf cavity (typically ≈ −2°C). Often this marine ice layer remains isothermal, steepening the thermal gradient of the ice above and raising the bulk temperature of the ice shelf. This has been observed repeatedly in borehole temperature profiles of the Amery Ice Shelf (hereafter Amery IS), as shown in Fig. 1a (Craven and others, Reference Craven, Allison, Fricker and Warner2009; Wang and others, Reference Wang, Zhao, Gladstone, Galton-Fenzi and Warner2022). As the strain rate of deforming ice in response to an applied stress is highly dependent on temperature, a layer of warm ice at an ice shelf base will lead to enhanced flow in the ice shelf and higher local stresses in the shallower meteoric ice. In the case of the Amery IS, Craven and others (Reference Craven, Allison, Fricker and Warner2009) estimate (using the ice deformation temperature-dependence relationship outlined by Budd and Jacka, Reference Budd and Jacka1989) that ice in the marine ice layer may be 15% less viscous as a direct result of this temperature perturbation. However, a temperature perturbation is not always present. Eicken and others (Reference Eicken, Oerter, Miller, Graf and Kipfstuhl1994) identified a ~ 70 m-thick layer of marine ice (comprising ~ 30% of the total shelf thickness) at the base of the Filchner-Ronne Ice Shelf (hereafter Filchner-Ronne IS) through electrical conductivity measurements, but the corresponding borehole temperature profile was close to linear, as would be expected in an ice shelf without a marine ice layer (Fig. 1b). This is likely due to a slower rate of marine ice accretion at that location, which is not the case throughout the entire shelf. Grosfeld and Thyssen (Reference Grosfeld and Thyssen1994) used a 2-D thermal model, in combination with a partial temperature profile collected at a different location, to investigate how the shape of thermal profiles may vary throughout the ice shelf. They concluded that an ‘S-shaped’ thermal profile associated with a isothermal basal layer, as seen in the Amery IS, likely forms where basal accretion rates are high (up to 2 m a−1 in the Filchner-Ronne IS), and the ice shelf is not able to thermally equilibrate as large amounts of warm ice are added to the base. Where the accretion rate is lower (closer to zero or even lower as melting dominates), the thermal profile tends towards the other extreme, a ‘cubic’ profile where only the ice at the base is close to ocean temperature, and the ice above has equilibrated to give a linear temperature gradient across the marine–meteoric ice transition. A range of thermal profiles between these two extremes may exist throughout an ice shelf, depending on the local marine ice accretion rate.
Ice shelf temperature profiles: (a) from borehole measurements on the Amery IS, with an isothermal basal marine ice layer (Craven and others, Reference Craven, Allison, Fricker and Warner2009; Wang and others, Reference Wang, Zhao, Gladstone, Galton-Fenzi and Warner2022); (b) from borehole measurements in the Filchner-Ronne IS (Eicken and others, Reference Eicken, Oerter, Miller, Graf and Kipfstuhl1994), where the basal layer is not isothermal; (c) an example of those imposed during the simulations in this paper. Depth values are normalised as a proportion of the total shelf thickness. The approximate thickness of the marine ice layer is shaded in grey.
Porosity and tortuosity are also likely to affect the distribution of heat through the ice shelf. Grosfeld and Thyssen (Reference Grosfeld and Thyssen1994) suggest that a ‘slushy’ layer of higher porosity has an insulating effect in the marine ice, preventing thermal equilibration and allowing the influx of heat from direct contact with sea water. However, ice with a high level of tortuosity and interconnectivity in its pore spaces will experience a greater influx of heat from the ocean. This is a likely contributor to the isothermal basal layer seen in the Amery IS, as the lower part of the marine ice layer, comprising nearly half of its volume, is hydraulically connected (Craven and others, Reference Craven, Allison, Fricker and Warner2009).
It has been shown in modelling of the Ross Ice Shelf (hereafter Ross IS) that spatial variations in basal melting and freezing rates have a significant impact on the effective viscosity of the ice (Rommelaere and MacAyeal, Reference Rommelaere and MacAyeal1997), suggesting that a better understanding of these processes is important to model ice shelf behaviour effectively.
Chemical composition
Marine ice has a salinity intermediate between those commonly measured in sea ice and meteoric ice (typically <0.15% (Eicken and others, Reference Eicken, Oerter, Miller, Graf and Kipfstuhl1994; Khazendar and others, Reference Khazendar, Tison, Stenni, Dini and Bondesan2001)), as it is formed from a mixture of sea water and fresh water, and brine is rejected from the ice crystals as they form. It is generally accepted that the presence of chemical impurities such as dissolved salt in polycrystalline ice alters its rheological behaviour (Stoll and others, Reference Stoll, Eichler, Hörhold, Shigeyama and Weikusat2021), however the precise magnitude and nature of this effect is not well understood.
Trace amounts of H2SO4 and NH3, when added at concentrations commonly seen in meteoric ice, have been observed to weaken synthetic ice deformed in the laboratory (Li and others, Reference Li, Iliescu and Baker2009; Hammonds and Baker, Reference Hammonds and Baker2018), whereas a similar concentration of Ca2+ has been observed to strengthen it (Hammonds and Baker, Reference Hammonds and Baker2016). It is unclear how these effects might scale to the larger concentrations likely to be present in marine ice. Observations of ice cores in Greenland suggest that ice with higher impurity content is generally softer (Dahl-Jensen and Gundestrup, Reference Dahl-Jensen and Gundestrup1987), but where the impurities have other controls on microstructural evolution, for example by shifting the balance of deformation mechanisms to result in a stronger CPO as observed by Dahl-Jensen and others (Reference Dahl-Jensen, Thorsteinsson, Alley and Shoji1997), they can in fact make ice harder to deform in some stress configurations.
It is often assumed that marine ice has a weakening effect in studies of ice shelves where marine ice exists as a heterogeneous mélange inside a suture zone or infilled rifts (e.g. Khazendar and others, Reference Khazendar, Rignot and Larour2009; Kulessa and others, Reference Kulessa, Jansen, Luckman, King and Sammonds2014; King and others, Reference King, Rydt and Gudmundsson2018). However, in these scenarios it is difficult to separate any effects of chemistry from those of temperature, and rheological heterogeneity within the mélange. In fact, Dierckx and Tison (Reference Dierckx and Tison2013) performed deformation experiments on marine ice from the Nansen Ice Shelf, East Antarctica, and observed no obvious weakening as a result of salt content, compared with impurity-free samples from other studies deformed at the same temperature. It is important to note that those experiments were stopped at the secondary creep stage, the point at which the deformation rate in response to a constant applied stress has reached a minimum. This minimum can be understood as the point where strain-hardening processes that lead to a decreasing strain rate during primary creep are balanced by the onset of deformation and recovery processes during the transition to tertiary creep (Budd and Jacka, Reference Budd and Jacka1989). It is possible that chemical impurities may have a greater effect on rheology during tertiary creep (after a constant deformation rate has been established), when dynamic recrystallisation mechanisms are more active. Craw (Reference Craw2018) found that while impure meteoric ice was slightly weaker than pure laboratory-made ice at the secondary creep stage, it was even weaker in tertiary creep. This may be due to effects such as impurity drag and heightened dislocation density affecting grain size and CPO strength more as deformation progresses (Alley and others, Reference Alley, Perepezko and Bentley1986; Obbard and Baker, Reference Obbard and Baker2007; Eichler and others, Reference Eichler2017).
Overall there is no clear consensus on the effect that chemical impurities have on deformation, but in this study going forward we will assume that they cause a small amount of mechanical weakening during the later stages of creep, consistent with the findings of Craw (Reference Craw2018).
Microstructural characteristics
Due to differences in formation mechanisms, in situ temperature and possibly chemical content, marine ice has a distinctly different microstructure to meteoric ice found in the same geographic location. Most relevant here are differences in grain size and CPO, both of which can have a significant effect on rheological behaviour (Budd and Jacka, Reference Budd and Jacka1989).
Meteoric ice is typically composed of crystals on the scale of centimetres, whereas marine ice crystals are smaller by comparison, on the scale of millimetres, sometimes increasing with depth (Moore and others, Reference Moore, Reid and Kipfstuhl1994; Treverrow and others, Reference Treverrow, Warner, Budd and Craven2010). This difference may be a function of chemical content, as higher concentrations of impurities in natural ice are associated with smaller grain sizes (Dahl-Jensen and others, Reference Dahl-Jensen, Thorsteinsson, Alley and Shoji1997; Obbard and Baker, Reference Obbard and Baker2007). Bubble content can also limit the rate of grain growth in ice by inhibiting grain boundary migration (Azuma and others, Reference Azuma, Miyakoshi, Yokoyama and Takata2012). Finer grain sizes in ice generally lead to rheological weakness, due to an increase in grain-size sensitive creep mechanisms (Paterson, Reference Paterson1991; Goldsby and Kohlstedt, Reference Goldsby and Kohlstedt1997). In glaciers and ice sheets, Behn and others (Reference Behn, Goldsby and Hirth2021) have shown that while the traditional Glen flow relation using n = 3 may apply on a large scale, the appropriate parameters on a small scale (including through the thickness of the ice shelf) can vary substantially from n ≈ 2 to n ≈ 4 due to variations in grain size. The appropriate value for n can be related to the balance of deformation and recovery mechanisms, which is controlled by temperature and grain size; in warm ice with larger grain sizes, ice is likely to be less viscous due to the dominance of dislocation creep (Ranganathan and others, Reference Ranganathan, Minchew, Meyer and Peč2021).
Generally in natural ice which has been flowing for some time, the geometry of the CPO directly reflects the stress orientation and temperature. However, in samples from the Amery IS, significantly different CPO geometries have been observed in samples of meteoric and marine ice taken from the same depth of core, which cannot be explained by temperature differences alone (Treverrow and others, Reference Treverrow, Budd, Jacka and Warner2012). Frazil crystals are commonly preferentially oriented during accretion resulting in a CPO (Langhorne and Robinson, Reference Langhorne and Robinson1986; Wongpan and others, Reference Wongpan, Langhorne, Dempsey, Hahn-Woernle and Sun2015), however the vertical large-circle girdle fabrics observed in parts of the AM01 and G1 cores measured by Treverrow and others (Reference Treverrow, Budd, Jacka and Warner2012) are not consistent with this mechanism, suggesting there are additional mechanisms governing the microstructural evolution of marine ice. The presence of a strong CPO can increase strain rates in deforming ice by an order of magnitude in some orientations (Budd and Jacka, Reference Budd and Jacka1989; Budd and others, Reference Budd, Warner, Jacka, Li and Treverrow2013). In fact, Dierckx and others (Reference Dierckx, Peternell, Schroeder and Tison2014) examined the rheological behaviour of marine ice deformed in different orientations to the existing fabric and found that appropriate values of the stress exponent n varied from 2.1 and 4.1 with orientation. Differences in CPO between ice layers in an ice shelf are likely to have a significant effect on rheology and strain localisation in an ice shelf.
There is no clear consensus in the literature on the scale by which these differences in microstructure might affect the rheology of marine ice. In this study, we assume that the presence of a strong CPO and finer grain sizes in marine ice have a small weakening effect when compared with meteoric ice.
Marine ice in ice shelf models
Modelling investigations into the sensitivity of ice shelves to marine ice temperature and physical properties have been sparse. Lange and MacAyeal (Reference Lange and MacAyeal1986) conducted finite-element simulations of the Filchner-Ronne IS, in an attempt to determine if a scenario with an additional 10% weaker layer of marine ice at the base of the ice shelf (represented as an increase in thickness and depth-averaged flow enhancement) would agree with observations of the flow regime. They found that a scenario with a thicker shelf but no weakening had the best agreement, implying that the marine ice rheology itself had little effect, although the sparsity of observational data made it difficult to be conclusive.
Some models attempt to incorporate heterogeneous ice properties through varying viscosity, by spatially varying the rate factor B(T′). B(T′) is dependent on temperature, and can be related to A(T′) through:
meaning that a higher value for B(T′) represents a higher viscosity. Hulbe and others (Reference Hulbe, Johnston, Joughin and Scambos2005) varied the rate factor to incorporate varying marine and meteoric ice rheological properties into a 2-D numerical model of the Brunt Ice Shelf (hereafter Brunt IS). In the Brunt IS, marine ice exists not as a basal layer, but as a matrix connecting discrete blocks of meteoric ice. After inverting for B(T′) across the ice shelf based on expected temperature values, they also manually scaled it in some areas by a factor of 0.5–0.9 to represent the contrast in properties of marine and meteoric ice. However, they found that in the right-lateral shear margin of the Stancomb-Wills Ice Stream a factor of 10 reduction in B(T′) was required to match observations. This suggests that either the marine and meteoric ice are poorly connected in that area (in stark contrast to the rest of the ice shelf), or that the marine ice in that region is rheologically weak. However, the authors concluded that a range of combinations of marine ice thickness and enhancement factor would produce a result which agrees with observations, highlighting the need for better measurements of ice shelf geometry and rheological behaviour.
Khazendar and others (Reference Khazendar, Rignot and Larour2009) confirmed the rheological differences between marine and meteoric ice on the Brunt IS by using an inverse modelling method to estimate the value of B(T′) spatially from surface velocity data. They found that the predicted values of B(T′) were consistently $\sim \!30\percnt$
lower for areas filled with marine ice, compared to areas filled with meteoric ice. This can largely be explained by differences in temperature between the two ice types, assuming little thermal exchange between them. However, they acknowledge that crystallographic structure, salinity and impurity content of marine ice may also play a part which could not be separated from thermal effects in their study. Although they are a powerful tool, inverse modelling methods are unable to differentiate between the multiple factors which contribute to variations in the rheological properties of ice and do not capture how rheological and thermal properties can change over time.
All of the models mentioned here are depth-averaged, meaning that the effects of vertical variations in rheology are only captured in a bulk sense. All models also use a version of Glen's flow relation throughout the entire shelf, although Hulbe and others (Reference Hulbe, Johnston, Joughin and Scambos2005) theorise that the use of different flow relations for marine and meteoric ice may be more accurate.
The focus of this study is to separate the effects of temperature from other factors that influence marine ice rheological behaviour, and examine their separate and combined effects on the dynamics of a 3-D idealised ice shelf. We use a Blatter–Pattyn approximation to capture horizontal and vertical variations in rheology and temperature and use the ESTAR flow relation (Graham and others, Reference Graham, Morlighem, Warner and Treverrow2018) to incorporate the effects of changing stress configurations throughout the ice shelf.
Methods
We developed an idealised ice shelf domain to test the sensitivity of ice shelf dynamics to changes in enhancement factor and temperature in a basal layer, implemented using the Ice-sheet and Sea-level System Model (ISSM) (Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012). ISSM is a finite-element model which can solve a range of simplifications to the full-Stokes equations. All model runs in this study use a higher-order approximation to the full-Stokes equations, specifically the Blatter–Pattyn approximation (Blatter, Reference Blatter1995; Pattyn, Reference Pattyn2003), without evolving the thermal model and assuming no surface accumulation or basal melting.
A schematic diagram of the model initial conditions is shown in Fig. 2. The domain was a rectangular ice shelf measuring 200 km wide and 500 km long, with an initial thickness of 500 m at the grounding line, and 180 m at the calving front. These dimensions were chosen to roughly represent the Amery IS, with thickness reduced to achieve a steady-state thickness throughout the experiments. Horizontal and vertical mesh resolutions were chosen based on grid sensitivity studies (described in the Appendix), with areas of higher resolution localised at the calving front corners to minimise mesh anomalies at those points.
Schematic diagram of the modelled ice shelf. Solid black lines connect mesh vertices, and dashed lines show possible locations of the meteoric–marine ice interface. The ice shelf is embayed, with no-slip boundary conditions imposed on all nodes in the lateral walls (grey dots). Not to scale.
The inflow velocity was fixed across the ice shelf inflow boundary using the below equation :
where x is the horizontal distance across the front of the ice shelf, x max is the maximum shelf width and x mid is the midpoint of the ice shelf (x max/2). This gives a maximum inflow velocity at the centre of 600 ma−1, decreasing to zero at the sides of the ice shelf. The model is an embayed shelf, with no slip on the lateral boundaries and free flux at the calving front, and is assumed to be in hydrostatic equilibrium.
As discussed earlier, the ESTAR flow relation assumes that the microstructure has evolved to be compatible with the stress configuration at all points (i.e. the ice is in tertiary creep, with a steady-state CPO and grain size which will not change as long as the stress and temperature conditions remain the same). When that assumption is true, the effects of microstructure are implicitly captured in E(λS). Crucially for this study, the microstructure of marine ice can be different to what is typical for meteoric ice in a given stress configuration, as discussed in the Introduction. As ESTAR was developed to represent meteoric ice, rather than marine ice, we took account of the rheological differences in meteoric and marine ice through varying the enhancement factor as outlined below.
The assumption of microstructural compatibility can be inappropriate in scenarios where the stress regime is changing rapidly (as discussed at length by Graham and others (Reference Graham, Morlighem, Warner and Treverrow2018) and McCormack and others (Reference McCormack2022)). However, in an idealised ice shelf scenario similar to ours, Graham and others (Reference Graham, Morlighem, Warner and Treverrow2018) argue that due to gradual changes in the stress regime and contours of λS that are generally well aligned with the ice flow, compatible fabrics would likely be maintained, even in the progression towards the ice front.
A control experiment was defined using a constant value for the simple shear enhancement factor of E S = 8, with E C = (3/8) E S (Treverrow and others, Reference Treverrow, Budd, Jacka and Warner2012). The temperature in this control experiment was set to − 20°C everywhere at the surface and − 2°C everywhere at the base, with a linear temperature gradient between. The control simulation was initially run for 1000 years with time steps of 6 months, and the ice shelf velocity and geometry at the final time step from this spin-up run was used as the initial condition for all perturbation/sensitivity experiments.
We investigated the impact of changing temperature and enhancement factor using six sets of experiments, which we have assigned abbreviations as listed in Table 1. In some experiments the temperature was set to be isothermal at − 2°C in a marine ice basal layer, with a linear temperature gradient above, extending to − 20°C at the surface, and in others the temperature gradient was entirely linear, from − 2°C at the base to − 20°C at the surface (see Fig. 1). The basal layer was either weakened or strengthened by changing E S (see Eqns (4) and (5)) in increments of 0.3 by up to ± 2.4 (± 30%). In each set, a simulation was performed for each of 24 basal layer thickness values ranging from 0 to 120 m in increments of 5 m. These values are based on the range of possible thicknesses observed in marine ice basal layers (e.g. Eicken and others, Reference Eicken, Oerter, Miller, Graf and Kipfstuhl1994; Fricker and others, Reference Fricker, Popov, Allison and Young2001; Craven and others, Reference Craven, Allison, Fricker and Warner2009), and a range of enhancement factor values that we consider to be plausible given the limited available information on marine ice rheology outlined in the Introduction. Thermal profiles are based on those observed in major ice shelves (see Fig. 1).
List of all experimental sets and their conditions
Each set consisted of 24 simulations, covering a range of basal layer thicknesses from 5 and 120 m.
Our experimental design allows us to separate the effects of increased temperature in a basal marine ice layer from the effects of other ice properties such as microstructure and chemistry, through their implicit/assumed effect on the enhancement factor. All experiments were run with a time step of 6 months continuously for 100 years, by which time the velocity and thickness values had converged, such that they varied by $\ll \!1\percnt$
between time steps. Ice mass flux through the front of the ice shelf was calculated at the final time step of each experiment.
Results
Results from all model runs (final ice velocity, thickness, front flux and viscosity) can be found in Supplementary Tables S1 and S2. In many of the results presented here, we focus on scenarios where the basal layer is weakened or strengthened by ±15$\percnt$
(E S = 6.8 and E S = 9.2). We have chosen to focus on these values as they lie in the middle of a range of what could reasonably be expected based on our limited knowledge of marine ice rheology. Specifically, this was informed by the possible effects of CPO orientations (Duval and others, Reference Duval, Ashby and Anderman1983; Budd and Jacka, Reference Budd and Jacka1989; Dierckx and others, Reference Dierckx, Peternell, Schroeder and Tison2014), grain size (Cole, Reference Cole1987; Behn and others, Reference Behn, Goldsby and Hirth2021) and impurity content (Hammonds and Baker, Reference Hammonds and Baker2016; Craw and others, Reference Craw, Qi, Prior, Goldsby and Kim2018; Hammonds and Baker, Reference Hammonds and Baker2018) on the viscosity of polycrystalline ice, as well as the experiments of Dierckx and Tison (Reference Dierckx and Tison2013), as outlined in the Introduction.
Figure 3 shows final values for ice mass flux through the front of the ice shelf (a) and thickness at the centre front (b) for the control experiment, and scenarios with an isothermal basal layer, and $15\percnt$
weakening or strengthening in the basal layer. Changes in flux and thinning in the isothermal basal layer scenarios ($T_{\rm i}^{E + }$
and $T_{\rm i}^{E-}$
) are much greater than in the linear temperature gradient scenarios ($T_{\rm l}^{E + }$
and $T_{\rm l}^{E-}$
). For example, a simulation with a basal layer thickness comprising 20$\percnt$
of the ice shelf volume with a $15\percnt$
softened basal layer predicts an ice mass flux of 34.3 Gt a −1 and front thickness after 100 years of 192 m when the temperature gradient is linear ($T_{\rm l}^{E + })$
, compared with 40.8 Gt a −1 flux and 182 m front thickness when the basal layer is isothermal ($T_{\rm i}^{E + }$
). Thinning in the $T_{\rm i}^{E + }$
and $T_{\rm i}^{E-}$
scenarios increases mildly non-linearly with increased basal layer volume.
Final values of flux through the front of the ice shelf (a) and thickness at the centre front (b) after 100 years for a homogeneous shelf, and for scenarios with $15\percnt$
strengthening or weakening in the basal layer. Flux and thickness are given as a percentage change from the control experiment. The dashed grey line is the control experiment, the solid grey line is a shelf with homogeneous E S and an isothermal basal layer, and green and orange lines signify combinations of strengthening, weakening and temperature gradient as defined in Table I.
Final centre front thickness and front flux values for all simulations are shown in Fig. 4. Again, the greatest differences in flux and thickness are seen in the isothermal scenarios (panels a and c), with flux increasing by up to 46% and thickness decreasing by up to 14% compared with the homogeneous ice shelf case. This is a result of a decrease in average ice viscosity in the ice shelf of up to 45%. In the linear temperature gradient scenarios (panels b and d), flux increased by up to 4.1% in an $T_{\rm l}^{E + }$
scenario, or decreased by up to 6.2% in an $T_{\rm l}^{E-}$
scenario. Thickness at the front decreased by up to 1.1% in the $T_{\rm l}^{E + }$
scenarios, and increased by up to 1.4% in the $T_{\rm l}^{E-}$
scenarios. This is a result of a change in average viscosity in the ice shelf by up to ± 8%.
Change in final values of ice mass flux through the front of the shelf (top) and final thickness at the centre front (bottom) after 100 years for all simulations, when compared with a shelf composed on entirely meteoric ice with a linear temperature profile and homogeneous enhancement factor: (a, c) where the temperature is set to − 2°C throughout the basal layer; (b, d) where the temperature gradient is linear for all experiments.
The stress distribution in the ice shelf varies between simulations, as shown in Fig. 5. In the linear temperature gradient scenarios, an increase in enhancement factor of $15\percnt$
in the basal layer causes a slight ($\sim \!3\percnt$
) decrease in effective deviatoric stress at the base of the ice shelf, and a corresponding smaller ($\sim \!1\percnt$
) increase in effective deviatoric stress at the surface of the ice shelf. In isothermal basal layer scenarios the effective deviatoric stress remains constant through much of the basal layer, and is markedly (20–30%) greater at the surface than in the homogeneous scenario. Deviatoric stresses are greater at the edge of the ice shelf, where it has thinned and so the proportion of marine ice (which has a constant thickness throughout the ice shelf) is greater.
Vertical profiles of effective deviatoric stress in the ice shelf after 100 years, from the same scenarios plotted in Fig. 3: (a) averaged over an area in the centre of the ice shelf (compression-dominated), and (b) averaged over an area near the edge of the ice shelf (shear-dominated). The approximate position of the basal layer is shaded in grey.
Discussion
The effect of marine ice on ice shelf dynamics
Our results using a Blatter–Pattyn ice-sheet model demonstrate that differences in temperature and enhancement factor on the scale of those which have been observed or inferred in marine ice can have a significant effect on the dynamics of the ice shelf. Because we have varied the temperature profile and enhancement factor separately, we are able to assess the relative impact of changes in these properties on shelf dynamics. As the thermal and physical properties of an ice shelf with a marine ice layer can vary substantially, this is important for understanding how the specific physical and thermal characteristics of an ice shelf may control its future behaviour.
Altering the temperature profile of the modelled shelf to be isothermal in the basal layer, when compared to a shelf with a simple linear temperature gradient, has an order of magnitude greater effect on viscosity and consequently thickness and mass flux than changing the enhancement factor in the basal layer alone. This implies that the effect of marine ice on the thermal profile is a much greater control on ice shelf dynamics than the combined effects of marine ice microstructure and chemistry, which are represented here by changes in the flow enhancement factor.
Most ice shelf models treat all ice in the ice shelf as meteoric. In ice shelves with a marine ice layer, our results suggest that differences in viscosity in that layer may cause the ice velocity, and therefore flux and thickness, to diverge from those predicted by models which do not explicitly include marine ice. In ice shelves where the marine ice layer is not isothermal, as is the case in parts of the Filchner-Ronne IS (see Fig. 1), our results suggest that the physical properties of the marine ice can potentially contribute to a change in ice mass flux of ± 2–6% (the direction of this change depends on whether the marine ice is weakened or strengthened relative to the meteoric ice). To estimate the possible scale of this effect in reality we consider the Filchner-Ronne IS, which has an estimated ice front flux of − 237 Gt a−1 (Moholdt and others, Reference Moholdt, Padman and Fricker2015), and marine ice present beneath approximately one-quarter of the ice shelf (Lambrecht and others, Reference Lambrecht, Sandhäger, Vaughan and Mayer2007). If marine ice impacts the Filchner-Ronne IS in a similar way as modelled in this study, and its thickness does not change significantly over time, this could represent an uncertainty in model-derived estimates for future ice front flux on the order of ± 1–3 Gt a−1 (0.5–1%), and future thickness on the order of ± 0.5–1 m (0.1–0.4%) after 100 years.
In ice shelves where the marine ice layer is isothermal, such as the Amery IS, the temperature in ice above the basal layer is also increased relative to a linear gradient through the entire thickness. This change in thermal profile would dramatically reduce the stiffness of the ice, leading to much greater increases in flux and shelf thinning. The Amery IS has an estimated front flux of − 44 Gt a−1 (Fricker and others, Reference Fricker, Warner and Allison2000), and marine ice present beneath approximately one-third of the ice shelf (Fricker and others, Reference Fricker, Popov, Allison and Young2001). Our results suggest that future front flux in the Amery IS could be underestimated by up to 10 Gt a−1 (22%), and thickness could be overestimated by up to 35 m ($7\percnt$
) in areas underlain by marine ice after 100 years by models in which the marine ice layer is not accounted for. These differences could have far-reaching consequences for our predictions of ice dynamics in the polar regions, as the buttressing effect of the ice shelf on flowing ice upstream could in turn be overestimated.
In addition to its effects on ice mass flux and shelf thickness, our results predict that the presence of marine ice will alter the distribution of stress in the ice shelf, with potential implications for rift formation and calving. Marine ice accretion has been shown to have a stabilising effect on an ice shelf by filling basal rifts from below (Khazendar and others, Reference Khazendar, Tison, Stenni, Dini and Bondesan2001; Khazendar and Jenkins, Reference Khazendar and Jenkins2003), and arresting fracture propagation where it exists in suture zones (Rignot, Reference Rignot1998; Jansen and others, Reference Jansen, Luckman, Kulessa, Holland and King2013; Kulessa and others, Reference Kulessa, Jansen, Luckman, King and Sammonds2014). However, our results indicate that when vertical changes in rheology are taken into account, in a shelf with a basal layer which is less viscous than the ice above, stress can become more concentrated at the surface. This could potentially increase the probability of surface fracture.
Limitations and further work
Here we will consider the main limitations of this study, discuss how these simulations could be improved, and detail how the results of this and future studies can be used to improve the accuracy of ice shelf and ice-sheet models.
There are many assumptions underlying our model design, both in the way we incorporated marine ice physical properties, and in the set-up of the model itself. As outlined in the Introduction, our knowledge of the in situ temperature, chemical composition and microstructural characteristics of marine ice are limited. We have chosen two representative thermal profiles for a shelf with a marine ice layer, but as Grosfeld and Thyssen (Reference Grosfeld and Thyssen1994) suggest in their discussion of the Filchner-Ronne IS, it is possible for the thermal profile to vary within a shelf as the marine ice accretion rate changes. Our ‘isothermal’ profile is also a simplification of the more S-shaped profile seen in the Amery IS. In fact, a shelf with a high surface accumulation rate might have a near-isothermal profile close to the surface, lowering the bulk temperature of the ice shelf and negating some of the effects of the marine ice layer. As we have identified the thermal profile of the ice shelf as being the most important factor for predicting the influence of marine ice on its dynamics, further investigation of ice shelf temperature profiles and how they vary spatially should be a priority.
Measurements of marine ice microstructure and chemical content are also extremely sparse. We have attempted to investigate the effects of strengthening and weakening in the marine ice that we consider to be possible as a result of these properties based on existing observations and modelling studies (e.g. Budd and Jacka, Reference Budd and Jacka1989; Rommelaere and MacAyeal, Reference Rommelaere and MacAyeal1997; Craven and others, Reference Craven, Allison, Fricker and Warner2009; Dierckx and Tison, Reference Dierckx and Tison2013; Dierckx and others, Reference Dierckx, Peternell, Schroeder and Tison2014; Hammonds and Baker, Reference Hammonds and Baker2016, Reference Hammonds and Baker2018; Craw and others, Reference Craw, Qi, Prior, Goldsby and Kim2018; Behn and others, Reference Behn, Goldsby and Hirth2021) but these may vary widely across an actual ice shelf. Further laboratory studies to measure the microstructural characteristics, chemical content and mechanical response to stress of marine ice, particularly in tertiary creep, would be extremely valuable for improving the way ice flow relations capture those effects.
Our method for incorporating rheological differences between marine and meteoric ice as a result of microstructure and chemistry, by scaling the enhancement factor to incorporate both, cannot take into account any separate effects they might have. The ESTAR flow relation is based on the assumption that the microstructure of ice is compatible with the underlying stress conditions (meaning that the CPO and grain size are consistent with the expected steady-state microstructure for meteoric ice for that stress magnitude and configuration). We assume that by scaling the enhancement factor we are, in addition, incorporating the effects of any microstructural characteristics which are not compatible (i.e. the grain size and CPO of marine ice, which can differ from those in meteoric ice as outlined in the Introduction). This overlooks the fact that the new microstructure may have an anisotropic effect on strain rates. The ESTAR relation also does not incorporate any stress dependence of E S and E C, suggested by the results of experiments performed by Treverrow and others (Reference Treverrow, Budd, Jacka and Warner2012). Once the relationships between E S and E C and a range of combined stress configurations has been more precisely defined through experiments, this should be incorporated into the formulation of the flow relation.
Changing the enhancement factor of a marine ice layer can scale the strain rates predicted by the flow relation linearly, but cannot capture any non-linear effects. In the future when more is known about the physical properties of marine ice, it may be more appropriate to vary other components of the flow relation spatially, such as the stress exponent n. Appropriate values for n have been shown to diverge significantly from its commonly accepted value of 3 in different settings (Treverrow and others, Reference Treverrow, Budd, Jacka and Warner2012; Qi and others, Reference Qi, Goldsby and Prior2017; Bons and others, Reference Bons2018; Zeitz and others, Reference Zeitz, Levermann and Winkelmann2020). If the effects of different physical properties of ice such as grain size can be separated from other microstructural and chemical characteristics, appropriate values for n or different formulations of the flow relation can be inferred based on specific studies of that property (e.g Behn and others, Reference Behn, Goldsby and Hirth2021; Qi and Goldsby, Reference Qi and Goldsby2021).
Despite these limitations, the main findings from this study can still be incorporated into ice shelf models. We show that in order to understand the impacts of marine ice on ice shelf dynamics, a 3-D model is ideal to capture the effects of the spatial distribution of ice properties, particularly temperature. This is especially true in cases where stress distribution and fracture are of particular interest. It is important to consider the thickness of the marine ice layer, the accumulation rate and the thermal profile of the ice shelf and adjust these to reflect the most likely scenario present in the ice shelf based on the observational data available.
However, running a 3-D model can be computationally expensive. Where it is necessary to use a depth-averaged model, it is still important to consider the 3-D distribution of ice physical properties and temperature, and how those might affect bulk viscosity. For an ice shelf like the Amery, with areas where the marine ice layer comprises 30% of its thickness and is near isothermal, the depth-averaged temperature should be increased by ~ 30% and the surface deviatoric stress by ~20–30% in those areas. In an ice shelf such as the Filchner-Ronne, with areas where the marine ice layer is equally thick but the thermal gradient is close to linear, it may be appropriate to increase the enhancement factor or stress exponent, or decrease the rate factor slightly to provide a weakening effect. The magnitude of this weakening effect will depend on the known or estimated physical properties of marine ice at that location.
Although inversion may implicitly capture some of the effects of temperature and physical properties discussed here, it does not capture how these properties may vary in their spatial distribution over time. As the polar regions undergo increasingly rapid changes, it is crucial to understand the constantly evolving physical and thermal structures of ice shelves, in order to predict their effect on ice-sheet mass balance in the future.
Conclusions
We found that of the two major controls on marine ice rheology investigated here – (1) thermal profile and (2) mechanical properties controlled by microstructure and chemistry – the thermal profile has a much greater control on ice viscosity, and therefore ice velocity and mass flux. When parameterising an ice shelf model, it is therefore more crucial to properly constrain the thermal profile of the ice shelf than to focus on the spatial distribution of ice mechanical properties. However, in cases where the thermal gradient is not significantly affected by the marine ice layer, its mechanical properties still influence dynamics, although to a lesser extent.
In an idealised ice shelf model with a basal marine ice layer of constant thickness, the presence of the marine ice can cause an increase in flux of up to 46% and an increase in effective deviatoric stress of up to 30% at the surface when compared with a homogeneous shelf. This is an upper estimate; in a real scenario, where marine ice varies in thickness and spatial distribution, the effect is likely to be smaller. Although some of the effects of marine ice can be implicitly captured in models through inversion, changes in deviatoric stress with depth and changes in marine ice thickness over time, which are not well represented in models, will have an effect on ice mass flux and the probability of surface fracture.
Future work should prioritise observations of marine ice thickness across different ice shelves, and how its thermal, chemical and microstructural properties vary in time and space. These data will be essential in informing the development of ice flow relations, including appropriate parameter values, and reducing uncertainty about the evolution of ice shelf dynamics into the future.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2022.66
Acknowledgements
We would like to thank Matthew Siegfried, and two anonymous reviewers for their helpful comments which greatly improved the manuscript. Lisa Craw is supported by an Australian Government Research Training Program Scholarship at the University of Tasmania. Felicity McCormack is supported by an Australian Research Council (ARC) Discovery Early Career Research Award DE210101433. Sue Cook received grant funding from the Australian Government as part of the Antarctic Science Collaboration Initiative programme. Part of the research work undertaken was supported by the Australian Research Council's Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001).
Appendix
We performed grid sensitivity tests to assess the effects of the horizontal and vertical mesh resolutions on model output. Using the same initial conditions as the homogeneous spin-up scenario described in the Methods, 14 mesh sensitivity simulations were run across a range of resolutions, and the results of these are shown in Supplementary Table S3. We varied the maximum horizontal vertex edge length (h max) from 3000 m to 10 000 m in increments of 1000 m (excluding the front corners where it was reduced to 500 m in all experiments to minimise an anomaly at those points), in a mesh with 20 vertical layers. Vertically, the number of mesh layers was varied from 5 to 30 in increments of 5, with a horizontal resolution of h max = 5000 m. Each of these simulations was run for 100 years from the same initial conditions, by which time the values of ice mass flux across the front of the shelf had converged to within 2$\percnt$
per time step. The final front flux values for all of these simulations are shown in Fig. 6. A horizontal resolution of h max = 5000 m and a vertical resolution of 20 layers were chosen for all further experiments in this study. At the chosen resolution, the results of the model should be independent from resolution to within $\ll \!0.5\percnt$
.
Final front flux values from a series of mesh convergence tests used to select an appropriate horizontal and vertical grid resolution for the simulations. All horizontal resolution tests were run with 20 vertical layers, and all vertical resolution tests were run with h max = 5000 m. The values chosen for use in this study are outlined in orange.
Open access
