initMIP-Antarctica: An ice sheet model initialization experiment of ISMIP6

L10 Budd friction law


General comments
In this study, the authors have investigated the coupled interactions between the subglacial hydrological system and the ice sheet through the basal friction coefficient -an important and challenging tuning parameter used in friction laws -and its dependence on the form of the effective pressure. They also proposed a new empirical formulation of the effective pressure.
Moreover, they highlighted the importance of subglacial processes on ice-sheet dynamics and conclude the need for geophysical observations of these processes would improve models. This paper matches the quality criteria required by The Cryosphere. I have only a few comments and recommendations to improve the quality and scientific rigour of the manuscript as well as its understanding for readers that are less familiar with the subject and the tools used.
Thank you for your comments on the paper, which we have addressed below.
1 My first general comment concerns the choice of the Budd sliding law. Knowing that the value of m has a major influence on ice dynamics, why hasn't the Budd sliding law with the exponent m = 3, which is the most commonly used value instead of m=1? Wouldn't it make it easier to compare to the Schoof sliding law?
Note: we have chosen to use the same convention as Brondex et al. (2019) (their equation 6), such that m = 1/n where n is the exponent in Glen's flow law, so that the m value can be more easily compared to that of the Schoof friction law. Hence, here m = 1/3 corresponds to the value of 3 which the reviewer suggests.
Thanks for this comment. We chose the value m = 1 because it has been used in many previous modeling studies (e.g. Åkesson et al., 2021;Åkesson et al., 2022;Yu et al., 2018;Choi et al., 2021;Baldacchino et al., 2022) that use ISSM as well. Using m = 1/3 with the Schoof law would lead to the same issue as the m = 1 case, of not having the same units of the basal friction coefficient. However, as the reviewer points out, the m = 1/3 case will lead to a similar dependence on the basal sliding speed as the Schoof friction law.
To address this, we ran an additional simulation using the Budd friction law with m = 1/3. We found strong negative correlations between the Budd basal friction coefficient and the effective pressure (-0.466 for N G and -0.592 for N O ) compared to the correlations when m = 1 is used (-0.304 for N G and -0.0260 for N O ). When m = 1/3 is used, the Budd basal friction coefficient counteracts the effects of the effective pressure more than in the m = 1 case, with higher values in areas of low effective pressure and lower values in areas of high effective pressure. This is particularly noticeable in the region of low/negative effective pressure centered around 2200 km northing and -400 km easting, where a large area of high basal friction coefficient develops for m = 1/3. The solution to this would be to raise the cap on effective pressure above 1 % of ice overburden pressure, but as we have already shown in Appendix A3 this would mean replacing a significant portion of the GlaDS effective pressure data with an effective pressure that is linearly proportional to the overburden pressure. We have prepared a description of this analysis for a new Appendix in the updated manuscript.
The text from this appendix is as follows: Here, we compare the impact of using a nonlinear exponent m = 1/3 in the Budd friction law (Eq. 6; as per Brondex et al., 2017,0;Kazmierczak et al., 2022) to our results with the linear exponent m = 1 (as perÅkesson et al., 2021;Åkesson et al., 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 m = 1/3 case had much smaller variance in the normalized basal friction coefficient compared to the m = 1 case (Fig. A3). That is, using N G , the variance of the normalized basal friction coefficient was 0.232 for m = 1/3 and 0.385 for m = 1; using N O , the variance was 0.532 for m = 1/3 and 0.628 for m = 1.
Despite the smaller normalized variance in the basal friction coefficients for the m = 1/3 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 N G , the correlation was −0.402 for m = 1/3 and −0.304 for m = 1; using N O , the correlation was −0.528 for m = 1/3 and −0.0260 for m = 1 (Fig. A4). This strong negative correlation for the m = 1/3 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 m = 1/3. 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 N O run, the large increase in effective pressure upstream resulted in a strong decrease in the basal friction coefficient upstream (Fig. A3d).
My second general comment concerns the calculation time required for a more complex hydrological model. At what time and space scale does a more complex model (such as GlaDS) make a significant and crucial difference in glacial dynamics and does it compensate for the additional time of calculations used ?
The space and time scale depends on the types of questions that are being answered. For short time scales (weeks to months), hydrology is likely to be fairly static. For very long time scales (multiple centuries) the hydrology will change but it will take a lot more computational resources to run a fully coupled ice sheet-subglacial hydrology model on these timescales, so if alterations in the the basal boundary conditions are not important to the ice dynamics modeling, then coupled variable hydrology will not be as important.
The area of interest is also an important factor when considering what spatial and temporal 3 scales are important for ice-hydrology interactions. For example, the dynamics of ice streams in the Siple Coast region are strongly dependent on ice-hydrology interactions, so coupling on shorter timescales (e.g. yearly to centennial timescales) would be important here (see e.g. Bougamont et al., 2015). Ice-hydrology coupling may also be essential for modelling the onset and location of ice streams (Kyrke-Smith et al., 2015). A recent study in Greenland showed that the impact of geothermal heat flow on basal ice temperatures is elevated when ice sheet and subglacial hydrology models are coupled (Smith-Johnsen et al., 2020), although this kind of analysis has not been tested more broadly across Antarctica or Greenland. In Antarctica, ice-hydrology coupled models may be essential in predicting the magnitude of basal melt and where it is injected into ice shelf cavities to accurately predict the impact on ice shelf melt; however, such analysis has not yet been conducted.
The end goal with ice-hydrology coupling is to have either an accurate parameterization of hydrology for ice dynamics modeling in situations that require too much computational time for hydrology modeling, or a coupled setup that can run efficiently. The former is what we're aiming for in this manuscript and the latter is something we're working towards. Because the latter hasn't been achieved yet (i.e. coupling of hydrology and dynamics in Antarctic) it's not yet possible to accurately define what the appropriate temporal and spatial scales would be, both for the Denman-Scott catchment, and more broadly across Antarctica.
It's likely that we have not yet answered well enough the question of when we need ice-hydrology coupling and what impact it makes, which is a strong motivator for further studies of its importance.
My third general comment concerns the figures. If a standardisation does not allow the scales to be the same, it must be stipulated in the text for all figures concerned so that there is no misunderstanding. Also, note that this difference in scale does not allow the same analysis and comparison resolution. Also, it is better to have complete captions (variable symbol-units) and the same than the legend (the text written next to the colorbar). It also is more readable if both limits of the scale are written.
To address this, we have altered sections of our analysis to consider normalised deviations from the mean of the basal friction coefficient across the grounded domain so that the same analysis can be performed for the Budd and Schoof friction laws. This is particularly the case when variances are considered, e.g. line 148 of the original manuscript. We have normalized the basal friction coefficients by their respective means in most updated figures. However, we keep some of the original -non-normalized -basal friction coefficients where we need to make comparisons, e.g. in comparing N G with N O in the Budd friction law where there is a substantial change in mean between the two. Limits of the colorbar have been added for all figures, and where they differ between subplots a note has been added to the caption to mention this. Complete captions have been added for all figures where applicable, though we do not add complete captions to the colorbar labels for figures where it will lead to awkward spacing. We agree that it would be best to define N O and N G once and then reference them as such throughout the text. We chose to look at N O as opposed to the version from Brondex et al. (2017) because N O has been used in many modeling studies and we wanted to determine whether this was a reasonable value to use for the effective pressure. We will keep the Brondex et al. (2017)    is the overburden hydraulic potential, but is often called the effective pressure, which is not good practice and something we wish to highlight with this manuscript. N O and N G are not the same quantity and that is why they differ so markedly. To clarify in the updated manuscript, at line 240-242 we modify: However, the definition for N O used here will produce high effective pressures in regions of thick ice grounded above sea level... to However, the definition for N O 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... I therefore propose some small changes in the text or the figure calls to improve clarity and understanding. The main thing is the insistence on the terms << basal >> and << subglacial >> which for me are important to keep throughout the text. Finally, I also propose to elaborate on more technical details with respect to the tools used and the choice of parameters.
We have added the words subglacial and basal throughout the manuscript where applicable.
Specific comments L37-L117. Please to specify whether these negative effective pressures are stable, at the steady state, seasonally dependent. . . Please add information on the stability of this case and also whether it varies over time (depending on the seasons or the tides). It could be interesting to add a sentence mentioning that the presence of zones with negative effective pressure may be persistent and does not lead to instability.
The model reaches near steady state; there are some small temporal changes in water depth in deep pockets on the order of mm per day. There is no seasonality, as there is no surface water input to the subglacial hydrology system and tides weren't considered in the model setup. Over longer timescales there may be changes due to subglacial lake drainage and/or changes in ice geometry, but these were not apparent in our models. To clarify, on lines 90-94 we have changed: We assume a temperate bed throughout, although regions with zero water input (in the interior or southernmost-regions of the domain) are essentially frozen to the bed. The model is run for 10,000 days, providing outputs including channel size and discharge, distributed system discharge, water depth, and effective pressure.
to 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 days to near steady state (there are changes in the water sheet thickness in deep pockets on the order of mm day −1 ), providing outputs including channel size and discharge, distributed system discharge, water depth, and effective pressure.
To address low and negative effective pressures persisting at steady state and not causing model instability we add on lines 182-183: These low and negative effective pressures persist at near steady state and do not cause instability in the GlaDS model. The order of the rigidity inversion and the basal friction coefficient inversions has been added in Section 2.2.1 as well as a reference to Appendix B has been added to line 114. We will add to section 2.2.1 lines 139-147, the following text: 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 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 initital 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.
Additionally, the cost function to be minimized during the rigidity inversions has been added to appendix B with an explanation of the various component cost functions (absolute velocity misfit, logarithmic velocity misfit, and Tikhonov regularization).
L140. As you explain that low effective pressures were associated with faster flow, give a brief explanation of this case.
The sentence in the original manuscript is not quite correct, we have changed it as follows and leave analysis of the role of hydrology in ice dynamics to the discussion section: We find a positive correlation (r 2 = 0.291) between N G and the basal friction coefficient  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). This has been changed, see response to above comment.
L196 : Add a reference on the alpine-like hydrology system.
Added reference Iken and Bindschadler (1986). Equation 1 has been changed to include m and on lines 120-122 we change: This is true for both N G and N O . We have changed on lines 369-370: That is, using the Schoof friction law, regions of lower effective pressures tend to also have lower simulated basal friction and faster flow -evidence for the controlling role of the hydrological system.
to That is, using the Schoof friction law, regions of lower effective pressures (both N G and N O ) tend to also have lower basal friction coefficients -evidence for the controlling role of the hydrological system. This is a good question, and there is no clear justification for one value of C max over another, including the value used in (Brondex et al., 2017). Hence, we test the sensitivity of our inversion results to various values of C max (we will add this sensitivity analysis as a subsection in an appendix of the revised manuscript). We find that a value of C max = 0.7 yields a basal friction coefficient with the smallest variance and we update the manuscript to use this new value. We will include the following text (lines 463-488; new appendix A5): Iken's bound, C max , 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. 7), τ b cannot exceed C max N , where C max represents 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 C max we test the effect of using C max = 0.5, 0.6, 0.7, and 0.8 on the Schoof basal friction coefficient using N G and N O .
We arrive at the same qualitative conclusions as in the main text for all four values of C max .
In all cases, the variance of the normalized basal friction coefficient is smaller for N G than for N O , and it is smaller than the variance of the normalized Budd basal friction coefficient.
Although the results for C max = 0.5, 0.6, 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 C max = 0.7 and C max = 0.8 for both N G and N O . For lower values of C max , a region of higher basal friction coefficient centered around 2200 km easting and -400 km northing (Fig. A5a) develops in the N G simulations. Solving Eq. (7) for C yields where we see that in the limit where τ b → C max N then C → ∞. The region of high basal friction coefficient for lower C max also has low effective pressure (Fig. 2a), resulting from τ b approaching C max N for larger values of N and the basal friction coefficient compensating 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 is 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 << C max N and a Coulomb sliding regime when τ b >> C max N . Hence, the value of C max will have an effect on what physical processes are being represented in the Schoof friction law. Figure F2 (Fig. F2) and Here, the choice of effective pressure will have minimal effect on the Schoof friction coefficient. In the Budd friction law α 2 = τ b /(N u b ), meaning that the effective pressure plays an important role in determining the basal friction coefficient throughout the entire catchment, which propagates the large discrepancy between N O and N G to the basal friction coefficients obtained from using these various effective pressures. The parameter values were obtained through sensitivity testing in the Antarctic tested against geophysical data Dow et al. (2020), and are generally assumed to be applicable for Antarctic glaciers. We have now expanded on the GlaDS methods which justifies the choice of parameters reported in Table 1. We have added the following text (lines 72-76): As discussed in Dow (2023), when the system is overconstricted, the pressures are unrealistically high and the model ceases to converge. When the system is underconstricted, the pressures are below overburden for much of the domain. 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 hydrology pressure in ice sheet dynamics equations. Future work with 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).
In GlaDS application to the Antarctic, the ice flow constant is the same for cavities and channels because they are of a similar size.  Changed 2.2.1 Solving for basal friction coefficients L85. Why eq 2 and eq C2 are not the same? If it's a mistake, the () m is missing in the denominator. Is it possible to isolate the coefficients of friction in order to make them more visible?
Thanks for catching this. Eq 2 is correct, Eq C2 has been changed to match. The basal friction coefficients have been bolded in the equations and we have clarified that the equations are scalar.
L87. Source of the Iken's bound missing and please refer to Appendix C.
Added reference (Appendix A5; Iken, 1981), the appendices have been rearranged and A5 addresses the impact of Iken's bound on the inversions.

