Sensitivity of the Lambert-Amery glacial system to geothermal heat flux

Abstract

Geothermal heat flux (GHF) is one of the key thermal boundary conditions for ice-sheet models. We assess the sensitivity of the Lambert-Amery glacial system in East Antarctica to four different GHF datasets using a regional ice-sheet model. A control solution of the regional model is initialised by minimising the misfit to observations through an optimisation process. The Lambert-Amery glacial system simulation contains temperate ice up to 150 m thick and has an average basal melt of 1.3 mm a⁻¹, with maximum basal melting of 504 mm a⁻¹. The simulations which use a relatively high GHF compared to the control solution increase the volume and area of temperate ice, which causes higher surface velocities at higher elevations, which leads to the advance of the grounding line. The grounding line advance leads to changes in the local flow configuration, which dominates the changes within the glacial system. To investigate the difference in spatial patterns within the geothermal datasets, they were scaled to have the same median value. These scaled GHF simulations showed that the ice flow was most sensitive to the spatial variation in the underlying GHF near the ice divides and on the edges of the ice streams.

INTRODUCTION

Numerical ice-sheet models are an important tool for understanding the contribution of the cryosphere to past, present and future sea level rise. The temperature of the ice sheet is an important control on the flow rate of ice, influencing both the rate of internal deformation, basal melt and subsequent basal sliding. The thermal boundary conditions of the ice sheet are the GHF at the base of the ice sheet, the air surface temperature at the exposed surface of the ice and the ocean temperature beneath floating ice. An additional control on the thermal boundary conditions is the surface vertical velocity, which in pseudo steady state is the same as the accumulation rate. This governs how quickly the surface temperature is advected into the ice sheet. GHF influences the temperature of the ice and in part controls the conditions at the base of the ice sheet. The temperature of the ice is also controlled by the deformational heat generated from strain within the ice, the advection of heat due to ice motion, the conduction of heat through the ice sheet and frictional heating from basal sliding. The temperature of the ice is a control on the rheology of the ice and subsequently the rate of its deformation (Budd and others, 2013), with temperate ice (ice at melting point, which may contain a small fraction of liquid water) enhancing the flow significantly (Lliboutry and Duval, 1985). In addition, basal melt can lead to the lubrication of the till, lowering the resistance of the bed and leading to basal sliding (Pattyn, 2010). The performance of ice-sheet models in modelling ice temperature is difficult to evaluate as only spatially limited observations of in situ ice temperatures exist, with most being situated either at ice divides or on ice shelves, which represent two extremes of ice flow. Ice divides have near zero flow rates, limiting the contributions of deformational heating and basal frictional heat, while ice shelves are dependant on properties of the underlying ocean, with no contribution from GHF.

The Antarctic GHF is difficult to observe due to the ice sheet itself, which impedes access to the bed to measure the GHF directly with the exception of some isolated coastal ice free regions and deep ice core drilling sites (Fisher and others, 2015). Limited crustal heat production measurements (Carson and Pittard, 2012) are also used to estimate the GHF in localised regions (Carson and others, 2014). The magnitude of the GHF depends on spatially varying geological conditions such as the mantle heat flux, crustal thickness, sediment deposits and local radiogenic crustal heat production (RCHP) (Sandiford and McLaren, 2002; Fox Maule and others, 2005; Carson and others, 2014).

Early ice-sheet models input GHF as a constant (Hansen and Greve, 1996; Kerr and Huybrechts, 1999) or a regionally varying constant based on the origin of the crust (Pollard and others, 2005). Later developments provided spatially variable fields of GHF using inference based on magnetic fields (Fox Maule and others, 2005) and seismic models (Shapiro and Ritzwoller, 2004). Both of these techniques make assumptions about the local RCHP and acknowledge that they will not capture local small-scale variations in the GHF. The two datasets, Fox Maule and others (2005) and Shapiro and Ritzwoller (2004) are significantly different (Fig. 1), and are used in a variety of ice-sheet studies, mainly being used for basal inversions and boundary conditions (e.g. Pattyn, 2010; Martin and others, 2011; Larour and others, 2012; Sato and Greve, 2012; Morlighem and others, 2013; Golledge and others, 2015). The large differences between the two datasets often means studies utilise both datasets, or simply choose one without discussing the implications of using that particular dataset (Martin and others, 2011; Sato and Greve, 2012; Morlighem and others, 2013; Golledge and others, 2015). This may be because validation of these datasets is limited and therefore determining which may offer better performance within an ice-sheet model difficult. Understanding the influence that variations in GHF has on ice flow is not fully understood.

