Basal conditions of Denman Glacier from glacier hydrology and ice dynamics modeling
Basal sliding in Antarctic glaciers is often modeled using a friction law that relates basal shear stresses to the effective pressure. As few ice sheet models are dynamically coupled to subglacial hydrology models, variability in subglacial hydrology associated with the effective pressure is often implicitly captured in the basal friction coefficient – an unknown parameter in the basal friction law. We investigate the impact of using effective pressures calculated from the Glacier Drainage System (GlaDS) model on basal friction coefficients calculated using inverse methods in the Ice-sheet and Sea-level System Model (ISSM) at Denman Glacier, East Antarctica, for the Schoof and Budd friction laws. For the Schoof friction law, a positive correlation emerges between the GlaDS effective pressure and basal friction coefficient in regions of fast ice flow. Using GlaDS effective pressures generally leads to smoother basal friction coefficients and basal shear stresses, and larger differences between the simulated and observed ice surface velocities compared with using an effective pressure equal to the ice overburden pressure plus the gravitational potential energy of the water. Compared with the Budd friction law, the Schoof friction law offers improved capabilities in capturing the spatial variations associated with known physics of the subglacial hydrology. Our results indicate that ice sheet model representation of basal sliding is more realistic when using direct outputs from a subglacial hydrology model, demonstrating the importance of coupling between ice sheet and subglacial hydrological systems. However, using our outputs we have also developed an empirical parameterization of effective pressure that improves the application of the Schoof friction law without requiring explicit hydrological modeling.
The health of Antarctic glaciers and their future susceptibility to climate-driven change is often assessed by the retreat rates of their grounding lines and melt rates of adjacent ice shelves. In the East Antarctic, Denman Glacier of the Denman–Scott catchment (Fig. 1) has seen some of the fastest grounding line retreat of the last 20 years, losing 5.4 km since 1996 (Brancato et al., 2020), and the highest ice shelf melt rates of 116 m a−1 from 2010–2018 (Adusumilli et al., 2020). Perched on a grounding line above a deep trough that descends 3.66 km below sea level (Morlighem, 2020), Denman Glacier is potentially at risk of marine ice sheet instability and rapid retreat, which is of significant concern given that the glacier drains a region containing 1.5 m of sea level equivalent (Brancato et al., 2020).
Determining the future of Denman Glacier, and others in the Antarctic, largely relies on ice dynamics models to capture the evolution of ice flow under changing climates. These models frequently use inversion techniques to constrain important parameters that control ice velocity, including those related to the basal environment beneath the ice sheet (Morlighem et al., 2013). One such parameter is the basal friction coefficient, which is a key component of friction laws including the Weertman (1957), Budd et al. (1979) and regularized Coulomb friction laws (Schoof, 2005; Gagliardini et al., 2007). Friction laws propose a functional dependency of the basal shear stress on the basal sliding velocity and in the case of the Budd and regularized Coulomb friction laws, the effective pressure , which is defined as the ice overburden pressure (pi) minus the subglacial water pressure (pw). However, if the basal shear stress does not actually have the functional dependence on the basal sliding velocity and effective pressure proposed in the basal friction law or if it has a functional dependence on other properties of the bed such as roughness, substrate, or temperature, then this will be implicitly captured in the basal friction coefficient. Therefore, a spatially variable basal friction coefficient suggests a friction law which either fails to capture the proper functional dependency on basal sliding velocity and effective pressure or omits the dependency of the basal shear stress on other quantities. By consequence, a basal friction coefficient that is both smooth – has little local variability – and has limited domain-wide trends is desirable. In the search for ways of describing basal friction, Beaud et al. (2022) proposed a generalized friction law to be used on both rigid and deformable beds. They suggest that the unknown parameters in their friction law can be found using ice surface velocity time series and inversion techniques on a variety of glaciers implemented at the global scale.
Antarctic subglacial hydrology is increasingly being shown to be varied and dynamic, with large channels that discharge into ice shelf cavities (Dow et al., 2020; Indrigo et al., 2021) and effective pressures that vary spatially both above and below zero (i.e., when subglacial water pressure exceeds the ice overburden pressure). Modeling of subglacial hydrology has shown close links between regions of low effective pressure and high ice surface velocity (Dow et al., 2018, 2020, 2022); this control of basal water pressure on sliding rates is also well known for Greenland and Alpine glacier systems (Nienow et al., 2017; Iken and Bindschadler, 1986).
Despite the key role of basal boundary conditions, and in particular subglacial water pressure, in ice dynamics, there has not yet been a systematic investigation of the impact of the effective pressure on basal sliding in different friction laws applied in the Denman–Scott catchment. Given the critical role of the friction law for accurate ice dynamics and sea level predictions in the future (e.g., Åkesson et al., 2021; Yu et al., 2018), this is an important missing link. Here, we address this by using the Glacier Drainage System (GlaDS; Werder et al., 2013) model to characterize the effective pressure and subglacial hydrology of Denman Glacier. We use these effective pressures as inputs in the Budd and Schoof friction laws implemented in the Ice-sheet and Sea-level System Model (ISSM; Larour et al., 2012) and calculate the resulting basal friction coefficients using inversion.
2.1 GlaDS setup
GlaDS is a 2D finite-element subglacial hydrology model that calculates the development of a distributed water system on the elements and channel growth, fed by the distributed system on the element edges (Werder et al., 2013). The model runs are initiated with no channels, and therefore channel formation is dictated entirely by the evolving pressure and flow conditions across the model elements. GlaDS model equations are described in full in Werder et al. (2013), and the model's application to other Antarctic systems is covered in detail in Dow et al. (2018, 2020, 2022). The model parameters used in our simulations are summarized in Table 1. The bedrock bump height, cavity spacing, and sheet conductivity act together in the distributed system equations to either constrict or open the system to water flow as show in
where is the rheological constant of the ice multiplied by an order 1 factor depending on the cavity geometry, h is the hydrology sheet thickness, N is the effective pressure, n is the exponent in Glen's flow law, he is the englacial storage, t is time, η is a prescribed source term shown in Fig. 1c, q is the hydrology sheet discharge, k is the sheet conductivity, ϕ is the hydraulic potential, and and are exponents in the Darcy–Weisbach law describing fully turbulent flow. w is the hydrology sheet opening rate equal to when h<hr and zero otherwise; here ub is the basal ice speed, hr is the typical bedrock bump height, and lr is the typical cavity spacing.
As discussed in Dow (2023), when the system is over-constricted (i.e., it is difficult for water to flow through the hydrologic system), the pressures are unrealistically high – much of the domain is above ice overburden pressure – and the model fails to converge. When the system is under-constricted (i.e., water flows through the hydrologic system with ease), the pressures are well below ice overburden pressure for much of the domain (much of the domain is below 50 % of ice overburden pressure), which is unrealistically low for steady-state Antarctic systems that are not driven by surface water input. The variables controlling the constriction of the hydrologic system are k, hr, and η, with a more constricted system arising from larger η and smaller k and hr. We test order of magnitude changes in k to determine a suitable level of constriction of the system. While there is some variation within the range of acceptable pressures, the output we present is the median and therefore is the most appropriate for representing the hydrologic pressure in ice sheet dynamics equations without further information from in situ measurements for example. Future work with a full coupling of hydrology and ice dynamics can explore sensitivity to different distributed system inputs. Channel conductivity is, similarly, a median value applied in GlaDS in Antarctica (Dow et al., 2020). The ice flow constant is set for an average ice column temperature of −10 ∘C.
The domain for the GlaDS model is a subset of the grounded portion of the Denman–Scott catchment based on the hydraulic potential (at overburden) catchment for this region (Fig. 1). The mesh includes 17 700 nodes and the domain-wide mean edge length is 3.1 km. However, we refine the mesh to a minimum edge length of 200 m in regions where basal sliding velocities are greater than 40 m a−1 and along the Denman Glacier grounding line. The basal and surface topography used in the GlaDS run is from BedMachine version 2 (Morlighem, 2020). Basal water and sliding velocity inputs are taken from the JPL_ISSM ISMIP6 Antarctic control run final time step (Seroussi et al., 2019), which was a thermal steady-state simulation using the enthalpy formulation implemented in ISSM (Seroussi et al., 2013) and the Blatter–Pattyn (BP) approximation to the full Stokes equations (Blatter, 1995; Pattyn, 2003). This model solved for ice viscosity of the floating ice shelf and the basal friction coefficient using data assimilation techniques (MacAyeal, 1993; Morlighem et al., 2010), which differs from our application of ISSM to assess basal boundary conditions. However, lacking an alternative starting point for basal sliding and velocity, this is the best available option prior to full coupling between hydrology and ice dynamics in ISSM. The sliding velocity acts to open up distributed system cavities and, at velocities greater than 800 m a−1, can cause model instabilities and so is capped at this value. Tests of similar caps for model runs of winter conditions at Helheim Glacier (Poinar et al., 2019) demonstrate it has little impact on the effective pressure. We assume a temperate bed throughout, although regions with zero water input (in the interior or southernmost-regions of the domain) are not an active part of the hydrological system. The model is run for 10 000 d to near steady state (there are changes in the water sheet thickness in deep pockets on the order of millimeters per day), providing outputs including channel size and discharge, distributed system discharge, water depth, and effective pressure.
2.2 ISSM setup
ISSM is a finite-element model that uses a non-uniform mesh to simulate ice dynamics. We employ the inverse capabilities within ISSM to estimate basal friction coefficients in the Denman–Scott catchment using various basal friction laws, as described below. The inverse model uses the shallow-shelf approximation (SSA; MacAyeal, 1989; Morland, 1987) to the full Stokes equations, described in Eq. (3) with the depth-averaged effective viscosity μ, which is given in Eq. (4). The SSA is described in full in Larour et al. (2012).
where u and v are the x and y components of ice velocity, respectively, H is the ice thickness, ρi is the density of ice, g is the absolute value of the gravitational acceleration, S is the surface elevation, is the effective strain rate tensor, n is the exponent in Glen's flow law, and is the ice rigidity.
The domain for the ISSM model is the full Denman–Scott catchment, extending the boundary used for the GlaDS simulations to include the floating ice shelves (Fig. 1). The ISSM mesh is comprised of 66 518 nodes, with non-uniform mesh refinement for faster flowing ice using the MEaSUREs v2 ice surface speed from (Rignot et al., 2011, 2017; Mouginot et al., 2012, 2017), described in Appendix A2. The minimum mesh edge length is 500 m. We note that the GlaDS mesh used for the subglacial hydrology simulations differs from that used in the ISSM simulations for the ice dynamics as a result of differences in ice and hydrology drainage catchments; the impact of this is discussed in Appendix A1. We test the sensitivity of the basal friction inversion simulations to various different meshes, finding that the results are less sensitive to the effect of the mesh configuration than to the specification of effective pressure and basal friction law used (Appendix A2). The bed topography and surface elevation are from BedMachine v2 (Morlighem, 2020). The ice rigidity is calculated using inverse methods (Appendix B).
Solving for basal friction coefficients
We investigate the impact of the prescription of the effective pressure in the Budd and Schoof friction laws by calculating basal friction coefficients using inversion. The scalar form of the Budd friction law is given by the following expression:
where τb (Pa) is the basal shear stress, α (s m) is the basal friction coefficient, N (Pa) is the effective pressure, ub (m a−1) is the basal sliding speed, and m=1 is a power law exponent taken to be linear as in Åkesson et al. (2021), Åkesson et al. (2022), Yu et al. (2018), Choi et al. (2021), and Baldacchino et al. (2022). The scalar form of the Schoof friction law can be written as
where C (kg m s) is the friction coefficient, m is the power law exponent (here given by ), and Cmax (m s) is Iken's bound, which is the effective cap on the basal shear stress (Appendix A5; Iken, 1981).
For each friction law, we calculate an inversion for the basal friction coefficients by minimizing a cost function that includes the contributions of absolute and logarithmic misfits between the observed and simulated ice surface speeds (Appendix B; Morlighem et al., 2013). Weighting parameters for the absolute and logarithmic misfit terms for the various friction laws are shown in Appendix B, Table B1. We also add a Tikhonov regularization term to the cost function, which regulates the local variability in the basal friction coefficient by penalizing large gradients in the basal friction coefficient, making the basal friction coefficient smoother. The weighting parameter of the Tikhonov regularization term is determined using an L-curve analysis (Appendix B; Hansen, 2000), which allows us to find a basal friction coefficient with small local variability while maintaining a good absolute and logarithmic velocity misfit. The final weighting parameters for the Tikhonov regularization term for each friction law are reported in Table B1.
Though it is not the primary focus of this work, we invert for ice rigidity as well, while initializing our model for the inversion for the basal friction coefficients. By inverting for the ice rigidity we capture ice rheological processes which are not explicitly accounted for in the model such as damage, anisotropy, chemical impurities, and liquid water.
We perform the following inversion procedure. First, we invert for the ice rigidity over the floating portion of the domain. Next, we invert for the Budd basal friction coefficient over the grounded portion of the domain, using an ice rigidity on grounded ice specified by the Paterson function from Cuffey and Paterson (2010) and surface temperatures from RACMOv2.3 (van Wessem et al., 2018). After the Budd inversion we invert for the ice rigidity over the entire domain. We next use the basal friction coefficient estimated using the Budd friction law to compute an initial estimate of the basal friction coefficient for the Schoof friction law. We perform inversions for the basal friction coefficients of the Budd and Schoof friction laws with the ice rigidity from the inversion prior; these are the main simulations discussed in the text that follows (it is worth noting that the Budd friction coefficient converges to the result of the Budd friction coefficient inversion). We perform a final rigidity inversion over the entire domain. The cost functions to be minimized for each inversion are described in detail in Appendix B.
We compare the difference in the basal friction coefficients when we use two different prescriptions for the effective pressure. (1) The first is an effective pressure given by assuming water pressure equals the ice overburden pressure plus the gravitational potential energy of the water: . Here ρi is the density of ice, ρw is the density of seawater, g is the absolute value of the gravitational acceleration, H is the ice thickness, and B is the bed elevation which takes negative values below sea level. (2) The second is the effective pressure taken directly from the GlaDS simulations, which we refer to as NG. We cap the effective pressure at 1 % of ice overburden pressure for Budd runs and 0.4 % of ice overburden pressure for Schoof runs, due to numerical artefacting that arises for values smaller than this. The impact of these choice of caps is discussed in Appendix A3. Inverting for the basal friction coefficient using a hydrology model output effective pressure and a friction law which explicitly depends on the effective pressure is similar to the work of Koziol and Arnold (2017).
We present the results of the GlaDS modeling followed by the inversion results using ISSM and compare the impact of using NG and NO for the Schoof and Budd friction laws.
3.1 Subglacial hydrology
The GlaDS modeling indicates that major subglacial hydrology channels form in the Denman–Scott catchment as seen in Fig. 2a, with significant discharge through both the Denman (Fig. 2a (i)) and Scott glaciers (Fig. 2a (ii)). For the former, the discharge is 15.8 m3 s−1, and for the latter, it is 6.0 m3 s−1. The Denman channel initiates as two branches that flow for 80 and 52 km through subglacial valleys before merging at the beginning of the significantly overdeepened trough. From here, the channel flows through the base of the trough, emerging 129 km downstream at the grounding line. The Scott channel, on the other hand, converges as a single channel 95 km from the grounding line. There is substantial water convergence with a maximum depth of 25 m in a basal depression (Fig. 2d (v)) that feeds the channels of both the Denman and Scott glaciers, although with the strongest flux towards the eastern branch of the Denman channel (Fig. 2a (iii)). The bed topography of this basin feature lies at 1900 m below sea level (Fig. 1a), and subglacial water flows upslope by approximately 1200 m to drain downstream. A second “lake-like” feature feeds the western branch of the Denman channel (Fig. 2a (iv)) and reaches a water depth of ∼8 m (Fig. 2d (vi)).
As effective pressure is the ice overburden pressure minus the subglacial water pressure, low effective pressure implies high subglacial water pressure, with negative effective pressure implying that the subglacial water pressure is greater than the ice overburden pressure. A low effective pressure corresponds to a high fraction of flotation, which is defined to be the subglacial water pressure divided by the ice overburden pressure. When the subglacial water pressure is equal to the ice overburden pressure; flotation is reached and the fraction of flotation is 100 %. Effective pressure in the GlaDS outputs is lowest in the basin feature (Fig. 2d (v)) and the lake-like feature (Fig. 2d (vi)), reaching −0.4 MPa in the former and −0.25 MPa in the latter (Fig. 2b). This translates to a maximum of 101 % fraction of flotation relative to the ice overburden pressure for both sites (Fig. 2c). Other regions of negative effective pressure occur in the troughs of Denman and Scott glaciers. Low effective pressures, close to zero (>92 % fraction of flotation), persist through the upper regions of the domain and through the subglacial valleys and troughs of both glaciers. Close to the grounding line, effective pressures rise to a 90 % fraction of flotation in the regions where the channels of both glaciers exit into the ice shelf cavity. These low and negative effective pressures persist at near steady state and do not cause instability in the GlaDS model.
3.2 Ice dynamics and inversion
We first consider the results for the Schoof friction law (Eq. 6). When using NG (Fig. 3a), the Schoof friction coefficient is estimated using inversion is relatively uniform across the basin, with over 50 % of the basin recording basal friction coefficients that range between 2910 and 3580 kg m s (Fig. 4a). Departures from this occur in isolated patches of very low and very high friction – corresponding to the basin-like feature (Fig. 2d (v)) feeding the eastern branch of the Denman channel and the Scott channel – and a region directly upstream of this feature where the ice surface speeds abruptly decrease. The basal friction coefficients are relatively lower in the Denman and Scott troughs, and alternating high and low “stripes” are evident in the region west of Denman Glacier and south of the Shackleton Ice Shelf (Fig. 1b (i)), where ice surface speeds are <50 m a−1 (Fig. 1a). We note that this region was excluded from the GlaDS model (see discussion in Appendix A1).
The Schoof basal friction coefficient estimated using NG is smoother compared with that using NO (Fig. 4a and d; Appendix B). Here, a smoother basal friction coefficient resulted in lower median differences and root mean square error (RMSE) and higher mean differences between the simulated and observed ice surface speeds for the NG simulation over the NO simulation (Table 2; Fig. 4c, f). We find a positive correlation (r2=0.291) between NG and the basal friction coefficient (Fig. 5a) when considering areas where ice surface speeds are ≥10 m a−1, in the region where ice is more dynamic closer to the grounding line. Similarly, there is a positive correlation between NO and the basal friction coefficient (r2=0.266, Fig. 5b).
We next compare the difference between basal friction coefficients with the Schoof and Budd friction laws and using NG. All basal friction coefficients are relatively lower in the Denman and Scott Glacier troughs and the lake-like and basin features that feed into these troughs and relatively higher in the region between the lake-like and basin features (Fig. 4b,h). However, the range of the basal friction coefficient across the catchment estimated with the Budd friction law is substantially greater than that with the Schoof friction law, with relatively high values in the interior of the catchment compared with towards the coast. This leads to a greater standard deviation in the Budd friction coefficient compared with that of Schoof when both are normalized to their respective means (0.228 for Schoof and 0.385 for Budd). Unlike the Schoof friction law, the choice of NO or NG has a significant impact on the distribution of the basal friction coefficient (Fig. B4) in the Budd friction law. This is because the upglacier portions of the catchment fall into a Weertman sliding regime in the Schoof friction law where τb≪CmaxN (Fig. A6) and . Here, the choice of effective pressure will have minimal effect on the Schoof friction coefficient. In the Budd friction law, meaning that the effective pressure continues to play an important role in the basal friction coefficient throughout the entire catchment, which propagates the large discrepancy between NO and NG to the basal friction coefficients obtained from using these various effective pressures. Despite the large range of the Budd friction coefficient across the catchment, the Budd friction coefficient is generally smoother than the Schoof friction coefficient, which may be a consequence of the choice of the Tikhonov regularization coefficient used in the inversion procedure. The Budd friction law also predicts lower mean and median differences between the simulated and observed ice surface speeds than the Schoof friction law (Table 2), for both NO and NG. Finally, in contrast with the Schoof result, there is a negative correlation between NG and the basal friction coefficient using the Budd friction law, with basal friction coefficients increasing for decreasing effective pressures (; Fig. 5c).
4.1 The GlaDS effective pressure
The GlaDS-calculated effective pressure in this region of the Antarctic reflects variability in ice thickness and basal topography, with regions of negative effective pressure concentrated in deep basins and subglacial valleys and troughs. The effective pressure in the remainder of the domain, including the interior of the catchment, is close to zero with minimal spatial variation. A region of high effective pressure is found to the north of the Denman trough. In this area, the bed topography is steep and the ice thin; this causes water to drain rapidly and lower the local water pressure. It is unclear whether this is a realistic representation of the basal hydrology in this region or whether the low pressures indicate a limitation of the hydrology model – in situ data collection is necessary to investigate this further.
4.2 Effective pressure and basal sliding laws
Using the Schoof friction law in the ice dynamics inversions, the resulting basal friction coefficient is smoother when NG is applied, compared to when NO is applied, irrespective of the Tikhonov regularization coefficient (Appendix B). When NO is used, large gradients in the basal friction coefficient appear that are not aligned with or perpendicular to flow, and naturally increase as the Tikhonov parameter decreases. Previous studies have linked flow-perpendicular gradients in the basal friction coefficient to topography or hydro-mechanical feedbacks (Sergienko et al., 2014); however, that is unlikely to be the case here given that our gradients are not aligned with or perpendicular to flow. Rather, it is likely that these large gradients in the basal friction coefficient are artefacts of the inversion procedure that arise only when NO is used.
Our parameterization for NO is a function of the ice overburden pressure and bed elevation (), which has been commonly used to define the effective pressure in ice sheet modeling (Åkesson et al., 2022; Yu et al., 2018; Åkesson et al., 2021). This parameterization is not the actual subglacial hydrology effective pressure, defined as , where pi and pw are the pressures of ice and water, respectively. Indeed, if at the overburden pressure, the effective pressure N should be zero everywhere. However, the definition for NO here is in fact the overburden hydraulic potential, not an effective pressure, and will produce high effective pressures (i.e., low basal water pressures) in regions of thick ice grounded above sea level and low effective pressure in regions of thick ice grounded below sea level (Fig. 3b), yielding high friction in interior regions of the Denman catchment. The water pressure in this parameterization has a dependence on bed elevation, resulting in areas of negative water pressure anywhere the bed is above sea level, which is not physical. Previous studies have investigated alternative parameterizations for the effective pressure due to the computational cost or numerical instabilities associated with coupling of an ice sheet model to a complex subglacial hydrology model. Full coupling between these models is a recent development in the field (Cook et al., 2022). For example, Brondex et al. (2017) use for areas grounded below sea level and NB=ρigH for areas grounded above sea level, which leads to zero water pressure above sea level. Bueler and van Pelt (2015) used a parameterization that incorporates the effects of till at the base. Our results suggest that ice sheet models would be more accurate with improved representations of the effective pressure in friction laws (e.g., Cook et al., 2022), highlighting the need for further investigation into appropriate parameterizations for the effective pressure when a subglacial hydrology model output is not directly available. To this end, we suggest an alternative parameterization for the effective pressure in Sect. 4.3 below.
Our finding of a positive correlation between NG and the basal friction coefficient for the Schoof friction law is unsurprising: effective pressures are lower where there is substantial basal lubrication, which generally leads to a lower basal friction coefficient and faster surface flow. This is evidence for subglacial hydrology being a key control on the surface speeds in regions of lower effective pressures for Denman Glacier. However, in the regions of fastest flow at Denman Glacier, particularly towards the grounding line, the model also suggests the presence of substantial channels, which have previously been associated with increasing effective pressures and slowing ice flow (Nienow et al., 2017). The simultaneous presence of channels, fast ice flow, and low effective pressures has also been observed at Aurora Subglacial Basin (Dow et al., 2020) and in the Weddell Sea region (Dow et al., 2022). This is due to the channels being near steady state, which means they operate at effective pressures closer to that of the surrounding distributed system, compared to alpine-like hydrologic systems (Iken and Bindschadler, 1986). In the inland regions of the Denman catchment, distributed system effective pressure is also near zero, yet ice flow speed in this region is low. Here, longitudinal friction plays a key role in restraining inland ice flow, with fast flow near the grounding line requiring a combination of high driving stress and low basal effective pressures.
The relationship between low effective pressures and low basal friction did not hold for the Budd friction law (Figs. 3a, 4e) despite a stronger negative correlation between the basal friction coefficient and ice surface speed for the Budd friction law compared to the Schoof friction law (Fig. 5a, c). Our study uses inverse methods to calculate the basal friction coefficient for a given τb and NG. For the Budd friction law, this yields a basal friction coefficient that adjusts to compensate for large variations in basal ice velocities – the basal friction coefficient decreases as basal ice velocities increase and vice versa, and there is a strong inverse relationship between these two fields. However, the requirement for such a strong inverse relationship between basal friction coefficient and basal ice velocity leads to a generally weaker relationship between the basal friction coefficient and the effective pressure, since the effective pressure depends on factors other than the basal ice velocity: the topography, the ice thickness, and water accumulation along with efficiency of basal channels to some extent. For the Schoof friction law there is no requirement for a strong correlation between the basal friction coefficient and basal ice velocities, and we find that spatial variations in the basal friction coefficient more closely mirror those of the effective pressure, resulting in the stronger correlation between these two fields.
The choice of friction law plays an important role in grounding line migration and mass loss, with studies showing that friction laws that incorporate effective pressure yield more realistic representation of ice sheet dynamics, grounding line retreat, and mass loss than those that do not (Åkesson et al., 2021; Brondex et al., 2019). In an idealized flow line model, Brondex et al. (2017) showed that when using a Weertman friction law, which does not consider the strong dependence of τb on the effective pressure, the grounding line can be stable on a retrograde slope, contrary to results from theoretical analysis (Schoof, 2007), and recommended the use of the Schoof friction law due to its strong physical basis. Regularized Coulomb friction laws have also demonstrated superior performance over the Budd friction law in capturing observed grounding line migration (Joughin et al., 2019) and more realistic estimates of future ice mass loss (Brondex et al., 2019). These results are largely a consequence of the cavitation effects captured in regularized Coulomb friction laws, which prevent unphysical, unbounded basal shear stresses for increasing ice velocities (Schoof, 2005). By contrast, the basal shear stress is unconstrained for finite values of m in other friction laws, including the Weertman or Budd friction laws. Consistent with Brondex et al. (2017, 2019), our results favor the use of a regularized Coulomb friction law combined with effective pressure output from the GlaDS model for more realistic representation of the basal friction coefficient calculated using inversion.
In this work we have used an SSA ice flow model, which fails to capture bed-parallel vertical shear deformations. This may affect the results of our inversions for the basal friction coefficient in areas of non-negligible vertical shear, such as at the onset of fast-flowing ice streams of the Denman and Scott troughs. However, the use of the Glen flow relation may also impact the capacity of even higher-order models to accurately capture bed-parallel vertical shear deformations. For example, McCormack et al. (2022) showed that even when the BP approximation to the full Stokes equations is used, the unenhanced Glen flow relation fails to capture the vertical shear profile expected in regions where ice is frozen to the bed of the glacier. In our simulations that use the SSA approximation and the Glen flow relation, it is possible that sliding is overestimated, and the basal friction coefficient underestimated, where vertical shear is an important deformation process. However, this is a common issue to all of our model runs and is therefore unlikely to alter the main conclusions of this work that compare how the form of the effective pressure impacts the basal friction coefficient.
Initializing our temperature field with surface temperatures from RACMOv2.3 (van Wessem et al., 2018) could have an impact on the rigidity and basal friction coefficient inversions we performed in this work. Zhao et al. (2018) showed that initializing a model with a colder temperature field resulted in a decrease in the basal friction coefficient computed from inversion, due to stiffer ice. It is possible that if we used a thermal model and computed depth-averaged temperatures to use in our model initialization – effectively increasing the initialization temperatures – we could similarly see an increase in the basal friction coefficients. Like with the SSA approximation, these are issues that would be common in all of our model runs and would be unlikely to alter the main conclusions we came to regarding the effect of the form of the effective pressure on the computation of the basal friction coefficient.
The GlaDS model used here has been demonstrated to accurately represent observed properties of the Antarctic subglacial hydrologic system (Dow et al., 2018, 2020, 2022). Nevertheless, it has not yet been coupled with an ice sheet model for prognostic simulations of Antarctic ice dynamics. Given the key role that effective pressure plays in basal friction laws, it is clear that future implementations of ice dynamics models should explicitly take relevant subglacial hydrology processes into account, such as spatially and temporally variable basal water pressure, if they are to accurately predict climate-driven changes in ice dynamics. A recent study by Kazmierczak et al. (2022) coupled simplified representations of subglacial hydrology to ice dynamics using the Budd friction law. This study by Kazmierczak et al. (2022) examined the impact of different representations of the effective pressure – approximated in turn by height above buoyancy (van der Veen, 1987; Huybrechts, 1990; Winkelmann et al., 2011; Martin et al., 2011), reduced by the subglacial water pressure (modified from Bueler and Brown, 2009) or sliding related to water flux (Goeller et al., 2013) from a simple subglacial water routing model (Le Brocq et al., 2009) and an effective pressure in a till (Bueler and van Pelt, 2015) – on ice mass lost over Antarctica over the 21st century. A key finding is that when the effective pressure is low in the grounding zone, the ice sheet is more sensitive to increases in climate forcing, ultimately increasing mass loss and resulting in sea level rise. The accuracy of these simplified subglacial hydrology models in representing the full physics of the ice–bed interface still needs to be tested. However, we found that the GlaDS model predicted low effective pressures in the grounding zone of the Denman–Scott catchment, which, following Kazmierczak et al. (2022), suggests this region might be sensitive to enhanced mass loss in the future. The next step is to couple the GlaDS model for a more complex representation of the basal hydrology with an ice dynamics model for application to Antarctic glaciers – an approach that has been applied to Greenland glaciers (e.g., Cook et al., 2022).
4.3 Empirical parameterization
We suggest here an alternative parameterization of effective pressure, NE, which is empirical in nature and based on spatial patterns in the GlaDS effective pressure output of Denman Glacier. Given the high correlation of water pressure to ice overburden pressure (r2=0.998 over the Denman domain) and the relatively constant effective pressure in the interior of Denman, we suggest a form of the water pressure proportional to the ice overburden pressure multiplied by a term (), where is a saturation term such that the water pressure reaches a maximum fraction of flotation in areas of high ice thickness. Here, is defined in Eq. (8) and x is defined in Eq. (9). The proposed parameterization is as follows:
Here, ρi is the density of ice, g is the gravitational acceleration, H is ice thickness, and the effective pressure is given by Eq. (7). γ is a constant representing typical effective pressure as a fraction of ice overburden in areas of high ice thickness, rl is the water pressure as a fraction of ice overburden pressure as the ice thickness goes to zero, Ht is a typical large ice thickness, Hs is a typical small ice thickness, and ϵ is a small constant taken so that in areas of low ice thickness (Hs) the water pressure is at a typically low value (rl+ϵ)ρigH. GlaDS output effective pressure and ice thickness data for Denman Glacier were used to find values for rl, γ, Hs, and Ht used in NE (ϵ was taken to equal 0.05). This resulted in rl=0.7, γ=0.96, Hs=500 m, and Ht=2800 m. The empirical nature of this parameterization means that these constants are not “optimized” and the use of NE on future glaciers may require an alteration of these constants, specifically Hs and Ht, which can be determined using ice thickness data. More details on the parameterization and its derivation are provided in Appendix C.
The proposed empirical effective pressure (NE), is compared with the GlaDS output effective pressure (NG), the typically prescribed effective pressure NO, and the Brondex et al. (2017) prescribed effective pressure (NB, Sect. 4.2) – all as a fraction of ice overburden pressure (Fig. 6). NE matches NG more closely than either NO and NB (note the different color bar scalings in Fig. 6a and b compared with Fig. 6c and d). This is further illustrated in Fig. 6e, which shows the difference between NE and NG as a fraction of ice overburden pressure: 72.9 % of the domain has a difference of <0.02 and 95.6 % of the domain has a difference of <0.05. Figure 6f shows the saturation term (gx(H)) as well as the physically equivalent scatter from the GlaDS output data, . The form of the scatter is qualitatively similar to the saturation term in areas of moderate to high ice thickness. The mean percent difference excluding the area of the domain where NG≤0, given by , is 118 %, though 45.4 % of this domain falls below 20 % and 77.0 % of the domain falls below 50 %, so the data are skewed from the areas of low effective pressure. Standing alone, this does not sound impressive, and it is a result of the variability in the scatter seen in Fig. 6f. This is still a substantial improvement over NO, which has a mean percent difference of 3688 %, where 57.7 % of the domain falls between 2000 % and 3000 % difference.
Though this parameterization is not derived from physical principles and it lacks complete hydrological connectivity to the ocean, it produces physically realizable effective pressures , which match closely to GlaDS output effective pressures over Denman Glacier. The success of NE over NO on Denman Glacier combined with its improved physical implications over NO warrants further investigation into its validity and potential use in future modeling studies, particularly on Antarctic glaciers.
The Denman–Scott catchment in East Antarctica is a region that is undergoing active glacial retreat and is vulnerable to ongoing mass loss as the climate warms. We have investigated the coupled interactions between the subglacial hydrological system and the ice sheet through the basal friction coefficient – an important tuning parameter used in basal friction laws – and the dependence of the basal friction coefficient on the form of the effective pressure. We find that when the Schoof friction law is used, there is a smoother basal friction coefficient, and a slightly improved match between the simulated and observed ice surface velocities is found when NG is used compared to when NO is used. Of the two friction laws tested here – the Schoof friction law and the Budd friction law – only the Schoof friction law captures the predicted relationship between the subglacial hydrological and the ice sheet dynamics. That is, using the Schoof friction law, regions of lower effective pressures (both NG and NO) tend to also have lower basal friction coefficients – evidence for the controlling role of the hydrological system.
In their application of the GlaDS model to the Aurora Subglacial Basin, Dow et al. (2020) showed the utility of geophysical observations in validating and constraining subglacial hydrology model parameterizations. By contrast with the Aurora Subglacial Basin, the Denman–Scott catchment is a relatively sparsely observed system, with few geophysical observations available to constrain our model simulations. Similarly to previous studies (Dow et al., 2020; Schroeder et al., 2013; Hager et al., 2022), we recommend the collection of high-resolution radar surveys (with a particular focus on the regions at the boundary of the distributed and channelized subglacial hydrologic systems and over the lake-like features upstream of the Denman and Scott channels) and direct measurement of subglacial hydrological flux across the grounding line to constrain simulated subglacial meltwater volumes. Such observations are vital in understanding the role of the subglacial hydrological system and its coupled feedbacks with the ice sheet.
We proposed a new parameterization for the effective pressure based on empirical data of the Denman–Scott catchment and the GlaDS results that more closely matches simulated effective pressures than suggestions from previous literature and which shows promise in capturing key features of the subglacial system in the absence of a complete subglacial hydrology model output. However, given the empirical nature of this parameterization, our results highlight the importance of simulating coupled interactions between the subglacial hydrology and ice sheet systems to accurately represent ice sheet flow, and future work should focus on the coupling of these two systems in transient simulations of the Antarctic Ice Sheet.
A1 Grid domains and effective pressure extrapolation
The domains of the ISSM and GlaDS models differ, with the GlaDS domain being a subset of the ISSM domain due to the differing subglacial hydrological and ice catchments and limits to the GlaDS domain size requiring restriction to the primary hydrological outlets of Denman and Scott glaciers. This means that there are regions within the ISSM domain for which the GlaDS effective pressure does not exist. We test the impact of two different methods to extrapolate the effective pressure in regions outside the GlaDS domain.
The first method relates the GlaDS effective pressure to the bed topography and ice surface speed. We partition the domain into three regions based on ice surface speed, with the first encompassing areas where v<500 m a−1, the second where m a−1, and the third where v≥1000 m a−1. We then calculate a linear regression between the bed topography and the GlaDS effective pressures in each region and linearly extrapolate the effective pressure onto the ISSM domain based on this relationship. The regression coefficients are shown in Table A1.
The second method uses NO in the areas of the ice sheet model domain where GlaDS was not run. The resulting effective pressure is shown in Fig. A1. The basal friction coefficient found using method 2 is generally more smooth than the basal friction coefficient found using method 1, with the exception of transition areas between the GlaDS and extrapolated effective pressures to the far west and south of the GlaDS domain. The match between the simulated and observed ice surface speeds is similar for the two extrapolations, though method 2 produces more variation in the region upstream of the grounding line than method 1.
In general, there is no reason that the hydrology and ice dynamics domains should match exactly, and where this is the case, we would always expect regions where we will need to extrapolate boundary conditions between the models. The fact that we have had to extrapolate effective pressure is not in itself a problem, and our results show the importance of taking into consideration the method of extrapolation. There may be other methods for extrapolating the boundary conditions that we did not consider, including parameterizing the effective pressure using different physical assumptions (e.g., Bueler and Brown, 2009; Bueler and van Pelt, 2015; Van der Wal et al., 2013; Huybrechts, 1990; Kazmierczak et al., 2022; van der Veen, 1987; Winkelmann et al., 2011), and future studies should establish best practice in this regard.
A2 Mesh sensitivity
Depending on the degree of regularization in the inverse method, small-scale features may be present in the basal friction coefficient. We perform a sensitivity analysis to examine the extent to which such features may arise due to the location of the nodes and vertices in the mesh. We generate five different meshes by adding a random perturbation in the surface velocity field of up to 20 % of the MEaSUREs velocities and then perform non-uniform mesh refinement on the velocity magnitude, as follows: the typical edge length of the elements is 3000 m; for regions where the ice surface speed is ≥500 m a−1, maximum edge lengths are 500 m; and for regions where the ice surface speed is between 100 and 500 m a−1, the maximum edge length is 2000 m. Elements close to the grounding line have a maximum edge length of 1000 m.
Using each of these five meshes, we calculate the Budd and Schoof basal friction coefficients. The probability distributions of the basal friction coefficients normalized to 3500 kg m s (half of the maximum allowed value of C) are shown in Fig. A2. The mean and standard deviation of the differences between both the Schoof and the Budd basal friction coefficients for the perturbed mesh simulations and that of the original mesh are smaller in magnitude than the mean and standard deviation of the difference between the basal friction coefficients using NG and NO. Hence, the details of the mesh have less of an impact on the basal friction coefficients than the form of the effective pressure. The accumulated probability at the lower and upper bounds in Fig. A2 corresponds to the regions of low and/or negative effective pressure near subglacial lakes and troughs and the area of negative effective pressure upstream of the subglacial lake feeding the Denman and Scott troughs (Fig. 2d (v)), respectively.
A3 Effective pressure cap
For regions of very low or negative effective pressure, the basal shear stresses are also very low or almost vanish. In these regions, the inverse method compensates by increasing the basal friction coefficient upstream and around the region of the anomalously low effective pressure, leading to an underestimate of ice surface speeds there compared with the observations. The ice surface speeds are also generally overestimated in the region of vanishing shear stresses.
To account for this, we adjust the effective pressure in regions of low or negative effective pressure, capping the minimum as a percentage of the ice overburden pressure. We test caps of 0.2 %, 0.4 %, 1 %, 2 %, and 4 %. For a cap of 0.4 %, approximately 2 % of the domain has an effective pressure that is linearly proportional to the ice overburden pressure, that is, 98 % of the effective pressure in the ISSM simulation is derived directly from the GlaDS simulated effective pressure. Increasing the cap to 1 % – which is used in the Budd runs – decreases the area over which the GlaDS effective pressure is used to 96 % and increasing the cap to 4 % decreases the area to 48 %. Hence, we use a cap of 0.4 % for our Schoof simulations and 1 % for our Budd simulations.
A4 Nonlinear Budd friction law
Here, we compare the impact of using a nonlinear exponent in the Budd friction law (Eq. 5; as per Brondex et al., 2017, 2019) to our results with the linear exponent m=1 (as per Åkesson et al., 2021, 2022; Yu et al., 2018; Choi et al., 2021; Baldacchino et al., 2022). We follow the same setup outlined in the main text and in Appendix B, capping the effective pressure at 1 % of ice overburden pressure. The L-curve analysis suggests a Tikhonov regularization value of 0.1 is optimal.
The case had much smaller variance in the normalized basal friction coefficient compared to the m=1 case (Fig. A3). That is, using NG, the variance of the normalized basal friction coefficient was 0.231 for and 0.385 for m=1; using NO, the variance was 0.508 for and 0.628 for m=1.
Despite the smaller normalized variance in the basal friction coefficients for the case, there was a stronger negative correlation between the basal friction coefficient and the effective pressure in areas where ice surface speeds are greater than 10 m a−1 (Fig. A4). Using NG, the correlation was −0.466 for and −0.304 for m=1; using NO, the correlation was −0.592 for and −0.0260 for m=1 (Fig. A4). This strong negative correlation for the case indicates that the basal friction coefficient counteracts the effects of the effective pressure, with high values in regions of low effective pressure and low values in regions of high effective pressure. This suggests that the dependency of the basal shear stress on the effective pressure is too strong when . This behaviour is particularly clear in the large area of high basal friction coefficient centered around 2200 km easting and −400 km northing in Fig. A3c, which corresponds to an area of low effective pressure (Fig. 2b). Here, the use of the effective pressure cap limits runaway values of the basal friction coefficient; increasing the cap to reduce this area of high basal friction coefficient would mean using less of the GlaDS data. For the NO run, the large increase in effective pressure upstream resulted in a strong decrease in the basal friction coefficient upstream (Fig. A3d).
A5 Iken's bound sensitivity
Iken's bound, Cmax, mathematically describes the idea that the bed can support a maximum stress (Iken, 1981; Schoof, 2005; Gagliardini et al., 2007). In the Schoof friction law (Eq. 6), τb cannot exceed CmaxN, where Cmax represents the rheological properties of the till (Brondex et al., 2019) and ranges between 0.17 and 0.84 (Cuffey and Paterson, 2010). To determine an appropriate value of Cmax we test the effect of using , and 0.8 on the Schoof basal friction coefficient using NG and NO.
We arrive at the same qualitative conclusions as in the main text for all four values of Cmax. In all cases, the variance of the normalized basal friction coefficient is smaller for NG than for NO, and it is smaller than the variance of the normalized Budd basal friction coefficient (Table A2). Although the results for , and 0.7 are similar, there is a comparatively large decrease in the mean of the basal friction coefficient and increase in the variance of the normalized basal friction coefficient between Cmax=0.7 and Cmax=0.8 for both NG and NO. For lower values of Cmax, a region of higher basal friction coefficient centered around 2200 km easting and −400 km northing (Fig. A5a) develops in the NG simulations. Solving Eq. (6) for C yields
where we see that in the limit where τb→CmaxN, then C→∞. The region of high basal friction coefficient for lower Cmax also has low effective pressure (Fig. 2a), resulting from τb approaching CmaxN for larger values of N and the basal friction coefficient compensating for this. To prevent potentially infinitely increasing values of the basal friction coefficient, it is possible to increase the cap on the effective pressure, but again, this would reduce the area over which the GlaDS data are used and the effect of using modeled hydrology will be less impactful.
The Schoof friction law is a regularized Coulomb friction law which tends towards a Weertman sliding regime when τb≪CmaxN and a Coulomb sliding regime when τb>>CmaxN. Hence, the value of Cmax will have an effect on what physical processes are being represented in the Schoof friction law. Figure A6 shows for various values of Cmax, where values close to 0 correspond to a Weertman sliding regime and values close to 1 correspond to a Coulomb sliding regime. The choice of Cmax appears to have little effect on where each of the Weertman and Coulomb sliding regimes occur, with distinct locations between the two for all values of Cmax. This suggests that the choice of Cmax will not have a significant impact on which physical processes are being represented throughout the domain and does not justify the use of one value of Cmax over another.
The inversion procedure described in Sect. 2.2 includes six inversions for the Budd and Schoof basal friction coefficients and the ice rigidity. Each inversion minimizes a cost function which penalizes absolute ice surface velocity misfit, logarithmic ice surface velocity misfit, and gradients in the basal friction coefficient/ice rigidity. The cost functions minimized in each inversion are described below.
Here, 𝒥a is the linear velocity misfit cost function; 𝒥l is the logarithmic velocity misfit cost function; 𝒥t is the regularization cost function; u is the modeled ice surface velocity; uobs is the observed ice surface velocity; ε is a small number (around machine precision) acting as the minimum observed velocity in Eq. (B2); S is the two-dimensional spatial domain of the model; and k is taken to be α in a Budd inversion, C in a Schoof inversion, and (the rigidity) in a rigidity inversion. The full cost function to be minimized is given in Eq. (B4) where ca, cl, and ct are the cost function coefficients for 𝒥a, 𝒥l, and 𝒥t, respectively. ct is referred to as the Tikhonov regularization coefficient.
The linear and logarithmic cost function coefficients (ca, cl) were chosen so that their contributions to the total cost function would differ by an order 1 factor in each inversion. An L-curve analysis (Fig. B1; Hansen, 2000) was used during the inversions for the basal friction coefficients to determine a suitable Tikhonov regularization coefficient (ct). Ideally ct is chosen such that any smaller value of ct will not have a significant effect on ca𝒥a+cl𝒥l. For the ice rigidity inversions, the Tikhonov regularization parameter was chosen so that its contribution to the total cost function differed by an order 1 factor from that of the linear misfit. The final basal friction coefficients are shown in Fig. B3, their distribution in Fig. B4, and the final rigidities in Fig. B5. The cost function coefficients for the various model runs are displayed in Table B1.
We propose a new effective pressure parameterization (NE) that is empirical in nature and based on the GlaDS effective pressure output in the Denman–Scott catchment. Effective pressure is defined as , where pi=ρigH is the ice overburden pressure and pw is the subglacial water pressure. In this equation, pi is well known; however, pw is generally poorly known in Antarctica due to few direct measurements. Hence, a more robust parameterization of the effective pressure requires a physically realistic prescription for this subglacial water pressure. The proposed parameterization uses the fact that subglacial water pressure has a high correlation with ice overburden pressure (r2=0.998 over the Denman domain) and that subglacial water pressure as a fraction of ice overburden pressure maintains some value for small ice thicknesses and saturates near 1 for large ice thickness. We suggest the following parameterization for subglacial water pressure:
where rl is the maximum allowable ratio of effective pressure to ice overburden pressure and gx(H) is a saturation term. We specify gx(H) as follows:
Here, is a constant thickness, which is chosen so that an area where water is expected to be pressurized with a fraction of overburden γ will reach this level of pressurization by ice thickness Ht. This yields an of the following form:
The only parameter left to be determined is x, and here we choose a value of x so that at a typically small ice thickness, Hs, the subglacial water pressure as a fraction of overburden pressure will be slightly above the minimum subglacial water pressure allowed: , where ϵ is a constant introduced with the purpose of having the saturation term concave down even in areas where there is low ice thickness on the glacier.
We rearrange Eq. (C3) to solve for x as follows:
The empirical nature of this parameterization brings into question its physical validity. Unlike NO and the Brondex et al. (2017) prescribed effective pressure (NB), the effective pressure proposed in Eq. (C5) does not have complete hydrological connectivity to the ocean at the grounding line. However, NO and NB can obtain values of negative or zero water pressure, which is nonphysical. The effective pressure proposed in Eq. (C5) can reach a minimum water pressure of pw=rlρigH meaning that it will always give a physically possible value of effective pressure.
Delving further into the implications of these effective pressure parameterizations, we see that the hydraulic potential is given by (Werder et al., 2013), which means that using NO gives a constant hydraulic potential and a stagnant subglacial hydrologic system. NB implies constant hydraulic potential for bed elevations below sea level and ϕ=ρwgB, where B is the bed elevation, for bed elevations above sea level, which means that water does not flow when the bed elevation is below sea level, and it flows entirely due to gravity and bed slopes when the bed elevation is above sea level. Taking the gradient in the hydraulic potential for NE yields the following equations:
Here, fB(H) is a dimensionless factor which describes the extent to which ice thickness gradients play a role in the hydraulic potential gradient. Cuffey and Paterson (2010) considered the gradient of the hydraulic potential with a spatially constant fraction of flotation and arrive at the form of Eq. (C6) with fB(H) replaced by the fraction of flotation, fw. Here, fB is not the fraction of flotation but is related to it via Eq. (C8) for fw≥rl.
fB is quadratic in fw and fB≥fw for . This means that for NE, gradients in surface elevation play a more important role in water routing than when constant fraction of flotation is considered. The importance of surface elevation gradients reaches a maximum when fw=0.924, corresponding to fB=1.04, ;values of fw above this threshold will actually result in smaller values of fB and a decreased role in the surface elevation gradient, which is at odds with the analysis of Cuffey and Paterson (2010). Modeling (Dow et al., 2014, 2020) outputs, including those from our Denman analysis, suggest that the fraction of flotation is not uniform (Fig. 6). Compared to using NO, this parameterization provides a more intuitive picture of the subglacial hydrological system, where water can flow throughout the entire domain and where flow is dependent on both the basal topography and the ice thickness/surface topography, as is expected.
Except for the lack of complete hydrological connectivity to the ocean, the proposed parameterization in Eq. (C5) agrees more closely with our understanding of subglacial hydrology than both NO and NB. The improved performance over NO and NB for Denman Glacier is encouraging and warrants further exploration, including exploring the validity of this parameterization compared with subglacial hydrology model outputs and its potential use in future modeling studies where full ice-sheet–subglacial hydrology model coupling is unfeasible. Improvements to this parameterization, such as adding in the effects of ice velocity, may improve upon the performance of this parameterization.
We use version 4.21 of the open-source ISSM software, which is freely available for download from https://issm.jpl.nasa.gov/download/ (ISSM, 2023). Basal and surface topography data can be found at https://doi.org/10.5067/E1QL9HFQ7A8M (Morlighem, 2020). Ice velocity data can be found at https://doi.org/10.5067/D7GK8F5J8M8R (Rignot et al., 2017). RACMO temperature data can be found at https://www.projects.science.uu.nl/iceclimate/models/racmo-data.php#totop (Utrecht University, 2021). Basal meltwater input can be found at https://doi.org/10.5281/zenodo.2651652 (Seroussi, 2019). GlaDS hydrology outputs and ISSM friction coefficients for the Denman–Scott catchment are available at https://doi.org/10.5281/zenodo.7709384 (McArthur et al., 2023).
KM designed and completed the ISSM model runs and produced figures; CFD ran the GlaDS model; FSM and CFD provided project direction. All authors wrote the paper.
At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We thank the Digital Research Alliance of Canada for access to supercomputer resources. We thank the editor – Nanna Bjørnholt Karlsson – and two reviewers – Thomas Zwinger and Elise Kazmierczak – for their constructive comments on this paper.
This research has been supported by the Australian Research Council (ARC) Discovery Early Career Research Award (grant no. DE210101433) and the ARC Special Research Initiative Securing Antarctica's Environmental Future (grant no. SR200100005).
This paper was edited by Nanna Bjørnholt Karlsson and reviewed by Elise Kazmierczak and Thomas Zwinger.
Adusumilli, S., Fricker, H. A., Medley, B., Padman, L., and Siegfried, M. R.: Data from: Interannual variations in meltwater input to the Southern Ocean from Antarctic ice shelves, UC San Diego Library Digital Collections [data set], https://doi.org/10.6075/J04Q7SHT, 2020. a
Åkesson, H., Morlighem, M., O'Regan, M., and Jakobsson, M.: Future Projections of Petermann Glacier Under Ocean Warming Depend Strongly on Friction Law, J. Geophys. Res.-Earth Surf., 126, e2020JF005921, https://doi.org/10.1029/2020JF005921, 2021. a, b, c, d, e
Åkesson, H., Morlighem, M., Nilsson, J., Stranne, C., and Jakobsson, M.: Petermann ice shelf may not recover after a future breakup, Nat. Commun., 13, 2519–2519, https://doi.org/10.1038/s41467-022-29529-5, 2022. a, b, c
Baldacchino, F., Morlighem, M., Golledge, N. R., Horgan, H., and Malyarenko, A.: Sensitivity of the Ross Ice Shelf to environmental and glaciological controls, The Cryosphere, 16, 3723–3738, https://doi.org/10.5194/tc-16-3723-2022, 2022. a, b
Beaud, F., Aati, S., Delaney, I., Adhikari, S., and Avouac, J.-P.: Surge dynamics of Shisper Glacier revealed by time-series correlation of optical satellite images and their utility to substantiate a generalized sliding law, The Cryosphere, 16, 3123–3148, https://doi.org/10.5194/tc-16-3123-2022, 2022. a
Blatter, H.: Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients, J. Glaciol., 41, 333–344, https://doi.org/10.3189/S002214300001621X, 1995. a
Brancato, V., Rignot, E., Milillo, P., Morlighem, M., Mouginot, J., An, L., Scheuchl, B., Jeong, S., Rizzoli, P., Bueso Bello, J. L., and Prats-Iraola, P.: Grounding Line Retreat of Denman Glacier, East Antarctica, Measured With COSMO-SkyMed Radar Interferometry Data, Geophys. Res. Lett., 47, e2019GL086291, https://doi.org/10.1029/2019GL086291, 2020. a, b
Brondex, J., Gagliardini, O., Gillet-Chaulet, F., and Durand, G.: Sensitivity of Grounding Line Dynamics to the Choice of the Friction Law, J. Glaciol., 63, 854–866, https://doi.org/10.1017/jog.2017.51, 2017. a, b, c, d, e, f, g
Brondex, J., Gillet-Chaulet, F., and Gagliardini, O.: Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law, The Cryosphere, 13, 177–195, https://doi.org/10.5194/tc-13-177-2019, 2019. a, b, c, d, e
Choi, Y., Morlighem, M., Rignot, E., and Wood, M.: Ice dynamics will remain a primary driver of Greenland ice sheet mass loss over the next century, Commun. Earth Environ., 2, 26, https://doi.org/10.1038/s43247-021-00092-z, 2021. a, b
Cook, S. J., Christoffersen, P., and Todd, J.: A Fully-Coupled 3D Model of a Large Greenlandic Outlet Glacier with Evolving Subglacial Hydrology, Frontal Plume Melting and Calving, J. Glaciol., 68, 486–502, 1727–5652, https://doi.org/10.1017/jog.2021.109, 2022. a, b, c
Dow, C., Werder, M., Babonis, G., Nowicki, S., Walker, R., Csatho, B., Morlighem, M., Dow, C. F., Werder, M. A., Babonis, G., Nowicki, S., Walker, R. T., Csatho, B., and Morlighem, M.: Dynamics of Active Subglacial Lakes in Recovery Ice Stream, J. Geophys. Res., 123, 837–850, https://doi.org/10.1002/2017JF004409, 2018. a, b, c
Dow, C., McCormack, F., Young, D., Greenbaum, J., Roberts, J., and Blankenship, D.: Totten Glacier Subglacial Hydrology Determined from Geophysics and Modeling, Earth Planet. Sc. Lett., 531, 115961, https://doi.org/10.1016/j.epsl.2019.115961, 2020. a, b, c, d, e, f, g, h, i
Dow, C. F., Kavanaugh, J. L., Sanders, J. W., and Cuffey, K. M.: A test of common assumptions used to infer subglacial water flow through overdeepenings, J. Glaciol., 60, 725–734, https://doi.org/10.3189/2014JoG14J027, 2014. a
Dow, C. F., Ross, N., Jeofry, H., Siu, K., and Siegert, M. J.: Antarctic Basal Environment Shaped by High-Pressure Flow through a Subglacial River System, Nat. Geosci., 15, 892–898, https://doi.org/10.1038/s41561-022-01059-1, 2022. a, b, c, d
Gagliardini, O., Cohen, D., Raback, P., and Zwinger, T.: Finite-Element Modeling of Subglacial Cavities and Related Friction Law, J. Geophys. Res.-Earth Surf., 112, 1–11, https://doi.org/10.1029/2006JF000576, 2007. a, b
Goeller, S., Thoma, M., Grosfeld, K., and Miller, H.: A balanced water layer concept for subglacial hydrology in large-scale ice sheet models, The Cryosphere, 7, 1095–1106, https://doi.org/10.5194/tc-7-1095-2013, 2013. a
Greene, C. A.: RAMP Radarsat Antarctic Mapping Project, https://www.mathworks.com/matlabcentral/fileexchange/52031-ramp-radarsat-antarctic-mapping-project (last access: 30 October 2022), 2022. a
Hager, A. O., Hoffman, M. J., Price, S. F., and Schroeder, D. M.: Persistent, extensive channelized drainage modeled beneath Thwaites Glacier, West Antarctica, The Cryosphere, 16, 3575–3599, https://doi.org/10.5194/tc-16-3575-2022, 2022. a
Hansen, P. C.: The L-Curve and Its Use in the Numerical Treatment of Inverse Problems, in: Computational Inverse Problems in Electrocardiology, edited by: Johnston, P., Advances in Computational Bioengineering, 119–142, WIT Press, Southampton, 2000. a, b
Iken, A.: The Effect of the Subglacial Water Pressure on the Sliding Velocity of a Glacier in an Idealized Numerical Model, J. Glaciol., 27, 407–421, https://doi.org/10.3189/S0022143000011448, 1981. a, b
Iken, A. and Bindschadler, R. A.: Combined Measurments of Subglacial Water-Pressure and Surface Velocity of Findelengletsher, Switzerland, J. Glaciol., 32, 101–119, https://doi.org/10.3189/S0022143000006936, 1986. a, b
Indrigo, C., Dow, C. F., Greenbaum, J. S., and Morlighem, M.: Drygalski Ice Tongue Stability Influenced by Rift Formation and Ice Morphology, J. Glaciol., 67, 243–252, 1727–5652, https://doi.org/10.1017/jog.2020.99, 2021. a
Joughin, I., Smith, B. E., and Schoof, C. G.: Regularized Coulomb Friction Laws for Ice Sheet Sliding: Application to Pine Island Glacier, Antarctica, Geophys. Res. Lett., 46, 4764–4771, https://doi.org/10.1029/2019GL082526, 2019. a
Kazmierczak, E., Sun, S., Coulon, V., and Pattyn, F.: Subglacial hydrology modulates basal sliding response of the Antarctic ice sheet to climate forcing, The Cryosphere, 16, 4537–4552, https://doi.org/10.5194/tc-16-4537-2022, 2022. a, b, c, d
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental Scale, High Order, High Spatial Resolution, Ice Sheet Modeling Using the Ice Sheet System Model (ISSM), J. Geophys. Res., 117, 1–20, https://doi.org/10.1029/2011JF002140, 2012. a, b
Martin, M. A., Winkelmann, R., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 2: Dynamic equilibrium simulation of the Antarctic ice sheet, The Cryosphere, 5, 727–740, https://doi.org/10.5194/tc-5-727-2011, 2011. a
McArthur, K., Dow, C., and McCormack, F.: “Basal conditions of Denman Glacier from glacier hydrology and ice dynamics modeling” datasets, Zenodo [data set], https://doi.org/10.5281/zenodo.7709384, 2023. a
McCormack, F. S., Warner, R. C., Seroussi, H., Dow, C. F., Roberts, J. L., and Treverrow, A.: Modeling the Deformation Regime of Thwaites Glacier, West Antarctica, Using a Simple Flow Relation for Ice Anisotropy (ESTAR), J. Geophys. Res.-Earth Surf., 127, e2021JF006332, https://doi.org/10.1029/2021JF006332, 2022. a
Morland, L. W.: Unconfined Ice-Shelf Flow, in: Dynamics of the West Antarctic Ice Sheet, edited by: Van der Veen, C. J. and Oerlemans, J., Springer Netherlands, Dordrecht, 99–116, ISBN 978-94-009-3745-1, https://doi.org/10.1007/978-94-009-3745-1_6, 1987. a
Morlighem, M.: MEaSUREs BedMachine Antarctica, Version 2, Boulder, Colorado USA, NASA National Snow and Ice Data Center Distributed Active Archive Center [data set], https://doi.org/10.5067/E1QL9HFQ7A8M, 2020. a, b, c, d, e
Morlighem, M., Rignot, E., Seroussi, H., Larour, E., Ben Dhia, H., and Aubry, D.: Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica, Geophys. Res. Lett., 37, L14502, https://doi.org/10.1029/2010GL043853, 2010. a
Morlighem, M., Seroussi, H., Larour, E., and Rignot, E.: Inversion of Basal Friction in Antarctica Using Exact and Incomplete Adjoints of a Higher-Order Model, J. Geophys. Res., 118, 1746–1753, https://doi.org/10.1002/jgrf.20125, 2013. a, b
Mouginot, J., Rignot, E., Scheuchl, B., and Millan, R.: Comprehensive Annual Ice Sheet Velocity Mapping Using Landsat-8, Sentinel-1, and RADARSAT-2 Data, Remote Sens.-Basel, 9, 364, https://doi.org/10.3390/rs9040364, 2017. a, b
Nienow, P. W., Sole, A. J., Slater, D. A., and Cowton, T. R.: Recent Advances in Our Understanding of the Role of Meltwater in the Greenland Ice Sheet System, Current Climate Change Reports, 3, 330–344, https://doi.org/10.1007/s40641-017-0083-9, 2017. a, b
Pattyn, F.: A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes, J. Geophys. Res., 108, 2382, https://doi.org/10.1029/2002JB002329, 2003. a
Poinar, K., Dow, C. F., and Andrews, L. C.: Long‐Term Support of an Active Subglacial Hydrologic System in Southeast Greenland by Firn Aquifers, Geophys. Res. Lett., 46, 4772–4781, https://doi.org/10.1029/2019GL082786, 2019. a
Rignot, E., Mouginot, J., and Scheuchl, B.: MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2, Boulder, Colorado USA, NASA National Snow and Ice Data Center Distributed Active Archive Center [data set], https://doi.org/10.5067/D7GK8F5J8M8R (last access: 6 October 2021), 2017. a, b, c
Schroeder, D. M., Blankenship, D. D., and Young, D. A.: Evidence for a Water System Transition beneath Thwaites Glacier, West Antarctica, P. Natl. Acad. Sci. USA, 110, 12225–12228, https://doi.org/10.1073/pnas.1302828110, 2013. a
Sergienko, O. V., Creyts, T. T., and Hindmarsh, R. C. A.: Similarity of Organized Patterns in Driving and Basal Stresses of Antarctic and Greenland Ice Sheets beneath Extensive Areas of Basal Sliding, Geophys. Res. Lett., 41, 3925–3932, https://doi.org/10.1002/2014gl059976, 2014. a
Seroussi, H., Morlighem, M., Rignot, E., Khazendar, A., Larour, E., and Mouginot, J.: Dependence of century-scale projections of the Greenland ice sheet on its thermal regime, J. Glaciol., 59, 1024–1034, https://doi.org/10.3189/2013jog13j054, 2013. a
Seroussi, H., Nowicki, S., Simon, E., Abe-Ouchi, A., Albrecht, T., Brondex, J., Cornford, S., Dumas, C., Gillet-Chaulet, F., Goelzer, H., Golledge, N. R., Gregory, J. M., Greve, R., Hoffman, M. J., Humbert, A., Huybrechts, P., Kleiner, T., Larour, E., Leguy, G., Lipscomb, W. H., Lowry, D., Mengel, M., Morlighem, M., Pattyn, F., Payne, A. J., Pollard, D., Price, S. F., Quiquet, A., Reerink, T. J., Reese, R., Rodehacke, C. B., Schlegel, N.-J., Shepherd, A., Sun, S., Sutter, J., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., and Zhang, T.: initMIP-Antarctica: an ice sheet model initialization experiment of ISMIP6, The Cryosphere, 13, 1441–1471, https://doi.org/10.5194/tc-13-1441-2019, 2019. a, b
Utrecht University: Ice and Climate: Polar climate modelling, https://www.projects.science.uu.nl/iceclimate/models/racmo-data.php#totop, last access: 28 October 2021. a
van der Veen, C. J.: Longitudinal Stresses and Basal Sliding: A Comparative Study, in: Dynamics of the West Antarctic Ice Sheet, edited by: Van der Veen, C. J. and Oerlemans, J., Glaciology and Quaternary Geology, 223–248, Springer Netherlands, Dordrecht, ISBN 978-94-009-3745-1, https://doi.org/10.1007/978-94-009-3745-1_13, 1987. a, b
Van der Wal, W., Barnhoorn, A., Stocchi, P., Gradmann, S., Wu, P., ry, M., and Vermeersen, L.: Glacial isostatic adjustment model with composite 3-D Earth rheology for Fennoscandia, Geophys. J. Int., 194, 61–77, 2013. a
van Wessem, J. M., van de Berg, W. J., Noël, B. P. Y., van Meijgaard, E., Amory, C., Birnbaum, G., Jakobs, C. L., Krüger, K., Lenaerts, J. T. M., Lhermitte, S., Ligtenberg, S. R. M., Medley, B., Reijmer, C. H., van Tricht, K., Trusel, L. D., van Ulft, L. H., Wouters, B., Wuite, J., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Antarctica (1979–2016), The Cryosphere, 12, 1479–1498, https://doi.org/10.5194/tc-12-1479-2018, 2018. a, b
Werder, M. A., Hewitt, I. J., Schoof, C. G., and Flowers, G. E.: Modeling Channelized and Distributed Subglacial Drainage in Two Dimensions, J. Geophys. Res., 118, 1–19, https://doi.org/10.1002/jgrf.20146, 2013. a, b, c, d, e
Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 1: Model description, The Cryosphere, 5, 715–726, https://doi.org/10.5194/tc-5-715-2011, 2011. a, b
Yu, H., Rignot, E., Seroussi, H., and Morlighem, M.: Retreat of Thwaites Glacier, West Antarctica, over the next 100 years using various ice flow models, ice shelf melt scenarios and basal friction laws, The Cryosphere, 12, 3861–3876, https://doi.org/10.5194/tc-12-3861-2018, 2018. a, b, c, d
Zhao, C., Gladstone, R. M., Warner, R. C., King, M. A., Zwinger, T., and Morlighem, M.: Basal friction of Fleming Glacier, Antarctica – Part 1: Sensitivity of inversion to temperature and bedrock uncertainty, The Cryosphere, 12, 2637–2652, https://doi.org/10.5194/tc-12-2637-2018, 2018. a