Ice dynamics

We use the BISICLES ice-sheet model in this study (Cornford and others, Reference Cornford2013). The momentum equation is based on the vertically integrated ‘L1L2’ approximation (Schoof and Hindmarsh, Reference Schoof and Hindmarsh2010) and is suited to fast flowing ice streams and ice shelves. Numerically, the model employs a finite volume method with AMR (adaptive mesh refinement), which allows the use of non-uniform, time-dependent meshes. Constructing meshes with coarse resolution in the slowly flowing inland areas and fine meshes around the grounding line and fast flowing areas allows us to capture the behaviour of outlet glaciers without wasting computational resources on the vast inland ice sheet. In this study we implement three refinement levels above the base resolution of 4 km with the mesh gradually refined toward the grounding line, where resolution is 0.5 km (Appendix). Time step varies to satisfy the Courant-Friedrichs-Lewy condition everywhere, meaning ~32 time steps a⁻¹ in Arura Basin. The mesh is updated at each time step.

The 2-D stress-balance equation is

(1) $$\nabla \cdot \left[{\rm \phi} h{\bar {\rm \mu}} \left(2 {\dot {\epsilon}} + 2tr({\dot {\epsilon}}){\bi I} \right) \right] + {\rm \tau} ^b = {\rm \rho}_igh\nabla s$$

where h is the thickness, $\dot{\epsilon} $ is the strain rate tensor, ρ _i is the ice density, ϕh is vertically integrated effective viscosity, s is surface elevation and g is gravitational acceleration. We use two formulations for τ ^b: linear and nonlinear forms of Weertman friction law in the experiments