Fig. 1. The difference between seismic sourced GHF dataset (Shapiro and Ritzwoller, 2004) and a magnetic sourced GHF dataset (Fox Maule and others, 2005).

A number of studies have investigated the effect of differences in the GHF in Antarctica, and while local ice flow has been shown to be sensitive to variations in GHF (Pittard and others, 2016), the overall effect on ice volume has been found to be small, at least when compared with errors in other components of the model such as ice thickness (Pollard and others, 2005; Larour and others, 2012). Larour and others (2012) comment that slower flowing ice in the interior of the ice sheet will be more sensitive to the GHF, but that the frictional heating dominates GHF heating in regions of fast ice flow. This suggests that GHF in these regions are unimportant, although this contrasts with northern hemisphere ice sheets, with both the Greenland ice sheet (Rogozhina and others, 2012) and the Fennoscandian ice sheet (Näslund and others, 2005) found to be sensitive to the chosen GHF in numerical studies. Näslund and others (2005) found that regions with elevated GHF could lead to an increase in localised basal melt, which lubricated the bed and led to faster ice flow, which in turn increased the frictional heating, further increasing basal melt. The average basal melt expected under Antarctica ranges from 1 to 5.3 mm a⁻¹ (Llubes and others, 2006; Bell and others, 2007; Bell, 2008; Pattyn, 2010). While this basal melt rate is relatively low, under ice streams such as Pine Island Glacier, the basal rates could be as high as 600 mm a⁻¹ (Bell, 2008). An investigation by Hansen and Greve (1996) shows that as the GHF increases, it changes the basal properties of the ice sheet significantly, with higher GHF leading to greater areas of the base of the ice sheet reaching melting point, and the formation of temperate ice layers. Pattyn (2010) found that the GHF required to reach the melting point of ice is a function of ice thickness and surface temperature, with GHF as low as 40 m Wm⁻² at the base of the ice allowing for basal melting in low accumulation regions.

Part of the difficulty determining the importance of GHF is that the conditions at the base of the ice sheet are mostly unknown. Ice-sheet models may take different approaches to estimating basal properties and utilise the GHF differently. The GHF may be used as a thermal boundary condition to many ice-sheet models, with basal properties such as the strength of the till parametrised or estimated to allow for an evolving base of the ice sheet (Bueler and others, 2007; Bueler and Brown, 2009; Winkelmann and others, 2011). Another method of estimating the properties at the base of the ice is inverting for basal friction coefficients and/or viscosity of the ice (Morlighem and others, 2013; Gong and others, 2014) with GHF being used to generate a temperature profile through the ice to be used in the inversion.

The Lambert-Amery glacial system in East Antarctica (Fig. 2) has been observed to be relatively stable since 1968 (King and others, 2007, 2012; Yu and others, 2010; Wen and others, 2010; Pittard and others, 2015). The Lambert basin drains into the Amery Ice Shelf, which is long, relatively narrow and its grounding line is one of the deepest in Antarctica with an ice depth up to 2500 m (Fricker and others, 2002). The stability of this region allows models to be evaluated against a steady-state benchmark, which is difficult in many regions of Antarctica due to the localised rapid changes in ice dynamics. (Rignot and others, 2008; Pritchard and others, 2012; Shepherd and others, 2012). The Lambert-Amery glacial system has been difficult to model, with the location of the modelled grounding line advancing relative to observations, which leads to the over-estimation of the ice volume in the region (Martin and others, 2011; Golledge and others, 2012).