L98. For the No equation complete that B takes negative values below sea level and defines the variables.
We have changed the paragraph from lines 148-154: We compare the difference in the friction coefficients when we use two different prescriptions for the effective pressure: (1) an effective pressure given by assuming water pressure equals the ice overburden pressure plus the gravitational potential energy of the water N = ρ i gH + ρ w gB, which we refer to as N O ; and (2) the effective pressure taken directly from the GlaDS simulations, which we refer to as N G .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 D.

to:
We compare the difference in the friction coefficients when we use two different prescriptions for the effective pressure. (1) An effective pressure given by assuming water pressure equals the ice overburden pressure plus the gravitational potential energy of the water N O = ρ i gH + ρ w gB. Here ρ i is the density of ice ρ w , is the density of water, 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 effective pressure taken directly from the GlaDS simulations, which we refer to as N G . 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.
Results 3.1 Subglacial hydrology L105 and L106. Even if it's clearly written in L103, mention that the data indicated (length of the channel and the flow) come from the GlaDs modelling.
We have changed on lines 161-162: Major subglacial hydrology channels form in the Denman-Scott catchment, with significant discharge through both the Denman (Fig. 2a (i)) and Scott Glaciers (Fig. 2a (ii)). to 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)). Changed. Fig.2d(v). I know that the choice of the limits of the legend is there to allow a better reading, but I find it strange to mention 25 m of thickness whereas the legend stops at 10 m.