(2) $${\rm \tau} ^{\rm b} = \left\{ {\matrix{ { - C{\left \vert u \right \vert} ^{\left( {\displaystyle{1 \over m} - 1} \right)}u\quad \;if\; \displaystyle{{{\rm \rho}_i} \over {{\rm \rho}_{\rm w}}}h \gt - r} \cr {\quad\quad\enspace 0\quad \quad \ \qquad {\rm otherwise}} \cr}} \right.$$

in which ρ _i and ρ _w are ice and ocean water density, respectively, r is bed elevation, h is ice thickness, u is horizontal ice velocity and the exponent, m = 1 for the linear friction law and m = 3 for the nonlinear friction law. The friction coefficient C that contributes to basal traction τ ^b and a stiffening factor ϕ that contributes to the viscosity are calculated by solving an inverse problem, minimizing the difference between modelled velocity and the observed surface velocity (e.g. MacAyeal, Reference MacAyeal1993; Joughin and others, Reference Joughin2009; Morlighem and others, Reference Morlighem2010). The inverse method is described in detail in the study of Cornford and others (Reference Cornford2015).

We carry out each simulation on a 1792 × 1792 km² domain (Fig. 1) taking bedrock topography and ice thickness at 1 km resolution from Bedmap2 (Fretwell and others, Reference Fretwell2013). The sea bed topography data beneath the Totten and Dalton ice shelves is poorly constrained in Bedmap2. We therefore modify the bathymetry following Gwyther and others (Reference Gwyther, Galton-Fenzi, Hunter and Roberts2014). This means we deepen the cavity by a depth R _C. According to Greenbaum and others (Reference Greenbaum2015), it is an average of 500 m deeper along the long axis of Totten cavity. So instead of setting R _C = 300 m as Gwyther and others (Reference Gwyther, Galton-Fenzi, Hunter and Roberts2014), we set R _C = 500 m along the cavity centreline, while along the grounding line, R _C = 0. Within the cavity, R _C is linearly interpolated between the cavity centreline and grounding line. In this study, parameterized melt rate depends only upon the depth of ice shelf bottom, so R _C will not affect the sub-ice shelf melt, but deepening the sea bed removes some very thin cavities beneath the ice shelf that would otherwise lead to immediate grounding. The 3-D temperature field used in the simulations is taken from Pattyn (Reference Pattyn2010) and is fixed over time. The surface mass-balance field is taken from Arthern and others (Reference Arthern, Winebrenner and Vaughan2006). The calving front is fixed in this study. Finally, the basal traction and stiffening factor are determined from surface velocity data (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011) with one notable exception: there is no data in the eastern part of Law Dome, so we set the basal traction coefficient there to the similar value as in the western half so that the bedrock elevation and surface slopes are similar (Fig. 2b).

Fig. 1.

The simulation region with bed elevation (colour bar) showing Totten Glacier (TG), Dalton Glacier (DG), Vanderford Glacier (VG), Law Dome (LD) and Reynolds Trough (RT). The insert shows a map of Antarctica with the study region outlined in black.

Fig. 2.

The initial conditions for experiments after 50 a of surface relaxation around Totten, Dalton and Vanderford glaciers, (a) μ_e is the effective viscosity coefficient, (b) β_e is the logarithm of basal friction coefficient log10 C, (c) M ₀ is the initial ice shelf basal melt rate, (d) $\left \vert\vec u\right \vert$ is the magnitude of surface velocity.

We find high-amplitude short-wavelength thickening or thinning at the initial stage of simulation. It presents a strong signal of change early in a simulation that results from non-physical causes and can mask the response to future climate anomalies. The magnitude of fluctuations is as large as 300 m a⁻¹ in places. We assume this unrealistic behaviour is due to discrepancy between the initial thickness and velocity, caused by the different epochs for measurements of geometry and ice velocity, interpolation method or undersampling. We attempt to address this issue by relaxing the geometry (Gong and others, Reference Gong, Cornford and Payne2014; Wright and others, Reference Wright2014). After 15 a almost all the grounded ice has thickness changes <10 m a⁻¹, and the grounding line has stabilized, and the ice thickness differs from bedmap2 by <100 m. A major exception is the southern flank of Totten glacier, which thins throughout the relaxation and indeed through every subsequent experiment. The bedrock in this region is only sparsely covered by airborne radar (Fretwell and others, Reference Fretwell2013), and so we assume the initial ice thickness is incorrect. The knock-on effect on other parts of the basin appears to be small – for example, as we will show, retreat along the trough between the Totten and Vanderford glaciers is conditional on the future ocean forcing, and does not take place without elevated melt rates. In common with Cornford and others (Reference Cornford2015) and Gong and others (Reference Gong, Cornford and Payne2014), we begin simulations with variable forcing after 50 a of relaxation with present-day climate, by which time the mean thickness change in the domain is 2.4 m.

Melt-rate parameterizations

Climate forcing from atmosphere and ocean impact ice-sheet dynamics by surface snow accumulation and sub-ice shelf melt, and the ice-sheet model responds to the two perturbations independently (Cornford and others, Reference Cornford2015). Dynamical thinning is the major mechanism by which mass is lost, and that is projected to be primarily forced by changes in the ocean for the next 2 centuries. So, in this study we keep the surface mass balance constant and consider only variation in the ice shelf basal melt rate. The ‘melt rate’ refers to ‘ice shelf basal melt rate’, and upper surface melt rates are not discussed further in the current study.

Sub-ice shelf melt rate is influenced by ocean conditions and the ice draft. Ideally, the simulation of melt rate would be carried out within a coupled ocean/ice shelf/ice-sheet model, but lacking such a model we examine the response of the ice-sheet model to simple melt-rate formulae constructed in light of the results of standalone ocean models. Each melt-rate formula is composed from a depth-dependent part based on a model of the present day, denoted by M, and a time-dependent part derived from emission driven climate projections, denoted by M _a.

Present-day melt rates are parameterized based on an ocean model of the Totten ice shelf region (Gwyther and others, Reference Gwyther, Galton-Fenzi, Hunter and Roberts2014), which was re-run with forcing conditions updated to cover the 1992–2012 period. The time-averaged melt rates are consistent with a parameterization of the melt rate as a function of ice shelf thickness. Figure 3a shows the fit between the parameterized melt rate and the modelled melt rate. Figure 2 shows the reconstructed melt rate based on the average modelled melt rates between 1992 and 2012 in the Aurora coastal area. We parameterize this data with a simple, depth dependent formulation (Asay-Davis and others, Reference Asay-Davis2015):

(3) $$M = Ad(T - T_{\rm f})\tanh \left( {\displaystyle{w \over {H_{\rm c}}}} \right),$$

where d is the depth of ice shelf bottom, T _f is the freezing point. T is the far field ocean temperature. w is the thickness of water column beneath ice shelf. H _c = 75.0 m is the reference ocean cavity thickness. A satisfactory fit for (3) is found setting T = −1.5°C, and A = −7.4 × 10⁻⁷ a⁻¹°C⁻¹, of which RMSE between the parameterized and mean melt rate is 3.8 m a⁻¹ (Fig. 3). In a previous study the parameterisation method is calibrated using the POP ocean model with idealized geometry (Asay-Davis and others, Reference Asay-Davis2015). However, the melt-rate parameterisation does not include ice/ocean interface slope, or any changes in hydrography resulting from external forcing, such as climate change. Therefore parameters used in this study are not likely to be appropriate for other ice shelves, or even this ice shelf for long duration simulations. Figure 3 shows the spatial variability of melt rates (blue dots). The parameterized melt rate (black line) fits the mean melt rates (red dots) of ice shelf thickness, this implies the parameterization will not well capture 2-D variability in melt rates. Considering that the highest melt rate is around the grounding line, which is important to ice dynamics and might be underestimated here, we also consider an experiment with higher melt-rate anomaly.

Fig. 3.

Parameterization of melt rate. (a) Simulated average melt 1980–2012 from ROMS and its parameterization (Eqn. (3)) as a function of ice shelf thickness. Mean melt rate at a given thickness is the average over all grid cells where the ice thickness is between thickness −Δd and thickness +Δd. 210 values of thickness between 86 m and 2644 m are taken. (b) Time evolution of melt rate anomalies produced by temperature anomalies from FESOM $(M_{\rm a} = 16{\rm \; m\;} {\rm a}^{ - 1}{\rm ^\circ} {\rm C}^{ - 1} \times \Delta T)$ driven by different climate models under different scenarios and the parameterization of Eqns (4) and (5).

Future anomalies are derived from the Ice2Sea FESOM simulations (Timmermann and others, Reference Timmermann2009; Timmermann and Hellmer, Reference Timmermann and Hellmer2013), which were driven by the two global climate models (with 1.25°–1.5° resolution for the ocean): HadCM3 and ECHAM5 in turn driven by two emissions scenarios: A1B and E1. Projections span 1980–2100 for ECHAM5 under both A1B and E1, 1980–2150 for HadCM3 under E1, and 1980–2200 for HadCM3 under A1B. FESOM is not able to resolve the small ice shelves along the Sabrina Coast, much as it could not resolve the Amundsen Sea ice shelves, as pointed out in Cornford and others (Reference Cornford2015). So rather than employ FESOM melt-rate anomalies directly we follow Cornford and others (Reference Cornford2015), who assumed a simple relationship between temperature anomalies ΔT in the nearby ocean and melt-rate anomalies M _a = 16 m a⁻¹°C⁻¹ × ΔT at the upper end of the range of observations and model studies (Rignot and others, Reference Rignot2002; Holland and others, Reference Holland, Jenkins and Holland2008). This linear relationship between melt-rate anomaly and ΔT might cause some errors especially when the warming trend is significant. Holland and others (Reference Holland, Jenkins and Holland2008) pointed out that, for ice shelf with high melt rate, the relationship between melt rate and thermal driving is nonlinear. Ocean temperature is characterized by averaging FESOM temperatures over a volume bounded laterally between the present day ice front, the 1000 m ocean depth contour and the extreme longitudes of the Aurora Basin and vertically between 200 and 800 m below sea level. This volume samples water masses that can both cross the continental shelf and access most of the sub-ice shelf cavity. FESOM temperatures differ considerably between the two climate models, but for each model, independent of climate scenarios. When driven by HadCM3, FESOM temperature rise 1°C by 2150 then remain almost constant, while under ECHAM5 they steadily decrease. Therefore we parameterize melt-rate anomalies over time as:

(4) $$M_{\rm a} = \left\{ {\matrix{ {7.2 \times {10}^{ - 4}t^2 - 0.027t + 0.7\; \; {\rm m}\; \,{\rm a}^{ - 1},} & {\; t \le 170\; \,{\rm a}} \cr {16\; \,{\rm m}\,{\rm a}^{ - 1},} & {{\rm otherwise}} \cr}} \right.,$$

for HadCM3, and

(5) $$M_{\rm a} = - 0.01t - 0.2\; {\rm m\;} {\rm a}^{ - 1},$$

for ECHAM5, where t is the time in years since 1980 in Eqns (4) and (5) (Fig. 3b). Figure 3b shows the fit of parameterized melt rate (lines) to melt rates calculated based on temperature anomalies produces by FESOM ocean model (scatters). The RMSE for ECHAM5 is 0.4 m a⁻¹ under the E1 scenario and 0.6 m a⁻¹ under the A1B scenario. The RMSE for HadCM3 is 1.3 m a⁻¹ for E1 scenario and 1.2 m a⁻¹ under the A1B scenario. There are significant differences between the projections from HadCM3 and ECHAM5. There are number of ways such differences can propagate through the regional ocean to the ice shelf. Higher air temperature over the southern ocean or higher ocean temperature at the FESOM northern boundary could cause universally higher temperatures in FESOM, causing increased shelf melting. Changes in prevailing wind over the southern ocean could potentially re-route water masses, affecting whether or not warmer water masses (such as CDW) come in contact with the ice shelf. Difference in sea ice extents and concentrations could greatly affect the air/sea exchange of heat and therefore water temperature. In the current study, the difference in the temperature anomaly simulated by FESOM forced by HadCM3 versus ECHAM5 results from the different surface heat flux simulated by each GCM (Timmermann and Hellmer, Reference Timmermann and Hellmer2013).

We carry out a total of seven experiments, summarized in Table 1. Two (Ctrl/m1 and Ctrl/m3) are forced for 200 a by applying the present day melt rate (3) with M _a = 0 to BISICLES model with the linear and nonlinear Weertman basal friction laws. Two (Had/m1 and Had/m3) come from applying the HadCM3 anomaly (4) to BISICLES model with the linear and nonlinear Weertman basal friction laws. A fifth (EC5/m1) applies the ECHAM5 anomaly to a BISICLES model with the linear basal friction law, a sixth (2 × Had/ml) doubles the melt rate of Had/m1, and the seventh (ka Had/m1) extends the Had/m1 simulation to t = ka.

Table 1.

Experiment settings used