Fig. 2. The regional domain with the initial ftt_mask (green) and the final ftt_mask (blue) indicated. The ice shelf extent from bedmap2 is indicated in black. Inset: location of the Lambert-Amery glacial system within Antarctica, showing the square region (black) that encompasses the regional model.

This study utilises a regional ice-sheet model of the Lambert-Amery glacial system used to investigate the influence of both the magnitude and variability of GHF on the ice sheet. We first outline the regional domain and then detail the optimisation of a number of ice-sheet model parameters until a steady-state solution which approximately estimates the current configuration of the region is found. We then test the sensitivity of the optimisation to different GHF datasets. Finally, we scale each dataset to the control GHF to assess the impact of spatial differences between the various datasets on the thermal regime of the region.

MODEL DESCRIPTION

Ice-sheet model

The ice-sheet model utilised by this study is the Parallel Ice Sheet Model (PISM) version 0.6.2 (Bueler and others, 2007; Bueler and Brown, 2009; Winkelmann and others, 2011). PISM is a 3-D thermodynamically coupled model with a shallow ice approximation (SIA)/shallow shelf approximation (SSA) hybrid scheme that utilises a structured finite difference discretization. The SIA approximates ice flow for grounded ice where vertical shear dominates and the SSA approximates ice flow for floating ice where horizontal shear dominates. In grounded regions that are sliding, part of the ice flow is calculated from the SSA and part from the SIA (Bueler and Brown, 2009). PISM utilises an enthalpy scheme for its thermal model component and is both mass and energy conserving (Aschwanden and others, 2012). The calving law options utilised by this study is an eigen calving law (Levermann and others, 2012) combined with a minimum thickness calving law. When the principal strain rate of a region of ice shelf exceeds a threshold set (eigen_calving_k), or the region drops below a set ice thickness (thickness_calving_threshold), the region will calve. The regional model used herein is described by Pittard and others (2016).

Regional domain

The Lambert-Amery glacial system is identified through the PISM drainage basin delineation tool (included in the PISM regional-tools), which determines the drainage basin by using the gradient of the surface elevation to determine the maximum source point of ice from a terminus specified by the user. The calculated basin was enlarged slightly to capture the ice divides more accurately. The drainage basin outline is shown in Figure 2, and within this basin the full PISM model applies. The region outside the basin has an adaptive surface mass-balance mechanism (using PISMO executable and the force to thickness mechanism), which forces the ice thickness within this region to match the initial ice thickness, ensuring that the boundary conditions at the edge of the domain, and the region outside the drainage basin of interest will minimally impact the solution within the domain itself (See Supplementary Information for full details).

Input datasets