L111 &
The upper limit of Fig. 2d has been increased to 25 m. L112. We cannot really see on the figure that the strongest flux is toward Denman, so do not mention the figure but rather the data.
We have changed (lines 166-171): There is substantial water amalgamation 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 Denman (Fig. 2a). 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 northern branch of the Denman channel and reaches a water depth of ∼8 m (Fig.   2d (vi)). to 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 125 (vi)).
These are already shown in Fig 2d. They have now been referenced in the text as follows on lines 176-178: 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) 3.2 Ice dynamics and inversion L127. Schoof friction law (éq. 2) Changed L132. slow ice flow (maybe mention the Fig. 1b with the surface speed) This sentence was a bit unclear so we have changed lines 189-192: The friction coefficients are relatively lower in the Denman and Scott troughs, and alternating high and low "stripes" are evident in the region of slow flow (ice surface speeds < 50 m a −1 ) to the west of the Denman Glacier and south of the Shackleton Ice Shelf (we note that this region was excluded from the GlaDS model; see discussion in Appendix A).

to
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 the 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).
L133. Locate the Shackleton Ice Shelf in one of the figures (e.g. Figure 1).
The Shackleton Ice Shelf was added into Fig. 1b, as per the response to the above comment.
L135. a space is missing between Fig. and  The paragraph on lines 135-141 of the original manuscript was unclear and has been changed to: The Schoof basal friction coefficient estimated using N G is smoother compared with that using N O (Fig. 4a and Fig. 4d; 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 N G simulation over the N O simulation (Table 2; Fig. 4c,f ). We find a positive correlation between N G and the basal friction coefficient (correlation coefficient r 2 = 0.291, 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 N O and the basal friction coefficient (r 2 = 0.281, Fig. 5b).
See response to line 138 L133-138-140. The limits given to consider an ice flow faster or slower are not the same, why?
On line 133 of the original manuscript, we used ice surface speeds to describe a specific region of the domain but on lines 138 and 140 we were referring to ice surface speeds throughout the whole domain. The words fast and slow have been taken out of these sections to avoid confusion.

L143. Glaciers
We are referring to the troughs so the s should be at the end of troughs instead of Glaciers.
We now have: Denman and Scott Glacier troughs. The reference has been changed L200. : The relationship between low effective pressures and low basal friction did not hold for the Budd friction law, (Fig 3 (a), Fig 4 (e)) despite a stronger negative correlation between the friction coefficient and surface speed for the Budd friction law compared to the Schoof friction law (Fig 5(a),(c)).    Table A1 : 500≥ v < 1000 m a −1 . The sign is not correct ?
This was in reference to Table B1, however, this table has been removed from the manuscript as the relevant data can be summed seen Fig. B1 and C4 of the original manuscript.

L337. : effective pressure calculation/representation
We have changed effective pressure used to form of the effective pressure.   The limits have been added and the period has been added as well.  This is now in the main text, we have added the reference (Iken, 1981).  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 the cap to 4 % decreases the area to 48 %. from Bechmachine v2 (Morlighem, 2020); (b) Ice urface speed (m a −1 ) from MEaSUREs v2 (Rignot et al., 2011(Rignot et al., , 2017Mouginot et al., 2012Mouginot et al., , 2017, using a logarithmic color scale;  -modify overburden pressure by ice overburden pressure. -On L114, you mention << the nothern branch >> of the Denman glacier, but in the caption you mention in Fig2a (iii) a western branch and in (iv) an eastern branch. I'm confused. Please, is it possible to have a clarification/harmonisation between the figure and the text and also to say which one is longer compared to what is written on L108.
It now reads ice overburden pressure. On line 170, northern branch has been changed to western branch.    As figure f is explained before firgure e in the text, I would switch them.
They have been switched. In the caption, use the same formulation than in the fig. 3