Modelling the influence of marine ice on the dynamics of an idealised ice shelf

Abstract

Understanding the dynamic behaviour of ice shelves, specifically the controls on their ability to buttress the flow of ice into the ocean, is critical for predicting future ice-sheet contributions to sea level rise. Many large ice shelves, which are predominantly composed of meteoric ice, have a basal layer of marine ice (formed from accumulated platelets at the ice–ocean interface), comprising up to 40% of their thickness locally. Differences in temperature, chemistry and microstructure between marine and meteoric ice mean the rheological properties of the ice vary throughout the ice shelf. These differences are not explicitly accounted for in ice-sheet modelling applications, and may have an important influence on ice shelf dynamics. We tested the sensitivity of a model of an idealised ice shelf to variations in temperature distribution and flow enhancement, and found that incorporating a realistic thermal profile (where the marine ice layer is isothermal) had an order of magnitude greater effect on ice mass flux and thinning than incorporating the mechanical properties of the marine ice. The presence of marine ice at the ice shelf base has the potential to significantly increase deviatoric stresses at the surface and ice mass flux across the front of an ice shelf.

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, 1986; Jenkins and Doake, 1991; Jenkins and Bombosch, 1995; Bombosch and Jenkins, 1995). This marine ice layer commonly comprises between a sixth and a third of the total ice shelf thickness (Oerter and others, 1992; Eicken and others, 1994; Treverrow and others, 2010), with its thickness varying spatially across the base of the ice shelf (Fricker and others, 2001). 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, 1994; Treverrow and others, 2010; Dierckx and others, 2014).

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, 2005), cause stresses to localise in the upper part of the shelf (e.g. Lange and MacAyeal, 1986) and affect rifting behaviour (e.g. Kulessa and others, 2014). In turn, changes in the ice shelf properties can have far-reaching consequences for buttressed ice upstream (Fürst and others, 2015; Pegler, 2018; Reese and others, 2018). 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, 1952, 1955, 1958):

(1)$$\dot{\varepsilon} = A( T') \tau_{\rm e}^{n-1}\sigma' ,$$

where $\dot {\epsilon }$ is the strain rate tensor (s⁻¹), 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, 2012; Qi and others, 2017; Bons and others, 2018; Millstein and others, 2022). In this form, the law assumes that ice is mechanically isotropic. More recent studies (e.g. Treverrow and others, 2012; Budd and others, 2013) have considered the addition of an enhancement factor, E, to Eqn (2), as follows:

(2)$$\dot{\varepsilon} = EA( T') \tau_{\rm e}^{n-1}\sigma',\; $$

where E is defined as the ratio of tertiary to secondary creep rates

(3)$$E = {\sigma_{\rm ter}\over \sigma_{\rm sec}}.$$

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 (2013) 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 (2018):

(4)$$\dot{\varepsilon} = E( \lambda_{\rm S}) A( T') \tau_{\rm e}^{2}\sigma' ,\; $$

where E(λ_S) is defined as:

(5)$$E( \lambda_{\rm S}) = E_{\rm C} + ( E_{\rm S} - E_{\rm C}) \lambda_{\rm S}^{2} .$$

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 (2012) 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:

(6)$$\lambda_{\rm S} = {\Vert \tau'\Vert \over \tau_{\rm e}}.$$

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 (2022).

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 (2018) and McCormack and others (2022).

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, 1972; Durham and others, 2001; Goldsby and Kohlstedt, 2001; Qi and others, 2017), or on meteoric ice (e.g. Gao and Jacka, 1987; Dahl-Jensen and others, 1997; Treverrow and others, 2012; Craw and others, 2018). 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, 2009; Wang and others, 2022). 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 (2009) estimate (using the ice deformation temperature-dependence relationship outlined by Budd and Jacka, 1989) 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 (1994) 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 (1994) 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⁻¹ 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.

Fig. 1. Ice shelf temperature profiles: (a) from borehole measurements on the Amery IS, with an isothermal basal marine ice layer (Craven and others, 2009; Wang and others, 2022); (b) from borehole measurements in the Filchner-Ronne IS (Eicken and others, 1994), 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 (1994) 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, 2009).

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, 1997), 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, 1994; Khazendar and others, 2001)), 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, 2021), however the precise magnitude and nature of this effect is not well understood.