The regional model requires a set of boundary, initial and forcing conditions. The boundary conditions of the ice sheet are the bed topography and GHF. The initial condition is the ice thickness. The bedrock topography and initial ice thickness for the regional model is given by a modified bedmap2 dataset (Pittard and others, 2016). The GHF dataset used was created by using the Fox Maule and others (2005) methodology on an updated magnetic field model 7 (http://websrv.cs.umt.edu/isis/index.php/Antarctica_Basal_Heat_Flux).

The forcing conditions used by the regional model are surface mass balance, ice surface temperature, oceanic forced basal melt rate and ocean temperature. The surface mass-balance and ice-surface temperature are from the 1979–2013 ANT27\2 RACMO2.3 (van Wessem and others, 2014) dataset, with minor modifications over nunataks to reduce ice growth on these regions (See Supplementary Information for full details). The ice shelf basal mass balance is controlled by a parametrisation for oceanic basal melt rates with ocean temperature held at a constant 271.45 K. The oceanic basal melt rate parametrisation is calculated following the modifications outlined in Pittard and others (2016). The initial oceanic melt rate is approximately the same as that in Galton-Fenzi and others (2012) with an average melt of 0.8 m a⁻¹ from our parametrisation compared with 0.78 m a⁻¹ from the ocean cavity model (See Supplementary Information for full details).

Model parameters

The PISM ice-sheet model is controlled by a range of input parameters (See PISM Users Manual, date accessed 17 June 2014). The regional model used within this study uses a horizontal resolution of 5 km, with the exception of thermal only simulations, which are simulated at 10 km horizontal resolution. The vertical resolution is 15 m for all simulations. The domain edge boundary conditions are derived from a low resolution full Antarctic domain model (see Supplementary Information). The PISM model variables bmelt, tillwat, enthalpy and velocity for the dirichlet velocity boundary, u_ssa_bc,v_ssa_bc, are applied in a 10 km strip outside the domain (See Supplementary Information for full details). Six of the model input parameters were chosen to vary through an optimisation process to create a realistic simulation of the Lambert-Amery glacial system. All other input parameters were held at the defaults (See PISM Users Manual) (See Supplementary Information for full details).

The four variables, which the model was firstly optimised for were the shallow ice approximation enhancement factor (sia_e), the shallow shelf approximation enhancement factor (ssa_e), the quotient in the pseudo plastic sliding law (pseudo_plastic_q, (See PISM Users Manual) and the parametrisation of till strength (topg_to_phi). These variables were chosen guided by the previous experiments of Martin and others (2011); Golledge and others (2015) and initial experiments testing the relative importance of each variable. The final two variables which are optimised vary the calving front location and calving rate within the model, with the threshold of the principal strain rate for eigen calving (eigen_calving_k) (Levermann and others, 2012) and threshold where the ice shelf is considered too thin to be realistic and is automatically calved (thickness_calving_threshold). The primary criteria for a stable solution was the grounding line being situated on the same topographic sill as observations, with secondary criteria being how accurately the simulated ice thickness and velocities matched observations. This secondary criteria was assessed by comparing the misfit between the observed and simulated ice sheet for ice thickness and surface velocity, calculating the mean and standard deviation of both the simulated and observed ice sheet, and finally computing the root-mean-square error between the two values. Each of the variables were iteratively varied and assessed using the two criteria, until a final set of variables, which most accurately matched observations were found.

The final parameters from this optimisation process were sia_e = 1.8, ssa_e = 1.6, pseudo_plastic_q = 0.5 and topg_to_phi = 10,30,−1500,−500, eigen_calving_k = 1.9e15 and thickness_calving_threshold = 225. The sia_e is lower than expected from laboratory experiments, which indicate that the enhancement due to anisotropy should be at least 3 on average across the domain (Budd and Jacka, 1989; Budd and others, 2013). This indicates that there is another factor, which is being convolved into the sia_e optimisation, such as basal resistance. PISM utilises a parametrisation, which aims to reflect the reduction in flow due to the roughness of the bed topography (See PISM Users Manual). This parametrisation is limited by the interpolation and smoothing of the bed topography datasets, which causes the reduction in flow as the bed is relatively rough in regions with high-resolution data and relatively smooth where the topography is under-sampled. The final regional model solution was simulated for a 200 000 years thermal simulation (-no_mass turned on which holds the ice thickness constant and evolves only the thermal ice sheet), followed by a 45 000 years simulation (henceforth the control solution).

METHODS

Geothermal flux datasets

Three geothermal flux databases were chosen to investigate their influence on the ice configuration of the regional domain compared with the GHF chosen in our regional domain (Fig. 3). The GHF used in the control solution utilises the Fox Maule and others (2005) methodology, but using an updated magnetic field model, FM7 (http://geomag.org/models/MF7.html) (henceforth fm_2012). The three databases, which we will compare our control solution to are Fox Maule and others (2005) (fm_2005), Shapiro and Ritzwoller (2004) (sr_2004) and An and others (2015) (an_2015), as shown in Figure 3. The fm_2005 and sr_2004 datasets were accessed through the ALBMAP compilation (Le Brocq and others, 2010). The sr_2004 and an_2015 datasets utilise a seismic model, while fm_2005 and fm_2012 are based on magnetic models.

Fig. 3. GHF over the domain from (a) fm_2012 (b) fm_2005 (c) sr_2004 (d) an_2015 (e) fm_median (f) fm_scaled (g) sr_scaled (h) an_scaled. Ice shelf mask from bedmap2 shown in red.

The fm_2005 dataset shows spatial similarities to the fm_2012 dataset, with similar features including a relatively higher GHF beneath the Lambert Glacier, the elevated region north west of the ice shelf and the relatively cold region beneath the Fisher Glacier. Overall, fm_2005 has a much higher GHF, with a median GHF of 59.1 m Wm⁻² compared with just 40.8 m Wm⁻² for the fm_2012 dataset. The sr_2004 dataset has significantly less spatial details than the other datasets, with a very small gradient in GHF from south-east to north-west. The median GHF in sr_2004 is 52.6 m Wm⁻². The an_2015 dataset shares similar features to the magnetic field datasets, with a higher GHF in the north west, but in contrast to the magnetic datasets the region beneath the Lambert Glacier is relatively cooler than the background field. The median GHF in an_2015 is 53.9 m Wm⁻².

To test the differences of the spatial variability of the GHF datasets without being influenced by the elevated GHF, four additional GHF datasets are constructed. The first dataset constructed was created by using the median of the fm_2012 dataset, 40.8 m Wm⁻², as a constant region wide value (labelled as fm_median, Fig. 3e). The three other datasets are created by scaling the fm_2005, sr_2004 and an_2015 datasets by the median of fm_2012 datasets. The fm_scaled, sr_scaled and an_2015 datasets were scaled by multiplying the GHF dataset by the median of fm_2012 divided by each respective dataset (40.8/59.1, 40.8/52.1, 40.8/53.9), forcing the median of each dataset to match that of the fm_2012 dataset. These datasets are labelled fm_scaled, sr_scaled and an_scaled respectively (Figs 3f–h). The median was chosen over the mean as there is a region of high GHF in the north eastern corner of the fm_2012 dataset, which skewed the mean to a much higher value of 49.1 m Wm⁻², which would have caused the constructed datasets to have an elevated GHF relative to fm_2012 in the regions beneath the major active glaciers.

Experimental design

Each of the eight different GHF datasets are used in experimental model runs (summarised in Table 1) with 10 km horizontal resolution, a constant ice thickness (-no_mass), simulated until the enthalpy is close to thermal quilibrium for the given ice thickness of the control solution. This step was conducted on the control solution GHF (fm_2012) as well, to measure any lingering transient thermal effects from the control solution. Following the thermal equilibrium runs, each GHF dataset is run at 5 km horizontal resolution for another 2000 years yielding a pseudo-steady-state solution for ice thickness. However, any significant grounding line migration or changes in surface velocities should be evident over this time period.

Table 1. List of experimental runs

CONTROL SOLUTION

Comparison with observations

The control solution (Fig. 4) meets our primary goal of a stable grounding line position along the same topographic sill as the observed grounding line. The velocities of the control solution are characterised by faster ice flow through the main trunks of the glaciers compared with observation and slower velocities adjacent to the main glaciers. This characteristic could be caused by the satellite footprints, which are lower-resolution than the numerical model. The ice thickness of the control solution is thicker in the drainage basin of the Fisher and Charybdis Glaciers than observations. These glaciers flow through narrow channels between nunataks, which will likely require higher horizontal resolution to better resolve. Conversely, the Lambert Glacier, and to a lesser extent the Mellor Glacier, show ice thickness that is thinner than observations. This is likely due to the deep topographic troughs that exist within these basins, which will lead to the topg_to_phi parametrisation enforcing a very weak till in these regions and allowing faster flow and easier sliding. The control solution's grounding line has advanced slightly compared with the observations, however, given the uncertainty in the bedrock elevation and the modifed bedmap2 dataset used it is considered an accurate representation. The calving front of the control solution is further north along the western edge of the embayment and further south to the eastern edge. The northward position of the western ice front could be due to the lack of a pinning point in the topography, which restricts and shifts the flow to the eastern edge in observations. The surface velocities are slower towards the deep grounding line of the AIS, but slightly faster towards the middle of the AIS before slowing towards the calving front. The slower velocities at the deep grounding line could be due to horizontal resolution of the model, or potentially over-buttressing from the side wall drag. The model reaches close to observed ice thickness and velocities through the centre of the ice shelf, which is due to the thickness based ice shelf parametrisation we apply. As the ice flow is restricted, the ice gets thicker at the deep grounding line leading to the basal melt rate increasing. The ice thickness at the calving front is very similar to the observed ice thickness, however, with the lower velocities relatively less ice mass is being calved within the model. Overall, this means that while the combination of calving and basal melt from the ice shelf preserves a similar ice thickness within our control solution, we are preferentially losing more ice loss from basal melting than what would be expected by combining the MEaSUREs velocities (Rignot and others, 2011) and the bedmap2 (Fretwell and others, 2013) ice thickness at the calving front.

Fig. 4. (a) The MEaSUREs surface velocities (Rignot and others, 2011). L = Lambert Glacier, M = Mellor Glacier, F = Fisher Glacier, C = Charybdis Glacier. (b) The bedmap2 ice thickness (Fretwell and others, 2013). (c) The difference between the control solution and the MEaSUREs velocities. (d) The difference between the control solution and the bedmap2 ice thickness. (e) The percentage difference between the control solution and the MEaSUREs velocities (Rignot and others, 2011). (f) The percentage difference between the control solution and the bedmap2 ice thickness. The bedmap2 ice shelf and coastline is outlined in black, the control solution's ice shelf and coastline is shown in green.

Thermal regime of the Lambert-Amery glacial system

The thermal regime of the Lambert-Amery glacial system is displayed in Figure 5. The average ice hardness (describes viscosity) varies from 5 × 10⁸ Pa s^1/3 towards the ice divides, softening to 2 × 10⁸ Pa s^1/3 towards the grounding line with the exception of the three major ice streams. The Mellor and Lambert Glaciers are harder than the surrounding ice, with the faster velocities from basal temperate ice and basal sliding allowing for the transport of harder ice towards the grounding line. The regions that are at the basal melting point temperature are primarily the major ice streams and the region with significantly elevated GHF north of the Charybdis Glacier. The spatial distribution of the temperate ice matches the distribution of the regions at the basal melting point temperature. The temperate ice thickness varies from up to 150 m in the Lambert Glacier, 100 m in the Mellor Glacier and as thin as 30 m in the Fisher Glacier. These layers are likely formed due to the relatively small width of the ice streams compared with the drainage region, with compression leading to higher internal heating and thicker layers of temperate ice. The saturated till shows the regions which will slide relatively easily, with an average basal melt rate within the area of the saturated till of 17.4 mm a⁻¹. This converts to an entire glacial system basal melt rate of 1.3 mm a⁻¹. The maximum basal melt rate is 504 mm a⁻¹. The basal melt rates fall within the expected range for ice sheets, with the range of basal melt rates calculated by previous models ranging from 1 to 5.3 mm a⁻¹. The maximum basal melt rate is relatively high, but is lower than the expected 600 mm a⁻¹ estimated for the faster flowing Pine Island Glacier. The thermal regime described falls with the bounds for our limited knowledge of the thermal properties of the ice-sheet base.

Fig. 5. Thermal properties of the control solution: (a) the average ice hardness. (b) The temperature of the ice at the base of the ice sheet. (c) The thickness of the temperate ice layer at the base of the ice sheet. (d) The basal melt rate of the grounded ice sheet with the region within the red contour indicating the extent of the saturated till. The green line indicates the control solution's grounding line and the black the observed grounding line.

ICE-SHEET RESPONSE TO GHF

Thermal Equilibrium Experiments

Each GHF dataset was run until thermal equilibrium was reached. The three GHF datasets with a higher average GHF required 400 000 a to reach equilibrium. The overall change in enthalpy was <1% for all simulations, with the enthalpy of each simulation ~2 × 10²³ Joules. Even though this change is small, it is important as the spatial distribution varies. The volume and area of temperate ice increases in thermal_fm_2005, thermal_sr_2004 and thermal_an_2015 simulations (Table 2). The fm_2005 GHF dataset has the highest average heat flux, which leads to the simulation having twice as much temperate ice as the thermal_control in addition to an increase in over 50% in area. The thermal_sr_2004 and thermal_an_2015 simulations demonstrate that the datasets lead to different spatial distributions, with thermal_an_2015 showing more temperate ice in a smaller region relative to thermal_sr_2004. The area of temperate ice also represents the area within the ice streams at basal melting point. The scaled datasets show relatively small amount of change, with the control solution having 10–20% more temperate ice. The basal melt rates (bmelt) of the hotter GHF datasets is higher than the control and scaled datasets, however, most are lower than the control solution, which demonstrates the influence of the horizontal resolution change. Analysis of changes in the main ice streams is limited by the lower horizontal resolution, which does not resolve the Mellor or Fisher Glaciers sufficiently to create a smooth surface velocity. This impacts the formation of temperate ice within the ice streams, with an uneven layer forming based on the variable surface velocities.

Table 2. The thermal properties of the ice sheet across all thermal simulations

Comparing the GHF datasets

The three simulations using the original GHF datasets (exp_fm_2005, exp_sr_2004 and exp_an_2015) had higher volumes of temperate ice compared with the control experiment, which indicates the velocities in the experimental run will be higher. The initially higher velocities leads to the transfer of mass from interior of the ice sheet towards the coast, consequently causing the grounding line to advance (Fig. 6). Once the grounding line advances, the flow configuration of the region changes. The observed flow configuration, which is simulated in the control solution, shows three tributary glaciers flowing into the Amery Ice Shelf, each with its own grounding line. As the grounding line advances, the three glaciers now converge into a single outlet glacier, which has a smaller area of outflow, but is significantly deeper. The upstream velocities reduce as the thickness around the former grounding line increases. These changes make it difficult to assess the impact of the three geothermal datasets as the primary changes in the thermal regime are linked to the changes in flow configuration rather than the change in GHF. While this confirms that the ice sheet is sensitive to elevated GHF, the model was optimised for a GHF with a lower overall magnitude. Further optimisation would stabilise the grounding line with any GHF, however, this is computationally expensive and hence the scaled datasets will be used to assess the importance of the spatial variability of the GHF.

Fig. 6. (a) The surface elevation of the thermal_control simulation. The difference in ice thickness between the exp_control and (b) exp_fm_2005, (c) exp_sr_2004 and (d) exp_an_2015. (e) The surface velocity of the exp_control simulation. The difference in surface velocity between the exp_control and (f) exp_fm_2005, (g) exp_sr_2004 and (h) exp_an_2015.

Comparing the scaled GHF datasets

The thermal properties of the scaled dataset experiments are relatively similar, with the control maintaining the slightly elevated volume and area of temperate ice seen in the thermal experiments (Table 3). The average melt rates vary only slightly, although the exp_control has slightly lower average melt rate within the regions where the till is saturated. More importantly, the average melt rates have increased to be similar to that of the control solution, and now falls within the expected range for the Antarctic Ice Sheet.

Table 3. The thermal properties of the ice sheet across all scaled experimental simulations

The scaled datasets maintain a grounding line on the same topographic slope as the control (Fig. 7), which enables the comparison of the relative changes in GHF. There are some clear differences in the surface elevation and surface velocities between the four datasets. The exp_fm_median simulation has faster velocities through the Mellor and Fisher Glaciers than the Lambert Glacier compared to the control solution, which corresponds with the thickness change in these regions. The simulation using the older magnetic derived dataset, exp_fm_scaled, shows significant differences to the control solution using the updated magnetic derived dataset. All three major ice streams are relatively slower relative to the control solution, with the main increases in velocity and thickness being outside these main ice streams. The two simulations using the GHF derived from seismic methods, exp_sr_scaled and exp_an_scaled, have clear similarities, with both showing increases in velocities within the Fisher and one arm of the Mellor Glacier, while the Lambert Glacier is considerably slower. This leads to changes in the ice thickness. One characteristic which holds across all scaled datasets is that the region north of the Charybdis Glacier is affected by the relatively high GHF in the fm_2012 dataset.

Fig. 7. The difference in ice thickness between the exp_control and (a) exp_fm_median, (b) exp_fm_scaled, (c) exp_sr_scaled and (d) exp_an_scaled. The difference in surface velocity between the exp_control and (e) exp_fm_median, f) exp_fm_scaled, (g) exp_sr_scaled and (h) exp_an_scaled. The control solution grounding line is shown in black and the scaled datasets in green.

To assess the importance of these changes, and where the changes are directly related to the geothermal change, the ratio of change in surface velocities (Figs 7e–h) to the underlying GHF change (Figs 8a–d) is derived (Figs 8e–h). A positive ratio indicates that a higher relative GHF leads to an increase in velocity or a lower GHF leads to a decrease in velocity. We assume that a negative ratio indicates that there is no link between the underlying change in GHF and the surface velocity.

Fig. 8. The relative difference in GHF between fm_2012 and (a) fm_median, (b) fm_scaled, (c) sr_scaled and (d) an_scaled. The ratio between the change in velocity (Figs 7e–h) and the relative change in GHF between exp_control and (e) exp_fm_median, (f) exp_fm_scaled, (g) exp_sr_scaled and (h) exp_an_scaled. The control solution grounding line is shown in black and the scaled datasets in green.

Three main patterns can be discerned from observing the ratios. The first is that the ice divides are sensitive to changes in local GHF, with strong ratios linked across all four scaled simulations. The second is that the ice streams are insensitive to change in the underlying GHF. The final pattern is that the edges of the ice streams are sensitive to changes, with their widths and upstream extent both being modified by the underlying GHF. This could be important to correctly modelling climate feedbacks, as regions which are cold relative to warmer adjacent ice streams may resist the acceleration driven by changes in bordering ice shelves. Conversely, if the regions adjacent to the ice streams are already close to pressure melting point and therefore sliding due to enhanced lubrication at the base from meltwater, the system could respond far more rapidly. There also appears to be evidence of ice piracy in exp_an_scaled, with the change in velocities between the two arms of the inland extent of the Mellor Glacier linked to an increase in the GHF in one arm, which causes the decrease in velocity within the other arm, regardless of the underlying GHF in that region.

Assessing the GHF in terms of which dataset is most realistic is difficult based on this study, as the optimisation process creates a bias towards the fm_2012 dataset. It is evident that higher GHF leads to faster velocities, however, under an optimisation process the solutions could end relatively close to the same in all cases. The relative differences could be important, for example, our control solution had a relatively thick Fisher Glacier and a relatively thin Lambert Glacier compared with observations. The exp_an_scaled dataset would alleviate this discrepancy as it preferentially leads to a thicker Lambert Glacier and a thinner Fisher Glacier, but would probably lead to the Mellor Glacier also becoming too thick. The regional variations between each dataset likely will lead to each dataset excelling in different regions, and these relative differences could be used to choose an ideal GHF dataset for each region.

CONCLUSIONS

We present a realistic simulation of the Lambert-Amery glacial system, with the grounding line and calving front accurate with observations. The thermal regime of the control solution shows temperate ice layers up to 150 m thick in the Lambert Glacier and up to 100 m thick in the Mellor. The control solution's glacial system wide average melt rate is 1.3 mm a⁻¹, with a maximum basal melt rate of 504 mm a⁻¹. These numbers are consistent with previous modelling estimates from Antarctica. Three different GHF datasets are compared with the control dataset, with higher GHF leading to the formation of significantly more temperate ice. The increase in temperate ice leads to faster surface velocities, which causes an advance in the grounding line, changing the flow configuration of the region. To compare the spatial differences a set of scaled geothermal heat geothermal heat flux datasets were created. The regions which were most sensitive to changes in the underlying GHF were near ice divides and adjacent to the ice streams. The ice streams themselves were relatively insensitive to changes. Future studies should consider a robust evaluation of the effects of choosing one GHF dataset over another on the solution. Further direct observation of GHF are needed to evaluate the remote sensing derived products available to modellers, and additional techniques are needed to quantify GHF and impact on the ice-sheet flow.