Trace amounts of H₂SO₄ and NH₃, when added at concentrations commonly seen in meteoric ice, have been observed to weaken synthetic ice deformed in the laboratory (Li and others, 2009; Hammonds and Baker, 2018), whereas a similar concentration of Ca²⁺ has been observed to strengthen it (Hammonds and Baker, 2016). 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, 1987), 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 (1997), 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, 2009; Kulessa and others, 2014; King and others, 2018). 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 (2013) 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, 1989). 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 (2018) 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, 1986; Obbard and Baker, 2007; Eichler and others, 2017).

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 (2018).

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, 1989).

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, 1994; Treverrow and others, 2010). 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, 1997; Obbard and Baker, 2007). Bubble content can also limit the rate of grain growth in ice by inhibiting grain boundary migration (Azuma and others, 2012). Finer grain sizes in ice generally lead to rheological weakness, due to an increase in grain-size sensitive creep mechanisms (Paterson, 1991; Goldsby and Kohlstedt, 1997). In glaciers and ice sheets, Behn and others (2021) 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, 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, 2012). Frazil crystals are commonly preferentially oriented during accretion resulting in a CPO (Langhorne and Robinson, 1986; Wongpan and others, 2015), however the vertical large-circle girdle fabrics observed in parts of the AM01 and G1 cores measured by Treverrow and others (2012) 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, 1989; Budd and others, 2013). In fact, Dierckx and others (2014) 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 (1986) 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:

(7)$$B( T') = A( T') ^{{-1}/{n}} ,\; $$

meaning that a higher value for B(T′) represents a higher viscosity. Hulbe and others (2005) 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 (2009) 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 (2005) 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, 2018) 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, 2012). 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, 1995; Pattyn, 2003), 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.

Fig. 2. 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 :

(8)$$v = 600\times\exp{\left(-{x-x_{\rm mid}\over {x_{\max}}/{2.5}}\right)^4} ,\; $$

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⁻¹, 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 (2018) and McCormack and others (2022)). However, in an idealised ice shelf scenario similar to ours, Graham and others (2018) 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, 2012). 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, 1994; Fricker and others, 2001; Craven and others, 2009), 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).

Table 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, 1983; Budd and Jacka, 1989; Dierckx and others, 2014), grain size (Cole, 1987; Behn and others, 2021) and impurity content (Hammonds and Baker, 2016; Craw and others, 2018; Hammonds and Baker, 2018) on the viscosity of polycrystalline ice, as well as the experiments of Dierckx and Tison (2013), 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 ⁻¹ 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 ⁻¹ 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.

Fig. 3. 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%.

Fig. 4. 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.

Fig. 5. 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⁻¹ (Moholdt and others, 2015), and marine ice present beneath approximately one-quarter of the ice shelf (Lambrecht and others, 2007). 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⁻¹ (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⁻¹ (Fricker and others, 2000), and marine ice present beneath approximately one-third of the ice shelf (Fricker and others, 2001). Our results suggest that future front flux in the Amery IS could be underestimated by up to 10 Gt a⁻¹ (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, 2001; Khazendar and Jenkins, 2003), and arresting fracture propagation where it exists in suture zones (Rignot, 1998; Jansen and others, 2013; Kulessa and others, 2014). 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 (1994) 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, 1989; Rommelaere and MacAyeal, 1997; Craven and others, 2009; Dierckx and Tison, 2013; Dierckx and others, 2014; Hammonds and Baker, 2016, 2018; Craw and others, 2018; Behn and others, 2021) 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 (2012). 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, 2012; Qi and others, 2017; Bons and others, 2018; Zeitz and others, 2020). 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, 2021; Qi and Goldsby, 2021).

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.