Introduction
Basal sliding is a key component of ice dynamics and plays an important role in the motion of fast flowing regions of ice (Engelhardt and Kamb, Reference Engelhardt and Kamb1998). It is incorporated into ice flow models via the use of a sliding law or a friction law, which relates the basal shear stress to the basal sliding speed and other quantities thought to be relevant to the basal shear stress, such as the effective pressure (difference between ice overburden pressure and subglacial water pressure). Recent developments in subglacial hydrology modeling have opened the door for the improved representation of basal conditions in ice-sheet modeling through the development of two-way coupled models, where the effective pressure in the friction law is taken from the subglacial hydrology model. These models have been applied at the catchment scale in both the Greenland and Antarctic ice sheets (Cook and others, Reference Cook, Christoffersen and Todd2022; Pelle and others, Reference Pelle, Greenbaum, Ehrenfeucht, Dow and McCormack2024; Sommers and others, Reference Sommers2024; Meyer and others, Reference Meyer, Warburton, Sommers and Minchew2026). However, while the representation of a temporally variable subglacial drainage system rests on stronger theoretical ground than models that assume temporally constant effective pressure, it also results in amplification of inaccuracies in the subglacial hydrology model and the basal friction law, the latter of which was designed for temporally constant effective pressure. For this reason, the adoption of two-way coupled models by the glaciological community will be inhibited by flaws in both subglacial hydrology models and friction laws.
One significant limitation with two-way coupled models is that they use two independent descriptions of the subglacial system, one from the friction law and one from the subglacial hydrology model. For these two-way coupled models to become useful, it is necessary to develop a description of subglacial conditions unified between the friction law and the subglacial hydrology model. One approach to achieving this is through altering the subglacial hydrology model, applied by Gilbert and others (Reference Gilbert, Gimbert, Thøgersen, Schuler and Kääb2022) and Togaibekov and others (Reference Togaibekov, Gimbert, Gilbert and Walpersdorf2024). Here, we focus on an alternative approach, unifying the description of basal conditions between the subglacial hydrology model and friction law by directly changing the friction law.
Friction laws
One of the early sliding laws is the Weertman sliding law (Weertman, Reference Weertman1957) given by
\begin{equation}
\tau_b=C_w^2\left\lvert\overrightarrow{u_b}\right\rvert^{s-1}\overrightarrow{u_b},
\end{equation}where 
$\tau_b$ is the basal shear stress, 
$C_w$ is the Weertman friction coefficient, 
$\overrightarrow{u_b}$ is the basal ice velocity and 
$s$ is a friction law exponent. Weertman (Reference Weertman1957) developed this sliding law by considering the flow of ice due to regelation and enhanced creep over cubic bedrock protrusions, and determined that 
$C_w$ was dependent on the ratio of the size of the bedrock protrusions to their spacing.
Shortly after, Lliboutry (Reference Lliboutry1959) proposed a friction law that, along with regelation and viscous creep, considered water filled cavitation in the lee of bedrock obstacles. The inclusion of cavitation resulted in a basal shear stress with an effective pressure dependence. Lliboutry and Weertman would continue to refine their theories of sliding throughout the 1960s (Lliboutry, Reference Lliboutry1968), but from these advances, only the original law proposed by Weertman (Reference Weertman1957) is widely used today (Fowler, Reference Fowler2010). However, later studies would continue to find evidence supporting effective pressure dependence of the basal shear stress. Through laboratory experiments, Budd and others (Reference Budd, Keage and Blundy1979) proposed the Budd sliding law, given by
\begin{equation}
\tau_b=C_b^2N^r\left\lvert\overrightarrow{u_b}\right\rvert^{s-1}\overrightarrow{u_b},
\end{equation}where 
$N$ is the effective pressure, and 
$r$ is a positive friction law exponent. The Budd sliding law is a power law that is also commonly used in ice-sheet modeling studies today.
Iken (Reference Iken1981) found that, for a glacier bed modeled after a tilted staircase, there is a maximum ratio of the basal shear stress to the effective pressure. This ratio is equal to the maximum angle on the stoss side of bedrock protrusions and is known as Iken’s bound (
$C_\text{max}$). Schoof (Reference Schoof2005) extended this result to a general bedrock geometry and modeled sliding over a linked cavity system assuming a linear ice rheology. For low ice speed and high effective pressure, the basal shear stress was proportional to the ice speed, and for low effective pressure and high ice speed the bounded nature of the basal shear stress was observed. Schoof (Reference Schoof2005) developed a friction law with the qualitative aspects of his findings, acting like the Weertman sliding law for high effective pressure and low ice speed, and having 
$\tau_b\approx C_\text{max}N$ for low effective pressure and high ice speed. This type of friction law has a continuous transition between Weertman and Coulomb sliding and is known as a regularized Coulomb friction law.
Finally, based on results using a nonlinear ice rheology to model sliding over sinusoidal and ellipsoidal bedrock protrusions, Gagliardini and others (Reference Gagliardini, Cohen, Råback and Zwinger2007) suggested friction laws should follow the following set of conditions:
(1)
$(\frac{\tau_b}{CN})^n=f(\chi),\quad$ where
$\chi=\frac{u_b}{C^nN^nA_s}$,
$C$ is a parameter which is less than or equal to
$C_\text{max}$,
$A_s$ relates the basal shear stress to the basal sliding speed in the limit there is no cavitation (Weertman sliding), and
$n$ is the exponent in Glen’s flow law (Glen and Perutz, Reference Glen and Perutz1955);(2)
$\left\lvert\frac{\partial f}{\partial\chi}\right\rvert_{\chi=0}=1\quad$ (i.e. when
$\chi=0$,
$f(\chi)$ is tangent to
$\chi$);(3)
$f(\chi=\chi_m)=1$ where
$\frac{\partial f}{\partial\chi}=0$ (i.e.
$f(\chi)=1$ at the maximum, where
$\chi=\chi_m$);(4)
$\tau_b \leq NC_\text{max}$ for all
$\chi \gt \chi_m$.
The first condition specifies how the basal shear stress, the basal sliding speed and the effective pressure can be related to each other. The second condition ensures that in the limit that there is no cavitation at the glacier bed (when 
$\chi=0$, 
$u_b=0$, so there is no cavitation because there is no sliding), the friction law behaves like the Weertman sliding law. The third and fourth conditions enforce Iken’s bound. From these four conditions, they proposed a friction law given by
\begin{equation}
\frac{\tau_b}{N}=C\left(\frac{\chi}{1+\alpha\chi^q}\right)^{1/n},\quad \alpha=\frac{(q-1)^{(q-1)}}{q^q},
\end{equation}where 
$q$ is a parameter of the friction law that controls the rate at which 
$\tau_b/N$ decreases after reaching 
$\chi_m$.
Coupled models
To use a sliding law or a friction law in an ice flow model, values of the friction law exponents, the basal friction coefficient, the basal sliding velocity and in some cases the effective pressure must be determined. The friction law exponents are typically taken to be spatially uniform and temporally constant, and the basal sliding velocity comes from the ice flow model itself. When modeling real domains, the basal friction coefficient is typically solved for using inverse methods, so that the modeled ice velocity matches well with observations of surface velocity at the start of the simulation (MacAyeal, Reference MacAyeal1993; Morlighem and others, Reference Morlighem, Seroussi, Larour and Rignot2013). Observations such as ice surface elevation can be equivalently used in place of or in conjunction with the ice surface velocity (Joughin and others, Reference Joughin, Smith and Schoof2019). The effective pressure is difficult to measure using geophysical observations, and it changes over time (Covington and others, Reference Covington, Gulley, Trunz, Mejia and Gadd2020), particularly in regions where there is surface melt reaching the base of the ice. Therefore, the effective pressure used in ice flow models is either parameterized or taken from the output of a subglacial hydrology model (McArthur and others, Reference McArthur, McCormack and Dow2023), which leads to the possibility of coupled models.
Hard-bed subglacial hydrologic systems have been observed to exhibit a variety of drainage configurations, with the ability for self-organizing large drainage channels to form in some glaciers. The characteristics of the subglacial drainage system vary through space and time and have a history of being represented by connected efficient or inefficient components to approximate these changes in flow regime. Inefficient drainage is conceptualized as a system of cavities linked together through narrow orifices, known as a linked cavity system (Walder, Reference Walder1986; Kamb, Reference Kamb1987), while efficient drainage is represented by a network of linked R-channels (Röthlisberger, Reference Röthlisberger1972), which are incised into basal ice. The most sophisticated hard-bed subglacial hydrology models are solved over two-dimensional (2-D) domains and can model both efficient and inefficient drainage. The Glacier Drainage System (Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013, GlaDS) model is solved using finite element methods; each element surface represents a linked cavity system, which is connected to an R-channel that evolves on each element edge. Where drainage is inefficient, channel cross-sectional area will remain close to zero. The Subglacial Hydrology and Kinetic, Transient Interactions (Sommers and others, Reference Sommers, Rajaram and Morlighem2018, SHAKTI) model is also solved over two dimensions using finite element methods. Unlike GlaDS, SHAKTI does not directly model R-channels, but has a water discharge that is parameterized as a function of the Reynolds number, which leads to efficiently flowing elements with higher average gap height that represent R-channel behavior.
The friction laws and subglacial hydrology models discussed so far have focused on hard-bed systems. Although there has been substantial work regarding soft bed friction (Iverson and others, Reference Iverson, Hooyer and Baker1998; Zoet and Iverson, Reference Zoet and Iverson2020) and subglacial hydrology (Walder and Fowler, Reference Walder and Fowler1994; Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000; Flowers and Clarke, Reference Flowers and Clarke2002a, Reference Flowers and Clarke2002b; Bougamont and others, Reference Bougamont, Price, Christoffersen and Payne2011, Reference Bougamont, Christoffersen, Price, Fricker, Tulaczyk and Carter2015), ice-sheet models running over large domains often use hard bed sliding laws (Seroussi and others, Reference Seroussi2020), and hard-bed subglacial hydrology models have been applied at the ice-sheet scale (Ehrenfeucht and others, Reference Ehrenfeucht, Dow, McArthur, Morlighem and McCormack2025). Furthermore, soft-bed friction exhibits a range of behavior that is mathematically similar to that of hard-bed friction in both slow (Zoet and Iverson, Reference Zoet and Iverson2020) and fast (Iverson and others, Reference Iverson, Hooyer and Baker1998) flowing regions, suggesting the applicability of regularized Coulomb friction laws to soft-bed systems. Despite the differences in actual mechanisms by which water is routed and stored in a soft-bed system, soft-bed subglacial hydrology exhibits the same underlying principals as hard-bed hydrology, where ice geometry drives drainage and water is stored locally. For these reasons, we focus entirely on hard-bed friction and subglacial hydrology modeling and further discuss the applicability of these modeling choices to soft-bed systems in the Supplementary Material Text S1.
Effective pressure outputs from hard-bed subglacial hydrology models have been applied in the friction laws of ice flow models. In some cases, this has taken the form of one-way coupling in both Greenland and Antarctica, where GlaDS and SHAKTI are run to steady state, and then the effective pressure is used in an ice flow model (Ehrenfeucht and others, Reference Ehrenfeucht, Morlighem, Rignot, Dow and Mouginot2023; Zhao and others, Reference Zhao2025). However, two-way coupled models have also been applied to the ice sheets, where effective pressure from the subglacial hydrology model is input to the friction law of the ice flow model and the evolving glacier geometry and velocity from the ice flow model is, in turn, input into the subglacial hydrology model (Cook and others, Reference Cook, Christoffersen and Todd2022; Pelle and others, Reference Pelle, Greenbaum, Ehrenfeucht, Dow and McCormack2024; Sommers and others, Reference Sommers2024; Meyer and others, Reference Meyer, Warburton, Sommers and Minchew2026; McArthur and others, Reference McArthur, Dow and Ehrenfeucht2026a). There have also been several studies applying two-way coupling to synthetic domains (Hewitt, Reference Hewitt2013; Hoffman and Price, Reference Hoffman and Price2014; Gagliardini and Werder, Reference Gagliardini and Werder2018; Lu and Kingslake, Reference Lu and Kingslake2024), not necessarily with GlaDS or SHAKTI. However, most two-way coupled models to date use two independent descriptions of subglacial conditions, one from the friction law and one from the subglacial hydrology model. So, it is possible that for a friction law that describes fast flow, the subglacial hydrology model predicts low cavitation, and for a friction law that describes slow flow, the subglacial hydrology model predicts high cavitation. Instead, Gilbert and others (Reference Gilbert, Gimbert, Thøgersen, Schuler and Kääb2022) developed a ‘hydraulically consistent’ framework for describing subglacial conditions, which formulates the friction law and subglacial hydrology model (they used a version of GlaDS) in terms of a shared ‘cavitation state’. However, this framework has only been applied to alpine glaciers (Gilbert and others, Reference Gilbert, Gimbert, Thøgersen, Schuler and Kääb2022; Togaibekov and others, Reference Togaibekov, Gimbert, Gilbert and Walpersdorf2024), and the reformulated subglacial hydrology model strays from the original from Werder and others (Reference Werder, Hewitt, Schoof and Flowers2013) to better match observations of glacier velocity.
The goal of the work here is to analyze the impact of the friction law on two-way coupled ice flow/subglacial hydrology models, using a two-way coupled modeling setup with GlaDS and the Ice-sheet and Sea-level System Model (Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012, ISSM) applied to the Siple Coast, West Antarctica, as a case study to aid our analysis. We use various friction laws in the coupled model, including the Budd sliding law, a version of the friction law from Eqn (3) that is commonly known as the ‘Schoof friction law’, and two regularized Coulomb friction laws that we propose specifically for two-way coupled modeling, which like Gilbert and others (Reference Gilbert, Gimbert, Thøgersen, Schuler and Kääb2022) and Togaibekov and others (Reference Togaibekov, Gimbert, Gilbert and Walpersdorf2024), have a description of basal conditions unified between the friction law and subglacial hydrology model. We additionally test the impact of these friction laws for one-way vs two-way coupled modeling by running the corresponding ice flow models with each of the mentioned friction laws.
Methods
Study area
The Siple Coast in West Antarctica houses topographically unconstrained ice streams that stagnate and reactivate on timescales of multiple centuries (Catania and others, Reference Catania, Hulbe, Conway, Scambos and Raymond2012). These changes in ice motion are thought to be the result of internal ice-sheet dynamics, including changes in the basal conditions of the ice sheet such as the evolution of the subglacial hydrologic system (Beem and others, Reference Beem, Tulaczyk, King, Bougamont, Fricker and Christoffersen2014). The Siple Coast ice streams are underlain with a water-saturated till layer that is close to ice overburden pressure (Blankenship and others, Reference Blankenship, Bentley, Rooney and Alley1986; Engelhardt and others, Reference Engelhardt, Humphrey, Kamb and Fahnestock1990). A mix of fast and slow flowing regions, combined with the characteristic low effective pressure of the Siple Coast ice streams make the Siple Coast an ideal study site for investigating the impact of the friction law on coupled modeling. For example, low effective pressure ensures effective pressure relevance in regularized Coulomb friction laws.
The model domain consists of the combination of the Ross Ice Shelf and the Ross Sea subglacial hydrologic catchment, computed by assuming that water flows in the direction of the negative gradient of overburden hydraulic potential. This is equivalent to assuming flow with a spatially uniform effective pressure. As we do not have estimates of effective pressure prior to running a subglacial hydrology model, this assumption of spatial uniformity or another effective pressure parameterization is necessary to determine a domain boundary. The initial steady-state subglacial hydrology run described in the model initialization section does not have water amalgamation at the domain boundary while using a zero flux boundary condition, suggesting that the assumption of uniform effective pressure was a good approximation for determining the domain boundary. The computational mesh consists of 429 316 elements with a minimum edge length of 225 m. The mesh has a typical edge length of 10 km but is refined near the troughs draining through the Transantarctic Mountains, in the Ross Ice Shelf, in the Lake Vostok area and in areas where ice speed from MEaSUREs v2 (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017) is 
$ \gt 100$ m a
$^{-1}$.
Ice velocity from MEaSUREs v2 (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017), and bed topography and ice thickness from BedMachine Antarctica v3.4 (Morlighem and others, Reference Morlighem2020) over the model domain are shown in Fig. 1. The fast flowing ice streams of the Siple Coast are clearly visible in Fig. 1a, draining into the Ross Ice Shelf from the West. These ice streams have minimal topographic control, with bed elevation change on the order of a couple of hundred meters from trunk to ridge (Fig. 1b), with local deeper depressions within the beds of the ice streams not resulting in locally faster flow (Fig. 1a,b).
Ross Sea subglacial hydrologic catchment and Ross Ice Shelf ice speed, bed topography and ice thickness. (a) is ice speed from MEaSUREs v2 (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011, Reference Rignot, Mouginot and Scheuchl2017; Mouginot and others, Reference Mouginot, Scheuchl and Rignot2012, Reference Mouginot, Rignot, Scheuchl and Millan2017), with Ice Streams A, B, D and E, and Byrd Glacier labeled. (b) is bed topography from BedMachine Antarctica v3.4 (Morlighem and others, Reference Morlighem2020), and (c) is ice thickness from BedMachine Antarctica v3.4 (Morlighem and others, Reference Morlighem2020). The pink lines represent the grounding line. Northing and Easting in all figures are from WGS 84/Antarctic polar stereographic coordinates referenced to the 71
$^\circ$ S latitude.
Subglacial hydrology model
We use GlaDS, implemented in ISSM, with a small change to the sheet equations that describe the linked cavity system. The hydrology sheet evolves as
\begin{equation}
\frac{\partial h}{\partial t}=\omega_s-\nu_s
\end{equation}
\begin{equation}
\omega_s=\begin{cases}
u_b(h_r-h)/\ell_r & h \lt h_r\\
0 & h\geq h_r
\end{cases},
\end{equation}
\begin{equation}
\nu_s=\tilde{A} h\left\lvert N\right\rvert^{n-1}N,
\end{equation}
\begin{equation}
N=\begin{cases}\phi_0-\phi+\rho_wgh & h \lt h_c\\\phi_0-\phi+\rho_wgh\cdot\text{exp}\left(\gamma(\frac{h-h_c}{h_c})\right) & h\geq h_c\end{cases},
\end{equation} where 
$h$ is the hydrology sheet thickness, 
$\omega_s$ is the hydrology sheet opening rate, 
$\nu_s$ is the hydrology sheet closing rate, 
$u_b$ is the basal ice speed, 
$h_r$ is the typical bedrock bump height, 
$\ell_r$ is the typical cavity spacing, 
$\tilde{A}$ is the ice rheological constant multiplied by an order one factor (
$2/n^n$) related to cavity geometry (semicylinder), 
$N$ is the effective pressure, 
$\phi_0=\rho_{w}gB+p_i$ is the overburden hydraulic potential, 
$\rho_{w}$ is the density of freshwater, 
$g$ is the gravitational acceleration, 
$B$ is the bed elevation, 
$p_i$ is the ice overburden pressure, 
$\phi=\rho_{w}gB+p_w$ is the hydraulic potential, 
$p_w$ is the water pressure, and 
$h_c$ and 
$\gamma$ are constants controlling the growth of subglacial lakes. Water is routed in the direction of the negative hydraulic potential gradient; therefore, changes in water pressure drive the evolution of flow through the distributed system. Subglacial lakes form when 
$h\geq h_r$ and fill due to the closing rate term, which becomes negative when the effective pressure becomes negative. Lake depth reaches steady state when 
$N=0$. For 
$h \lt h_c$, the effective pressure is calculated by taking the subglacial water pressure at the water surface. However, GlaDS sometimes produces areas with large hydraulic potential lows, so to prevent the runaway growth of subglacial lake depth we allow the effective pressure to exponentially increase as a function of lake depth for 
$h\geq h_c$. In reality, a system where 
$\phi \gt \phi_0$ would result in a visco-elastic response of the ice sheet, which would spread water away from the subglacial lake. GlaDS lacks such a mechanism and our exponential increase in effective pressure is meant to spread water away from these hydraulic potential lows in the absence of such a visco-elastic mechanism. A detailed description of the full GlaDS equations, including the channel equations, which we do not alter, can be found in Werder and others (Reference Werder, Hewitt, Schoof and Flowers2013).
GlaDS requires inputs of ice velocity, glacier geometry and basal melt rate. For all runs, melt water is taken from the initMIP Antarctica ISSM Control run final time step (Seroussi and others, Reference Seroussi2019). During the initialization procedure (outlined in more detail below) when GlaDS is run uncoupled to steady state, ice velocity is taken from MEaSUREs v2 (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017) and glacier geometry is taken from BedMachine Antarctica v3.4 (Morlighem and others, Reference Morlighem2020). Following that, ice velocity and glacier geometry are taken from the ice flow model. A zero flux boundary condition is imposed on the domain boundary, and a 
$\phi=0$ Dirichlet boundary condition is imposed at the grounding line.
Table 1 provides a list of parameters used in GlaDS. GlaDS has a history of application in Antarctica (Dow and others, Reference Dow, McCormack, Young, Greenbaum, Roberts and Blankenship2020, Reference Dow, Ross, Jeofry, Siu and Siegert2022; McArthur and others, Reference McArthur, McCormack and Dow2023; Pelle and others, Reference Pelle, Greenbaum, Ehrenfeucht, Dow and McCormack2024; Ehrenfeucht and others, Reference Ehrenfeucht, Dow, McArthur, Morlighem and McCormack2025), and parameter values remain similar across these studies. This informs our choice of model parameters here, except for the sheet conductivity (
$k_s$), which is tuned to produce realistic water pressure (generally close to overburden pressure but without large areas above overburden pressure) for this domain. A value of 
$k_s=10^{-3}$ m
$^{7/4}$ kg
$^{-1/2}$ is used in our study, as lower values of the sheet conductivity lead to large areas of negative effective pressure in the initialization procedure. Effective pressure outputs from alternative sheet conductivity values as well as a more detailed discussion on GlaDS parameter choices can be found in the Supplementary Material Text S2. Further details on GlaDS parameters can be found in Werder and others (Reference Werder, Hewitt, Schoof and Flowers2013).
Table of relevant parameters.
Ice flow model
The ice flow model consists of a stress balance model and a mass transport model implemented in ISSM. For the stress balance model, we use the 2-D shallow shelf approximation (MacAyeal Reference MacAyeal1989, SSA), which includes the basal shear stress as its own term in the balancing equation, as opposed to imposing it as a basal boundary condition as in a three-dimensional (3-D) model. The basal shear stress is calculated from the basal ice speed and the effective pressure using a friction law. We use the depth averaged mass transport model described in Larour and others (Reference Larour, Seroussi, Morlighem and Rignot2012) to update the ice geometry at each time step. The ice flow model requires inputs of effective pressure, basal friction coefficients and ice rigidity. The effective pressure is taken from GlaDS, and the basal friction coefficient and ice rigidity are determined using inverse methods to match simulated and observed ice surface velocity during the model initialization procedure. The ice flow model is run over the Ross Sea subglacial hydrologic catchment and the Ross Ice Shelf. Dirichlet boundary conditions on the ice velocity and thickness are applied at the domain boundary for stress balance and mass transport respectively. The model domain does not evolve in time, which, in conjunction with the boundary conditions of the stress balance and mass transport models, results in a temporally constant calving rate. The ice mask is allowed to evolve if ice reaches zero thickness in the model.
Friction laws
The basal shear stress is calculated using one of four friction laws: (1) the Budd friction law (Budd and others, Reference Budd, Keage and Blundy1979); (2) the Schoof friction law (Schoof, Reference Schoof2005; Gagliardini and others, Reference Gagliardini, Cohen, Råback and Zwinger2007); (3) a version of the Schoof friction law that we have modified, which we refer to as the ‘hydrology-modified Schoof friction law’; and (4) a new friction law we propose that acts as a linear combination of Weertman and Coulomb sliding terms, which we refer to as ‘the linear combination friction law’. The Schoof friction law is given by
\begin{equation}
\tau_b=\frac{C_s^2\left\lvert \overrightarrow{u_b}\right\rvert^{m-1}}{(1+(C_s^2/(C_{\text{max}}N))^{1/m}\left\lvert \overrightarrow{u_b}\right\rvert)^m}\overrightarrow{u_b},
\end{equation}where 
$C_s$ is the Schoof friction coefficient, and 
$m=1/n$ is a friction law exponent. This is a special case of Eqn (3), where 
$q=1$, 
$C=C_\text{max}$ and 
$A_s=C_s^{-2/n}$. 
$\xi_s=(C_s^2/(C_{\text{max}}N))^{1/m}\left\lvert \overrightarrow{u_b}\right\rvert$ can be thought of as a dimensionless parameter in the denominator of the Schoof friction law that controls the ‘sliding regime’ of the friction law. When 
$\xi_s\ll1$ then 
$\tau_b\approx C_s^2\left\lvert \overrightarrow{u_b}\right\rvert^{m-1}\overrightarrow{u_b}$ and the friction law is in a Weertman sliding regime, and when 
$\xi_s\gg1$ then 
$\tau_b\approx C_\text{max}N$ and the friction law is in a Coulomb sliding regime. The sliding regime (
$\xi_s$) therefore depends on the ratio of 
$C_s^2/C_\text{max}$, which varies spatially over the model domain.
The hydrology-modified Schoof friction law is given by
\begin{equation}
\tau_b=\frac{C_h^2\left\lvert \overrightarrow{u_b}\right\rvert^{m-1}}{(1+(\zeta/N)^{1/m}\left\lvert \overrightarrow{u_b}\right\rvert)^m}\overrightarrow{u_b},
\end{equation}where 
$C_h$ is the hydrology-modified Schoof friction coefficient (equivalent to 
$C_s$ in the Schoof friction law and 
$C_w$ in the Weertman friction law), and 
$\zeta$ is the spatially uniform coupling strength parameter that ensures the sliding regime of the friction law is dependent on subglacial conditions related to basal sliding predicted from subglacial hydrology modeling.
We define the degree of cavitation at the glacier bed as
\begin{equation}
\theta=\begin{cases}
h/h_r&h \lt h_r\\
1&h\geq h_r
\end{cases}.
\end{equation} When 
$\theta\approx0$ there is no cavitation and we expect the friction law to be in a Weertman sliding regime, and when 
$\theta\approx1$ there is sliding over water filled cavities and we expect the friction law to be in a Coulomb sliding regime. Therefore, the degree of cavitation at the glacier bed acts as a proxy for sliding conditions predicted from subglacial hydrology modeling, and we would expect the sliding regime to be a monotonically increasing function of the degree of cavitation, going to 0 when 
$\theta=0$, and going to 
$\infty$ when 
$\theta=1$. The relationship between 
$\xi$ and 
$\theta$ should be the same throughout the model domain, so that the sliding regime is consistent with subglacial conditions predicted from the subglacial hydrology model. Throughout the manuscript, we will equivalently say that the sliding regime is ‘spatially consistent’, with the degree of cavitation. Assuming that the hydrology sheet is in steady state, 
$h \lt h_r$ and 
$N \gt 0$ we can set the left hand side of Eqn (4) equal to 0 and rearrange to find
\begin{equation}
\frac{u_b}{N^n}=\frac{\tilde{A}\theta\ell_r}{1-\theta},
\end{equation}where we see that since 
$\tilde{A}$ and 
$\ell_r$ are spatially uniform in this study, if we take 
$\zeta$ to be spatially uniform then the sliding regime of the hydrology-modified Schoof friction law (
$\xi_h=(\zeta/N)^{1/m}\left\lvert \overrightarrow{u_b}\right\rvert$) is spatially consistent with the degree of cavitation at the glacier bed, 
$\xi_h$ goes to 0 as 
$\theta$ goes to 0, and 
$\xi_h$ goes to 
$\infty$ as 
$\theta$ goes to 1. Equation (11) assumes that sliding is the only mechanism by which cavities can open. In reality, cavities also open by melt from turbulent dissipation of potential energy, and frictional and geothermal heat. However, sliding has been shown to be the dominant opening mechanism in a linked cavity system (Walder, Reference Walder1986; Kamb, Reference Kamb1987), and Eqn (11) is therefore a reasonable approximation. Furthermore, any region with sufficiently large melt will result in the formation of an R-channel (whose dominant opening mechanism is melt from dissipation of potential energy) and no longer be part of the linked cavity system.
We note that the hydrology-modified Schoof friction law is equivalent to the Schoof friction law when the ratio 
$C_s^2/C_\text{max}$ is spatially uniform. One way to accomplish this is by taking 
$C_s$ and 
$C_\text{max}$ individually to be spatially uniform. Such an approach is common when applying regularized Coulomb friction laws to synthetic domains with two-way coupling (Hewitt, Reference Hewitt2013; Hoffman and Price, Reference Hoffman and Price2014; Gagliardini and Werder, Reference Gagliardini and Werder2018; Lu and Kingslake, Reference Lu and Kingslake2024). However, the application of regularized Coulomb friction laws to non-synthetic domains typically requires spatially heterogeneous friction parameters to match modeled and observed ice velocity. The Coulomb friction parameter (Iken’s bound) is typically held spatially uniform, while the Weertman friction parameter (the friction coefficient) is varied. This is common in two-way coupled models (Cook and others, Reference Cook, Christoffersen and Todd2022; Pelle and others, Reference Pelle, Greenbaum, Ehrenfeucht, Dow and McCormack2024), as well as one-way coupled or uncoupled models (Brondex and others, Reference Brondex, Gillet-Chaulet and Gagliardini2019; Ehrenfeucht and others, Reference Ehrenfeucht, Morlighem, Rignot, Dow and Mouginot2023). However, this means that there is not a data-informed constraint on the Coulomb friction parameter, which should be the relevant parameter in fast flowing regions. By requiring 
$\zeta$ to be spatially uniform in the hydrology-modified Schoof friction law, varying 
$C_h$ to match modeled and observed ice velocities will simultaneously vary 
$C_\text{max}=C_h^2/\zeta$. Consequently, an appropriate value of Iken’s bound can be determined in fast flowing regions of the domain.
The functional form of the hydrology-modified Schoof friction law and the Schoof friction law does not allow for the ability to both choose 
$C_i$ (
$C_s$ for the Schoof friction law and 
$C_h$ for the hydrology-modified Schoof friction law) and 
$C_\text{max}$ and have a sliding regime dependent on the degree of cavitation at the glacier bed. The Schoof friction law allows for the choice of both 
$C_s$ and 
$C_\text{max}$ but has a sliding regime that is spatially inconsistent with the degree of cavitation at the glacier bed, whereas the hydrology-modified Schoof friction law has a sliding regime that is spatially consistent with the degree of cavitation at the glacier bed but does not allow for the choice of both 
$C_h$ and 
$C_\text{max}$. This is in fact true for all friction laws of the form proposed by Gagliardini and others (Reference Gagliardini, Cohen, Råback and Zwinger2007), through condition 1 outlined in the ‘Introduction’ section. Gagliardini and others (Reference Gagliardini, Cohen, Råback and Zwinger2007) developed their friction law criteria based on simulations of sliding over periodic sinusoidal and ellipsoidal bedrock obstacles. The periodicity of the bed ensured that 
$C$ and 
$A_s$ were spatially uniform in their analysis. However, for a nonuniform bed with varying 
$C$ and 
$A_s$, it is apparent that a functional form that differs to that proposed in condition 1 is necessary to independently choose Weertman and Coulomb friction parameters to match simulated and observed ice velocity, while having a sliding regime that is spatially consistent with the degree of cavitation at the glacier bed. As a simple example of a functional form that accomplishes this, we propose the linear combination friction law, which acts as a linear combination of Weertman and Coulomb sliding terms, with coefficients dependent on the explicit sliding regime, which sum to one and determine the relevance of each term. The linear combination friction law is given by
\begin{equation}
\tau_b=C_\text{max}N\left\lvert \overrightarrow{u_b}\right\rvert^{-1}\overrightarrow{u_b}\left(\frac{\xi_\ell}{1+\xi_\ell}\right)+C_\ell^2\left\lvert \overrightarrow{u_b}\right\rvert^{m-1}\overrightarrow{u_b}\left(\frac{1}{1+\xi_\ell}\right),
\end{equation}
\begin{equation}
\xi_\ell=\alpha\frac{u_b}{N^n},
\end{equation} where 
$C_\ell$ is the friction coefficient of the linear combination friction law, 
$\xi_\ell$ is the sliding regime of the linear combination friction law, and 
$\alpha$ is the coupling strength parameter of the linear combination friction law. 
$\alpha$ serves the equivalent purpose of 
$\zeta^{n}$ in the hydrology-modified Schoof friction law, in the sense that it relates 
$u_b/N^n$ to the sliding regime and controls the sensitivity of the basal shear stress to subglacial hydrologic effects. However, since the linear combination friction law allows for the simultaneous choice of 
$C_\text{max}$ and 
$C_\ell$, the physical meaning of 
$\alpha$ (with respect to bed roughness and angles at the stoss faces of bedrock obstacles) is lost. We choose a value of 
$\alpha$ similar to 
$\zeta^{n}$ in the hydrology-modified Schoof friction law. The linear combination friction law is designed to respond to changes in basal sliding speed and effective pressure in the same way that the Schoof and hydrology-modified Schoof friction law do, in the limit that 
$\xi_i$ goes to 0 or 
$\infty$, where 
$\xi_i$ is 
$\xi_s$ for the Schoof friction law and 
$\xi_h$ for the hydrology-modified Schoof friction law. We acknowledge that the functional form of the linear combination friction law is not derived from a theoretical analysis of basal sliding and first principals. Instead, it is a practical solution to model the evolution of friction on real domains in a way that is consistent with our subglacial hydrology modeling. Future work on specifying the functional form of the friction law while keeping the practical use of said friction law in mind is needed, but beyond the scope of this study.
Finally, for the Budd friction law, we take 
$r=1$ and 
$s=1/3$ so that we have the same effective pressure dependence as the regularized Coulomb friction laws in a Coulomb sliding regime and the same ice speed dependence as the regularized Coulomb friction laws in a Weertman sliding regime.
Model initialization
We aim to approximate the conditions of the ice sheet in 2010 (the start time of the transient simulations) by using an initialization procedure involving inversions for the basal friction coefficient and ice rigidity, model relaxations and running a thermal model to estimate the ice rigidity over the grounded portion of the domain using an empirical parameterization. The initialization procedure can be broken into four main steps described below (a flowchart to help understand the initialization procedure is given in the Supplementary Material Text S6).
(1) Estimate the effective pressure: We run GlaDS for 100 a to steady state using observed ice velocity from MEaSUREs v2 (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017) and ice geometry from BedMachine Antarctica v3.4 (Morlighem and others, Reference Morlighem2020). We start with adaptive timestepping as short as 0.1 days and gradually increase the time step to 0.5 days as the model gets closer to steady state. We test sheet conductivities of
$k_s=10^{-4}$ m
$^{7/4}$ kg
$^{-1/2}$,
$k_s=10^{-3}$ m
$^{7/4}$ kg
$^{-1/2}$ and
$k_s=5\cdot10^{-3}$ m
$^{7/4}$ kg
$^{-1/2}$ but choose to proceed with
$k_s=10^{-3}$ m
$^{7/4}$ kg
$^{-1/2}$, which produced the most realistic effective pressure. The results of these steady-state hydrology runs are discussed in the Supplementary Material Text S3.(2) Invert for the friction coefficient and ice rigidity: Using the effective pressure from steady-state GlaDS, we perform a series of inversions to determine the hydrology-modified Schoof friction coefficient and the ice rigidity. First, we invert for the ice rigidity over the floating portion of the domain and determine it over the grounded portion of the domain by using an empirical parameterization built into ISSM that relates ice temperature to ice rigidity based on Table 3.4 in Cuffey and Paterson (Reference Cuffey and Paterson2010). The ice temperature is obtained by solving a 3-D thermal model (Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012) with 59 6116 elements and five vertical layers, and depth averaging the results. Depth averaging the temperature and computing the ice rigidity will likely lead to a more rigid ice sheet, but as all our models use the same rigidity this is unlikely to affect the comparison of the various friction laws. Parameters for this thermal model are given in Larour and others (Reference Larour, Seroussi, Morlighem and Rignot2012). A Dirichlet boundary condition of modeled atmosphere temperature from RACMO v2.3 (van Wessem and others, Reference van Wessem, van de Berg and van den Broeke2023) is applied at the glacier surface, and geothermal heat flux from the Aq1 model (Stål and others, Reference Stål, Reading, Halpin and Whittaker2020; Reference Stål, Reading, Halpin and Whittaker2021) is applied as a Neumann boundary condition at the glacier base. Next, the basal shear stress is estimated using the Budd friction law and by inverting for the Budd friction coefficient. This initial inversion for the basal shear stress using the Budd friction law is necessary since the hydrology-modified Schoof friction law requires a good initial estimate of the basal shear stress but the Budd sliding law does not. The ice rigidity over the floating portion of the domain is then inverted for once again. Finally, we convert to the hydrology-modified Schoof friction law by solving for the hydrology-modified Schoof friction coefficient given the basal shear stress from the prior friction inversion, and then invert for the hydrology-modified Schoof friction coefficient.
(3) Relax the models individually: Due to feedbacks between GlaDS and the ice flow model leading to instabilities when relaxed together (e.g. negative effective pressure, hydrology sheet growth and ice speedup), the models must first be relaxed separately. We start by running a 20 a ice flow relaxation using effective pressure from the steady-state GlaDS run, the ice rigidity and friction coefficient from step 2, initial ice velocity from MEaSUREs v2 (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017), initial ice geometry from BedMachine Antarctica v3.4 (Morlighem and others, Reference Morlighem2020) and 1 day timestepping. We then use the ice geometry from this ice flow relaxation to invert again for the hydrology-modified Schoof friction coefficient and the ice rigidity over the floating portion of the domain. The idea is that, by the end of the ice flow relaxation, the ice velocity and geometry should be close to observed values. The ice flow relaxation is then rerun with the friction coefficient and ice rigidity from these inversions. The final ice geometry and velocity from this second relaxation are then used to initialize a 70 a steady-state GlaDS run, with a switch from 0.2 day timestepping to 0.5 day timestepping at 2 a, and a switch from
$h_c=100$ m to
$h_c=10$ m at 50 a (due to unrealistically large sheet thickness in the coupled models when
$h_c$ is not reduced).(4) Relax the coupled model: Once the models are relaxed individually they can be coupled and run for a final short relaxation. We take the effective pressure from GlaDS and the basal shear stress from the ice flow model in step 3 and we use them to calculate the friction coefficient of the friction law to be used in the transient simulation. We perform final inversions for the friction coefficient and the ice rigidity to ensure a close initial match to observed ice motion and run a 1 year relaxation of the model to be used in the transient simulation. This last relaxation and the transient simulations have 0.5 day timestepping for GlaDS and 14 day timestepping for the ice flow model.
The ice flow model is forced with mean surface mass balance between 1980 and 2010 (Kittel and others, Reference Kittel2021; Kittel, Reference Kittel2021), and basal mass balance from the initMIP-Antarctica control run final time step (Seroussi and others, Reference Seroussi2019) over the grounded domain and from Adusumilli and others (Reference Adusumilli, Fricker, Medley, Padman and Siegfried2020) over the floating domain.
Two-way coupled model initialization is more complicated than the initialization of one-way coupled or uncoupled models due to the strong dependence of ice flow on the basal shear stress and the strong dependence of the hydraulic potential field on ice geometry. If the ice geometry is not sufficiently relaxed at the start of the transient simulation, this could lead to a subglacial hydrologic response and feedbacks that are forced by the initialization procedure. By first relaxing the models separately, we ensure a relaxed ice-sheet geometry and a relaxed hydraulic potential field for that geometry.
List of experiments
To test the impact of the different friction laws, as well as the impact of two-way coupled vs one-way coupled modeling, we run nine transient simulations (Table 2) from 2010 to 2100 under an SSP5-8.5 surface mass-balance forcing scenario (Kittel and others, Reference Kittel2021; Kittel, Reference Kittel2021). SSP5-8.5 is a high emissions climate scenario resulting in significant global warming and associated climatic impacts. For each of the four friction laws, we run one ice flow-only model and one two-way coupled model. For model-IH and model-CH, we use 
$\zeta=10^7$ m
$^{-4}$ kg s, but as part of a sensitivity analysis, we run a second coupled model with the hydrology-modified Schoof friction law (model-CH2), with 
$\zeta=2\cdot10^7$ m
$^{-4}$ kg s. The results of model-CH2 are discussed in the Supplementary Material Text S4. A list of the experiments with their respective names is given in Table 2.
List of experiments.
Note: H-Schoof refers to the hydrology-modified Schoof friction law.
Results
Ice flow
All ice flow models had similar trends throughout the 90 a simulation (Fig. 2). The Siple Coast ice streams slowed down, except in model-IS, where the northern trunk of Ice Stream B (Whillans) accelerated during the simulation (Fig. 2c). The ice upglacier of Ice Streams A (Mercer), B, D (Bindschadler) and E (MacAyeal) thinned, with more significant thinning of 20–30 m upglacier of Ice Streams A and B, and the region upglacier of the Kamb Ice Stream thickened by 30 m in some places by the end of the simulation. Byrd Glacier accelerated, and the ice downglacier of Byrd, in the Ross Ice Shelf, sped up by 
$ \gt $100 m a
$^{-1}$ as a result of poorly constrained ice rigidity and friction. These changes in ice motion in the Siple Coast ice streams and Byrd Glacier propagated into the Ross Ice Shelf. However, the change in ice thickness of the Ross Ice Shelf did not always follow the trends in ice motion, indicating that the Ross Ice Shelf may be more sensitive to surface and basal mass balance forcing than to dynamic processes.
Changes in ice speed and ice thickness for the ice flow models between 2010 and 2100. (a), (c), (e) and (g) show changes in ice speed (m a
$^{-1}$) and (b), (d), (f) and (h) are changes in ice thickness (m). The model for each row is named on the left. The pink lines represent the grounding line at the start of the simulation and the black lines represent the domain boundary.
The Budd friction law
Model-CB showcased large speedup events of Ice Stream B and the area upglacier of Ice Stream B, resulting in unrealistically large changes in ice speed and thickness on timescales of decades (Fig. 3). Thinning of ice upglacier of Ice Stream B lowered the hydraulic potential, which drew more water into the region, and consequently, lowered the effective pressure. The reduced effective pressure triggered acceleration of inland ice between 2056.8 and 2059.9 (Fig. 3a,b). By 2059.9, this inland ice was moving 
$ \gt $1000 m a
$^{-1}$ faster than in 2058 (maximum speedup from 55.0 to 1494 m a
$^{-1}$). This ice speedup propagated through Ice Stream B, which was followed by a significant slowdown to 77.0 m a
$^{-1}$ by 2080, far below its starting value of 459 m a
$^{-1}$ (note that ice stream speeds, effective pressures and sliding regimes are calculated as an areal average with locations defining the ice streams given in Fig. S3). During this time, when Ice Stream B was decelerating, the subglacial hydrologic system upglacier of Ice Stream B was flooding (hydrology sheet thickness above the typical bedrock bump height) as a result of an extensive hydraulic potential low. The latter was caused by extreme dynamic thinning that occurred following the initial ice speedup event (Fig. 3k). The initial ice speedup event was confined to a faster flowing stream of inland ice, but the expansive hydraulic potential low that it triggered superseded all previous spatial patterns in ice motion (Fig. 3g). This resulted in a large upglacier area of fast moving ice, which funneled into Ice Stream B, increasing in speed to 3318 m a
$^{-1}$ (
$ \gt $5000 m a
$^{-1}$ in some places; Fig. 3d). By the end of the simulation, Ice Streams A and B experienced significant dynamic thickening of 
$ \gt $200 m in some places, and the region upglacier of Ice Stream B experienced significant dynamic thinning of 
$ \gt $100 m in some areas (Fig. 3l).
Outputs for model-CB at various points in time. (a)–(d) plot the change in ice speed from the start of the simulation (m a
$^{-1}$); (e)–(h) plot the hydrology sheet thickness (m); and (i)–(l) plot the change in ice thickness relative to the start of the simulation. Pink lines represent the grounding line at the start of the simulation and black lines represent the domain boundary. The simulation time is given above each column. The change in ice speed has a logarithmic scale in the positive and negative direction, with changes with absolute value 
$ \lt $10 m a
$^{-1}$ in white.
Trends between model-IB and model-CB diverged significantly in the portion of the Siple Coast consisting of Ice Streams A and B. However, these effects were somewhat localized, and trends in Ice Streams D and E were, by comparison, similar between model-IB and model-CB. In model-IB, both these ice streams gradually slowed over the course of the simulation, so that by the end of the simulation ice speed for both ice streams remained within 10 m a
$^{-1}$ of their starting values. In model-CB, Ice Stream E sped up by 24.8 m a
$^{-1}$ over the course of the simulation, with almost 20 m a
$^{-1}$ of this speedup coming within the first 15 a of the simulation. The ice speed variability of Ice Stream E was enhanced in model-CB compared to model-IB, due to changes in the subglacial hydrologic system. In model-CB, Ice Stream D gradually accelerated until mid 2061, where it rapidly slowed by 46.7 m a
$^{-1}$ over the next 3 a as a result of the changes to the ice sheet that took place further South in the Ice Stream B area. An ice speedup event of 78.2 m a
$^{-1}$ took place over 3 a in the early 2080s, which was also a consequence of changes in the Ice Stream B region.
The Schoof friction law
Much like model-CB, model-CS experienced large speedup events in Ice Stream B. In model-CS, these speedup events reached larger extremes than in model-CB but were confined within a 30 a period between 2040 and 2070 (Fig. 4). Like in model-CB, inland ice upglacier of Ice Stream B thinned at the start of the simulation, resulting in a reduction of effective pressure, particularly in regions of faster flowing ice. The ice stream onset of the southern tributary feeding Ice Stream B accelerated from 
$ \lt $232.8 m a
$^{-1}$ in 2047.6 to 8625 m a
$^{-1}$ (maximum speedup) over a span of 
$ \lt $3 a (Fig. 4b,c), resulting in Ice Stream B reaching a speed of 7102 m a
$^{-1}$ (
$10^4$ m a
$^{-1}$ in some regions), well above the starting value of 426.6 m a
$^{-1}$ (Fig. 4d). The speedup events of model-CS were more extreme than those of model-CB. They were initiated earlier in the simulation than in model-CB, and, unlike in model-CB, they terminated by the end of the simulation. The area upglacier of Ice Streams A and B thinned by 
$ \gt $100 m in some places, but the ice thickness of the ice streams themselves thinned by 
$ \lt $50 m, and downstream of Ice Streams A and B, the Whillans Ice Plain thickened by 
$ \lt $50 m (Fig. 4l).
Outputs for model-CS at various points in time. (a)–(d) plot the change in ice speed from the start of the simulation (m a
$^{-1}$); (e)–(h) plot the hydrology sheet thickness (m); (i)–(l) plot the change in ice thickness relative to the start of the simulation. Pink lines represent the grounding line at the start of the simulation and black lines represent the domain boundary. The simulation time is given above each column. The change in ice speed has a logarithmic scale in the positive and negative direction, with changes with absolute value 
$ \lt $10 m a
$^{-1}$ in white.
The hydrology-modified Schoof friction law
As with the coupled models already discussed—model-CB and model-CS—the ice upglacier of Ice Stream B thinned at the start of model-CH. This thinning resulted in a hydraulic potential low and flooding of the subglacial hydrologic system. However, there was no inland ice speedup, despite lower effective pressure. Ice Stream B sped up twice throughout model-CH. In 2063, it increased speed by 52.2 m a
$^{-1}$ in just over 4.5 a, and in 2080 it increased speed by 66 m a
$^{-1}$ in a similar time frame (Fig. 5d). These speedup events were followed by the slowdown of the ice stream, so that by the end of the simulation, Ice Stream B was 28.8 m a
$^{-1}$ slower than its starting speed. The lack of ice acceleration upglacier of Ice Stream B means that dynamic thinning did not occur, and the flooding of the subglacial hydrologic system upglacier in Ice Stream B was not as expansive as it was in model-CB and model-CS. As in model-CB and model-CS, Ice Stream D was affected by the motion of Ice Stream B. However, since the motion of Ice Stream B is much less affected in model-CH than in model-CB and model-CS, the response of Ice Stream D was reduced in model-CH. In fact, only one of the speedup events of Ice Stream B registered in the ice speed trend of Ice Stream D, and it did so as a 4.54 m a
$^{-1}$ increase in ice speed. Aside from this small pulse in ice motion experienced in Ice Stream D in model-CH, trends in ice motion for Ice Streams D and E for model-IH and model-CH stayed within 2 m a
$^{-1}$ of each other.
Time series of ice speed, sliding regime (variable controlling the behavior of the friction law) and effective pressure in Ice Stream B for the models with regularized Coulomb friction laws. The left column shows ice speed (m a
$^{-1}$) as a function of time for the coupled model; the middle column shows sliding regime as a function of time in red and the ratio of effective pressure to ice pressure as a function of time in blue for the coupled model; and the right column shows the sliding regime as a function of time in red and ice speed (m a
$^{-1}$) as a function of time in black for the ice flow model. These outputs are shown for the Schoof friction law (a)–(c), the hydrology-modified Schoof friction law (d)–(f) and the linear combination friction law (g)–(i). Please note the different scales used for plotting outputs from the different models.
As expected, the initial sliding regime of the hydrology-modified Schoof friction law in model-CH aligns with patterns in the initial hydrology sheet thickness. However, the sliding regime of the Schoof friction law in model-CS does not follow these patterns, and areas where hydrology sheet thickness is close to the typical bedrock bump height (implying low effective pressure) have Weertman sliding regimes (
$\xi_s \lt 10^{-10}$). Areas that have low hydrology sheet thickness match more closely between model-CS and model-CH, with model-CS typically classifying these regions as more Coulomb than model-CH. The higher sliding regime of model-CS in these slow moving regions allowed for large changes in the basal shear stress as the subglacial hydrologic system flooded, causing the inland speedup events in model-CS that did not occur in model-CH.
The trend in the logarithm of the sliding regime of Ice Stream B follows very closely to the negative of the trend of the ratio of effective pressure to ice pressure in model-CH (Fig. 5e). Specifically, the correlation between the ratio of effective pressure to ice pressure and the logarithm of the sliding regime is 
$-$0.992 in this region. When effective pressure does not vary, as in model-IH, the sliding regime remains nearly constant (Fig. 5f), showcasing the importance of the effective pressure in the evolution of the basal shear stress. In model-CS, there is still a high negative correlation between the logarithm of the sliding regime and the ratio of the effective pressure to ice pressure (
$-$0.861), but it is not as extreme as in model-CH. In model-CS, large changes in ice speed have an impact on the sliding regime, which is in part responsible for the less significant correlation.
The linear combination friction law
The ice speedup events of Ice Stream B in model-CLC were comparable in scale to model-CH, and there was no acceleration of ice upglacier of Ice Stream B despite reduced effective pressure in this area. Three speedup events occurred in model-CLC, one in 2058, one in 2081 and one in 2098 (Fig. 5g). The largest of these events sped up the ice stream by 47.1 m a
$^{-1}$. None of these ice speedup events registered in the ice speed trend of Ice Stream D. Ice Stream D and Ice Stream E followed similar patterns of slowdown in model-ILC and model-CLC, with ice speed from both models remaining within 10 m a
$^{-1}$ of each other in these regions. Like in model-CH, the initial sliding regime of the linear combination friction law follows patterns in the initial hydrology sheet thickness in model-CLC. The sliding regime in model-CLC is in general more Coulomb everywhere than in model-CH, even though the speedup events in model-CLC are of comparable size to those of model-CH. Although the sliding regime of model-CLC is more Coulomb than model-CH, the functional forms of the hydrology-modified Schoof friction law and the linear combination friction law are different, which may explain why the speedup events are of comparable size. As in model-CH, model-CLC has a large negative correlation between the logarithm of the sliding regime and the ratio of the effective pressure to the ice pressure in Ice Stream B (Fig. 5h, 
$-$0.989), and the sliding regime remains nearly constant in model-ILC.
Melt and channelization
Channels formed under the Siple Coast ice streams and in the Byrd trough but had little effect on the dynamics observed in this study. Low hydraulic potential gradients combined with low effective pressure in the Siple Coast resulted in the existence of channels with large area (500 m
$^{2}$) but limited discharge rate (
$ \lt 40$ m
$^{3}$ s
$^{-1}$) at the start of the transient simulations. Consequently, little channel evolution was observed. A discussion of channelization from the steady-state GlaDS simulation is given in the Supplementary Material text S3.
Discussion
Here, we analyze how the basal shear stress evolves in relation to changes in basal sliding speed and effective pressure. For an effective pressure dependent sliding/friction law, the basal shear stress evolves as
\begin{equation}
d\tau_b=\frac{\partial\tau_b}{\partial u_b}du_b+\frac{\partial\tau_b}{\partial N}dN,
\end{equation}where 
$d\tau_b$ is an infinitesimal change in the basal shear stress, 
$du_b$ is an infinitesimal change in the basal ice speed and 
$dN$ is an infinitesimal change in the effective pressure. For the Budd friction law, the partial derivatives of the basal shear stress are given by
\begin{equation}
\frac{\partial\tau_b}{\partial u_b}=\frac{s\tau_b}{u_b},\quad\frac{\partial\tau_b}{\partial N}=\frac{r\tau_b}{N}.
\end{equation}For the Schoof and hydrology-modified Schoof friction laws, the partials are given by
\begin{equation}
\frac{\partial\tau_b}{\partial u_b} = \frac{m\tau_b}{u_b(1+\xi_i)},\quad\frac{\partial\tau_b}{\partial N}=\frac{\tau_b\xi_i}{N(1+\xi_i)},
\end{equation}where 
$\xi_i$ is the sliding regime of the friction law, the parameter 
$\xi_s$ for the Schoof friction law and 
$\xi_h$ for the hydrology-modified Schoof friction law, defined in the ‘Methods’ section. The partial derivatives of the basal shear stress for the linear combination friction law are given by
\begin{equation}
\frac{\partial\tau_b}{\partial u_b}=\frac{C_\text{max}N}{u_b}\frac{\xi_\ell}{(1+\xi_\ell)^{2}}+C_\ell^2u_b^{m-1}\frac{m+(m-1)\xi_\ell}{(1+\xi_\ell)^{2}},
\end{equation}
\begin{equation}
\frac{\partial\tau_b}{\partial N}=C_\text{max}\frac{\xi_\ell(1-n+\xi_\ell)}{(1+\xi_\ell)^2}+\frac{nC_\ell^2u_b^m}{N}\frac{\xi_\ell}{(1+\xi_\ell)^{2}}.
\end{equation} For all the regularized Coulomb friction laws (the Schoof friction law, the hydrology-modified Schoof friction law and the linear combination friction law), in the limit that 
$\xi_i$ goes to 0, we have
\begin{equation}
\left.\frac{\partial\tau_b}{\partial u_b}\right\rvert_{\xi_i\to0}=\frac{m\tau_b}{u_b},\quad\left.\frac{\partial\tau_b}{\partial N}\right\rvert_{\xi_i\to0}=0,
\end{equation}and in the limit that 
$\xi_i$ goes to 
$\infty$ we have
\begin{equation}
\left.\frac{\partial\tau_b}{\partial u_b}\right\rvert_{\xi_i\to\infty}=0,\quad\left.\frac{\partial\tau_b}{\partial N}\right\rvert_{\xi_i\to\infty}=\frac{\tau_b}{N}.
\end{equation}The Budd sliding law
We find that the Budd sliding law is only applicable for running ice flow models in slow moving areas, and we suggest that it is not suitable for two-way coupled modeling. For 
$r=1$ and 
$s=m=1/n$, as is commonly chosen for the Budd sliding law, the Budd sliding law responds to changes in basal sliding speed as if the entire domain were in a Weertman sliding regime, and it responds to changes in effective pressure as if the entire domain were in a Coulomb sliding regime. If running an uncoupled or one-way coupled ice flow model, the Budd sliding law would behave exactly like the Weertman sliding law, given the same initial basal shear stress. In regions that should be in a Weertman sliding regime, the Budd sliding law will sufficiently describe the evolution of the basal shear stress. However, when two-way coupled to a subglacial hydrology model, the Budd sliding law will incorrectly describe the evolution of the basal shear stress everywhere, as changes with respect to effective pressure will be incorrect in areas that should be in a Weertman sliding regime, and changes with respect to basal sliding speed will be incorrect in areas that should be in a Coulomb sliding regime. When used for two-way coupled modeling, the Budd sliding law has the ability to produce fast moving ice and low effective pressure in regions that should be in a Weertman sliding regime. This is the consequence of a positive feedback mechanism between effective pressure and ice speed. A decrease (increase) in effective pressure will result in a decrease (increase) in basal shear stress, triggering ice acceleration (deceleration), which will cause dynamic thinning (thickening) of ice, reducing (raising) the hydraulic potential, which will draw more (less) water into the area and lower (raise) the effective pressure. Such a feedback mechanism is possible when trends in ice thickness either match or are overridden by the feedback, and when changes in effective pressure have a larger impact on the basal shear stress than changes in ice speed. That is,
\begin{equation}
\left\lvert\frac{\partial\tau_b}{\partial N}dN\right\rvert \gt \left\lvert\frac{\partial\tau_b}{\partial u_b}du_b\right\rvert,
\end{equation}which, for the Budd sliding law, implies 
$dN/N \gt sdu_b/(ru_b)$. With 
$r=1$ and 
$s=m=1/n$, this feedback can occur in relatively slow moving regions (
$ \lt $100 m a
$^{-1}$), as effective pressure is typically between 0 and 20% of ice overburden pressure (92.1% of the grounded domain in model-CH). This feedback was observed upglacier of Ice Stream B in model-CB, showcasing the risks of using the Budd sliding law in two-way coupled modeling. Therefore, the Budd sliding law only agrees with our understanding of sliding over cavities when running ice flow models in slow moving areas with limited cavitation at the glacier bed, and we suggest that it is not suitable for two-way coupled modeling.
The Schoof friction law
We find that the Schoof friction law should only be used for uncoupled or one-way coupled ice flow modeling in slow flowing regions and suggest that it should not be used for two-way coupled modeling on real domains. Under the assumption that the Schoof friction law has a sliding regime that is consistent with the degree of cavitation at the glacier bed, the Schoof friction law will describe the evolution of the basal shear stress in a Weertman sliding regime in a way that is consistent with our understanding of sliding over cavities, in an uncoupled or one-way coupled ice flow model. In a Coulomb sliding regime, where the basal shear stress has a strong effective pressure dependence, two-way coupled modeling is necessary to accurately describe the evolution of the basal shear stress. However, the assumption that the sliding regime is spatially consistent with the degree of cavitation only holds if the ratio of 
$C_s^2/C_\text{max}$ is spatially uniform. Schoof (Reference Schoof2005) and Gagliardini and others (Reference Gagliardini, Cohen, Råback and Zwinger2007) modeled ice flow over periodic bedrock protrusions on scales that are much smaller than those of most ice flow modeling studies. Therefore, the resultant regularized Coulomb friction laws were proposed without taking spatially heterogeneous bed conditions (variable 
$C_s$ and 
$C_\text{max}$) into account. On the small scales considered by Schoof (Reference Schoof2005) and Gagliardini and others (Reference Gagliardini, Cohen, Råback and Zwinger2007), and sometimes when modeling synthetic domains, 
$C_s$ and 
$C_\text{max}$ are spatially uniform and this is not an issue. However, when modeling real domains, 
$C_\text{max}$ is typically chosen to be spatially uniform, while 
$C_s$ is inverted for, meaning that in practice the sliding regime of the Schoof friction law is not spatially consistent with the degree of cavitation at the glacier bed.
In slow moving regions that should be in a Weertman sliding regime (
$\tau_b\approx C_s^2u_b^m$), the chosen value of 
$C_\text{max}$ will not affect the inversion for 
$C_s$. If 
$C_\text{max}$ is chosen to be larger than it should be, the sliding regime will be underestimated, and if 
$C_\text{max}$ is chosen to be smaller than it should be, the sliding regime will be overestimated. However, in both scenarios, the value of 
$C_s$ determined by the inversion should be accurate, and the sliding regime should in theory remain much 
$\ll$1 (Weertman). In fast flowing regions that should be in a Coulomb sliding regime (
$\tau_b\approx C_\text{max}N$), the basal shear stress is controlled by the value of 
$C_\text{max}$, not 
$C_s$. However, the inversion procedure alters the basal shear stress by adjusting 
$C_s$, meaning that significant changes in 
$C_s$ are required to find a basal shear stress that matches simulated and observed ice speed. If 
$C_\text{max}$ is smaller than it should be, then 
$C_s$ is substantially increased causing a large overestimation in the sliding regime of the friction law. If 
$C_\text{max}$ is larger than it should be, then 
$C_s$ is substantially reduced, causing a large underestimation in the sliding regime of the friction law. The inversion procedure for 
$C_s$ is therefore more accurate in regions that should be in a Weertman sliding regime. In model-CS, much of the Siple Coast ice streams were taken to be significantly more Weertman than the surrounding area, and areas that should be in a Weertman sliding regime were marginally overestimated compared to model-CH (Fig. 6).
(a) is the initial sliding regime of model-CS; (b) is the initial sliding regime of model-CH; (c) is the initial sliding regime of model-CLC; and (d) is the initial hydrology sheet thickness (m) of model-CH. The pink lines represent the grounding line at the start of the simulations and the black lines represent the domain boundary.
For the Schoof and hydrology-modified Schoof friction laws, Eqn (21) implies that 
$dN/N \gt mdu_b/(\xi_iu_b)$. This means that the positive feedback mechanism between effective pressure and ice speed experienced in model-CB is only possible for a sufficiently Coulomb sliding regime when using the Schoof or hydrology-modified Schoof friction laws. In model-CS, the overestimation of the sliding regime in areas that should be Weertman triggered this positive feedback mechanism, resulting in significant acceleration of inland ice. Model-CS illustrates the importance of correctly determining the sliding regime of the Schoof friction law, and the issues that arise from the standard method of determining the basal shear stress with said friction law.
These issues regarding inaccurate estimation of the sliding regime of the Schoof friction law are present not just in two-way coupled models, but in ice flow-only models as well. However, the consequences of such inaccuracies are much more prominent in two-way coupled models. The reasons for this are threefold: (1) the basal shear stress can have at most a cube root (
$1/n$) dependence on sliding speed but can have a stronger linear dependence on effective pressure; (2) the magnitude of change in effective pressure is typically comparatively larger than the magnitude of change in the sliding speed (
$\left\lvert dN/N\right\rvert \gt \left\lvert du_b/u_b\right\rvert$); and (3) in place of the positive feedback mechanism between effective pressure, basal shear stress and ice speed, there is potential for a negative feedback mechanism between sliding speed and basal shear stress since faster sliding increases the basal shear stress but a larger basal shear stress reduces sliding speed. For these reasons, the impacts of the incorrect sliding regime of the Schoof friction law are less obvious in ice flow-only models.
Given the impact of the inversion procedure for 
$C_s$ on the sliding regime of the Schoof friction law, we suggest that the Schoof friction law should only be used for uncoupled or one-way coupled ice flow modeling in slow flowing regions with low cavitation, and it should not be used for two-way coupled modeling on real domains. This is the same conclusion we arrived at for the Budd friction law, though for different reasons. When using a Schoof friction law in an uncoupled or one-way coupled ice flow model, the basal shear stress still responds to changes in basal sliding speed based on the sliding regime of the friction law; however, the sliding regime of the friction law is not spatially consistent with the degree of cavitation at the glacier bed. In contrast, the Budd friction law responds to changes in the basal sliding speed as if the entire domain were in a completely Weertman sliding regime. The Schoof friction law may therefore be more desirable than the Budd friction law if running an uncoupled or one-way coupled ice flow model in a slow flowing region.
New friction laws
We find that the hydrology-modified Schoof and linear combination friction laws are suitable for modeling slow flowing regions with a one-way coupled or uncoupled model and can be used to model any grounded region of a glacier with a two-way coupled model. Unlike the Schoof friction law, the hydrology-modified Schoof friction law ensures that the sliding regime of the friction law is spatially consistent with the degree of cavitation at the glacier bed. Therefore, with an appropriate effective pressure input, the hydrology-modified Schoof friction law will describe the evolution of the basal shear stress in a Weertman sliding regime, in a way that is consistent with our understanding of sliding over cavities, even in a one-way coupled model. In a two-way coupled model, the hydrology-modified Schoof friction law will describe the evolution of the basal shear stress in a way that is consistent with our understanding of sliding over cavities throughout the model domain. However, by fixing the ratio of 
$C_h^2/C_\text{max}=\zeta$, the inversion for 
$C_h$ also determines Iken’s bound. That is, 
$C_h$ and 
$C_\text{max}$ cannot be chosen independently. 
$C_h$ is dependent on the ratio of bedrock bump size to spacing (Weertman, Reference Weertman1957), and 
$C_\text{max}$ is the maximum angle on the stoss face of bedrock bumps (Iken, Reference Iken1981; Schoof, Reference Schoof2005). These two ideas are physically independent of each other, yet, in the hydrology-modified Schoof friction law, Iken’s bound is determined by the inversion for 
$C_h$. In areas with low cavitation, the inversion procedure chooses 
$C_h$ such that 
$\tau_b\approx C_h^2u_b^m$, and in areas with high cavitation, the inversion procedure chooses 
$C_h$ such that 
$\tau_b\approx C_h^2N/\zeta$. Therefore, the value of 
$C_h$ determined by the inversion procedure only remains accurate during transient model simulations while the current sliding regime remains similar to the initial sliding regime that the inversion was performed on. That is, large changes in the degree of cavitation at the glacier bed (switches between Weertman and Coulomb sliding) may lead to an evolution of the basal shear stress that is not consistent with our understanding of sliding over cavities. The Kamb Ice Stream did not reactivate in the timeframe of any of our models, but if it did, the hydrology-modified Schoof friction law may be insufficient to accurately describe such a reactivation where there would be a transition from Weertman to Coulomb sliding.
The different functional form between the linear combination friction law and the Schoof/hydrology-modified Schoof friction laws means that areas with sliding regimes between Weertman and Coulomb are dealt with differently. While we know how friction laws should behave in the limit of Weertman and Coulomb sliding, it is unclear how to deal with these areas in between, and future work on the development of basal friction laws should consider these regions. The primary benefit of the linear combination friction law is having the ability to choose Weertman and Coulomb friction parameters while having a sliding regime that is spatially consistent with the degree of cavitation at the glacier bed. However, the value of 
$C_\text{max}$ was chosen to be spatially uniform in the linear combination friction law, as there is no good alternative to choosing both 
$C_\ell$ and 
$C_\text{max}$ at the moment. Future work should focus on a methodology to constrain both 
$C_\text{max}$ and 
$C_\ell$ simultaneously. Such a methodology would likely need to go beyond a snapshot inversion, which for strong Weertman and Coulomb sliding regimes will only be able to accurately constrain 
$C_\ell$ or 
$C_\text{max}$ respectively. Likely, a relationship between Weertman and Coulomb friction parameters must be identified (if one exists) through studying areas with temporally variable bed conditions, such as surging glaciers, or glaciers with seasonal velocity patterns. If such a relationship can be identified, then an inversion could be used to constrain both parameters simultaneously. Alternatively, these parameters could be constrained through other means, such as using geologic information.
Although the linear combination friction law offers desirable properties that the Budd, Schoof and hydrology-modified Schoof friction laws lack, it is not derived from first principals and is not consistent with the criteria outlined by Gagliardini and others (Reference Gagliardini, Cohen, Råback and Zwinger2007) that a friction law should follow. In particular, Iken’s bound is not strictly enforced. When 
$\xi_\ell \gt \alpha(C_\text{max}/C_\ell^2)^n$, then 
$\tau_b \gt C_\text{max}N$. For large 
$\xi_\ell$, we can Taylor expand Eqn (12) and find
\begin{equation}
\tau_b\approx C_\text{max}N+\xi_\ell^{-1+1/n}\frac{C_\ell^2}{\alpha^{1/n}C_\text{max}}-\xi_\ell^{-1}.
\end{equation} So, for sufficiently large 
$\alpha$, and large 
$\xi_\ell$ where the bound is most important, the bound will very approximately hold. Therefore, the relationship between the basal shear stress and effective pressure not strictly obeying Iken’s bound has the largest potential impact on areas intermediate to Weertman and Coulomb sliding, which have not been a focus in the development of past friction laws. The new regularized Coulomb friction laws analyzed in this study balance desirable properties from theoretical analyses of basal sliding with the practical necessities required to model real domains. Additional investigation into a more appropriate functional form of the friction law that can encompass all of these properties is necessary.
The hydrology-modified Schoof friction law achieves a unified description of subglacial conditions between the friction law and subglacial hydrology model through constricting parameter values of the friction law. An alternative approach is to keep the friction law unaltered and constrict the parameter values of the subglacial hydrology model. This is similar to the approach of Gilbert and others (Reference Gilbert, Gimbert, Thøgersen, Schuler and Kääb2022) and Togaibekov and others (Reference Togaibekov, Gimbert, Gilbert and Walpersdorf2024), who reformulated the friction law and subglacial hydrology model in terms of a dimensionless ‘cavitation state’ variable. The cavitation state was defined so that the friction law would be unaltered, and the sheet conductivity and sheet thickness were reformulated in terms of a power law dependence on the cavitation state. The cavitation state evolved due to opening from sliding and closing from ice creep, written in terms of the friction law parameters. In this study, we have instead focused on the development of friction laws that are consistent with subglacial hydrology models. However, it is really a consistent description of basal conditions between the friction law and subglacial hydrology model that is important, and this can be achieved through alternative means, as done in Gilbert and others (Reference Gilbert, Gimbert, Thøgersen, Schuler and Kääb2022), Togaibekov and others (Reference Togaibekov, Gimbert, Gilbert and Walpersdorf2024). In the context of our study, the Schoof friction law could be made to have a sliding regime that is spatially consistent with the degree of cavitation through varying 
$\ell_r$ so that it is proportional to 
$C_\text{max}/(C_s^2\tilde{A})$. However, this would create a much more difficult inversion problem to solve when modeling real domains, as varying 
$\ell_r$ will vary the effective pressure, which the basal shear stress is dependent on. Furthermore, there is no inherent advantage to choosing a value of 
$C_s$ over choosing a value of 
$\ell_r$ (as done in this study), as there are, in general, limited geophysical observations of basal conditions to constrain either parameter.
The assumption of steady-state hydrology
The hydrology-modified Schoof friction law and the linear combination friction law assume that the hydrology sheet thickness is in steady state, so that 
$u_b/N^n$ is completely determined by 
$h$. This is in general a good assumption, as an order of magnitude approximation of the hydrology sheet opening rate gives 
$\omega_s\approx10^{-8}$ m s
$^{-1}$ (
$u_b\approx10$ m a
$^{-1}$, 
$h_r-h\approx0.1$ m, 
$\ell_r\approx$ 1 m) for a Weertman sliding regime, and for a Coulomb sliding regime 
$\omega_s\approx10^{-8}$ m s
$^{-1}$ (
$u_b\approx1000$ m a
$^{-1}$, 
$h_r-h\approx0.001$ m). 
$\frac{\partial h}{\partial t}$ typically changes on the order of 
$10^{-12}$ m s
$^{-1}$ meaning that 
$\omega_s$ and 
$\nu_s$ are the dominant terms in Eqn (4). Still, in regions where there is fast evolution of the hydrology sheet thickness (
$\frac{\partial h}{\partial t}$ balances the opening or closing term in Eqn 4), Eqn (11) will no longer be applicable. On a related note, the friction law criteria proposed by Gagliardini and others (Reference Gagliardini, Cohen, Råback and Zwinger2007) is based on a steady-state assumption, and the regularized Coulomb friction laws developed based on this work are also likely not applicable in these regions.
Effective pressure dependent friction laws calculate the basal friction of a steady-state hydrology system characterized by the effective pressure. Alternatively, Thøgersen and others (Reference Thøgersen, Schuler, Gilbert and Malthe-Sørenssen2019) and Gilbert and others (Reference Gilbert, Gimbert, Thøgersen, Schuler and Kääb2022) proposed regularized Coulomb friction laws that calculate the basal friction of a steady-state hydrology system characterized by subglacial cavitation. Similarly, we can use Eqn (11) to reformulate the hydrology-modified Schoof and linear combination friction laws in terms of 
$h$. For the hydrology-modified Schoof friction law, we have
\begin{equation}
\tau_b=\frac{C_h^2\left\lvert \overrightarrow{u_b}\right\rvert^{m-1}}{(1+ \zeta^{1/m} \tilde{A}\ell_rh/(h_r-h))^m}\overrightarrow{u_b},
\end{equation}and for the linear combination friction law we have
\begin{equation}
\tau_b=C_\text{max}N_h\left\lvert \overrightarrow{u_b}\right\rvert^{-1}\overrightarrow{u_b}\left(\frac{\xi_\ell}{1+\xi_\ell}\right)+C_\ell^2\left\lvert \overrightarrow{u_b}\right\rvert^{m-1}\overrightarrow{u_b}\left(\frac{1}{1+\xi_\ell}\right),
\end{equation}
\begin{equation}
\xi_\ell=\alpha\tilde{A}\ell_r\frac{h}{h_r-h},\quad N_h=\left(\frac{u_b(h_r-h)}{\tilde{A}\ell_rh}\right)^m.
\end{equation} We refer to Eqns (23) and (24) as the ‘cavity’ formulation of the hydrology-modified Schoof and linear combination friction laws, and we refer to the original Eqns (9) and (12) as the ‘pressure’ formulation. If we define 
$h_N$ to be the steady-state sheet thickness for a given effective pressure (calculated from rearranging Eqn 11), then the pressure formulation of the hydrology-modified Schoof friction law and the linear combination friction law can be written in the form of Eqns (23) and (24), respectively, with 
$h$ replaced by 
$h_N$, and 
$N_h$ replaced with 
$N$ for the linear combination friction law. Furthermore, for an evolving hydrology sheet, 
$h$ is transitioning toward 
$h_N$. So, the cavity formulation can be thought of as describing basal friction of the past system, and the pressure formulation can be thought of as describing basal friction of where the system is headed. Therefore, basal friction should be intermediate to that predicted from the cavity formulation and that predicted from the pressure formulation. This can be described using Eqns (23) and (24) with 
$\tilde{h}=\lambda h+(1-\lambda)h_N$ in place of 
$h$, where 
$\lambda\in[0,1]$ is a parameter that controls the importance of where the system was vs where the system is headed. For the hydrology-modified Schoof friction law, this gives
\begin{equation}
\tau_b=\frac{C_h^2\left\lvert \overrightarrow{u_b}\right\rvert^{m-1}}{(1+\zeta^{1/m} \tilde{A}\ell_r\tilde{h}/(h_r-\tilde{h}))^m}\overrightarrow{u_b},
\end{equation}and for the linear combination friction law this gives
\begin{equation}
\tau_b=C_\text{max}\tilde{N}\left\lvert \overrightarrow{u_b}\right\rvert^{-1}\overrightarrow{u_b}\left(\frac{\xi_\ell}{1+\xi_\ell}\right)+C_\ell^2\left\lvert \overrightarrow{u_b}\right\rvert^{m-1}\overrightarrow{u_b}\left(\frac{1}{1+\xi_\ell}\right),
\end{equation}
\begin{equation}
\xi_\ell=\alpha \tilde{A}\ell_r\frac{\tilde{h}}{h_r-\tilde{h}},\quad \tilde{N}=\left(\frac{u_b(h_r-\tilde{h})}{\tilde{A}\ell_r\tilde{h}}\right)^m.
\end{equation} If 
$\lambda=0$, we recover the pressure formulation, and if 
$\lambda=1$, we recover the cavity formulation. Equations (26) and (27) also provide a sliding regime that is intermediate to that of the cavity formulation and the pressure formulation, which is also likely more applicable in regions where there is a rapidly evolving hydrology sheet.
Uncoupled ice flow models
In ice flow modeling studies, the effective pressure is often not determined by running a subglacial hydrology model, with an effective pressure parameterization being used in the basal friction law instead. These parameterizations are often poor estimators of the effective pressure and are in some cases nonphysical (McArthur and others, Reference McArthur, McCormack and Dow2023). The standard parameterization used is actually the overburden hydraulic potential with the density of freshwater replaced with the density of salt water, 
$N=\phi_0=\rho_igH+\rho_{sw}gB$. This parameterization assumes perfect hydraulic connectivity to the ocean (subglacial water pressure is equal to hydrostatic pressure at the grounding line) but leads to unrealistically large values of effective pressure away from the grounding line. In fact, when the bed elevation is above sea level, this parameterization implies negative subglacial water pressure, which is not physical. This large overestimation of effective pressure will lead to a large underestimation in the sliding regime of the friction law. McArthur and others (Reference McArthur, McCormack and Dow2023) proposed an empirical parameterization of the effective pressure, motivated from subglacial hydrology modeling of Denman Glacier, East Antarctica. This parameterization is based on the idea that the flotation fraction (fraction of subglacial water pressure relative to ice overburden pressure) decreases as the ice thickness increases, an idea that does not work well for topographically unconstrained ice streams such as those in the Siple Coast. Still, this empirical parameterization ensures the effective pressure does not reach nonphysically large values (larger than ice overburden pressure) and resulting significant underestimations of the sliding regime of the friction law. Two-way coupled models are computationally more expensive than uncoupled ice flow models and have added modeling complexities (such as complicated initialization procedures) that will make two-way coupled modeling undesirable in some cases. In these situations, improved effective pressure parameterizations are necessary to describe basal conditions in fast flowing Coulomb sliding regimes, which are often of interest.
The impact of depth averaging on basal friction
In this study, we used a 2-D ice flow model, with depth-averaged ice velocity. Therefore, the basal ice speed is likely overestimated, particularly in areas where there should be Weertman sliding, and ice motion is primarily a result of internal deformation. An overestimation of the ice velocity implies there is an underestimation in the basal friction coefficient, so that the same basal shear stress is achieved. In a strong Weertman sliding regime, where an increase in ice surface speed may not reflect the onset of basal sliding, there should be an overestimation of basal friction in our models. However, in a transition from Weertman to Coulomb sliding, where Iken’s bound becomes the relevant friction parameter, friction is likely underestimated in the hydrology-modified Schoof friction law, where 
$C_\text{max}=C_h^2/\zeta$. In the Schoof and linear combination friction laws, this should be less of an issue as Iken’s bound is chosen independently of the Weertman friction parameter. Due to the dependence of the friction law on basal ice motion, a 3-D ice flow model is desirable for modeling slow moving ice in a Weertman sliding regime. However, even a 2-D model with a vertically heterogeneous ice velocity such as the MOno-Layer Higher-Order (Dias dos Santos and others, Reference Dias dos Santos, Morlighem and Brinkerhoff2022) model is preferable over the SSA. In a Coulomb sliding regime, a 2-D ice flow model with the linear combination friction law should more appropriately describe the evolution of the basal shear stress compared to a Weertman sliding regime.
Conclusions
In this study, we investigated the impact of the friction law on ice flow modeling using one-way and two-way coupled ice flow/subglacial hydrology models in the Siple Coast as a case study. Additionally, we presented two new regularized Coulomb friction laws, developed specifically for two-way coupled modeling.
It is clear from the results of this study that the present implementation of basal friction in ice flow modeling leaves room for improvement. We have found that one-way coupled or uncoupled ice flow models can sufficiently describe ice motion in a Weertman sliding regime. However, fast moving, highly cavitated areas that are in a Coulomb sliding regime require the use of two-way coupled modeling where the evolution of the basal shear stress is dominated by changes in effective pressure to capture the physical processes present in those regions. The Budd sliding law and the Schoof friction law, which are commonly used in ice flow modeling studies, were found to be unsuitable for two-way coupled modeling. When using the Budd sliding law, the basal shear stress responded to changes in ice speed as if the entire domain were in a Weertman sliding regime and responded to changes in effective pressure as if the entire domain were in a Coulomb sliding regime. The latter resulted in a positive feedback mechanism that significantly increased the speed of inland ice while simultaneously flooding the subglacial hydrologic system. The sliding regime of regularized Coulomb friction laws was found to be an important factor, which controls the evolution of the basal shear stress. The Schoof friction law has a sliding regime that is not dependent on the degree of cavitation at the glacier bed in a spatially consistent fashion, which led to the same feedback mechanism experienced with the Budd friction law.
In light of the limitations of the Budd and Schoof friction laws, we proposed and tested two additional regularized Coulomb friction laws that have a sliding regime that is spatially consistent with the degree of cavitation at the glacier bed. The first is a version of the Schoof friction law that we modified so that 
$C_s^2/C_\text{max}$ is spatially uniform, called the ‘hydrology-modified Schoof friction law’, and the second is a linear combination of Weertman and Coulomb sliding terms with coefficients that sum to one and are dependent on the sliding regime of the friction law, called the ‘linear combination friction law’. These friction laws did not lead to nonphysically large changes in ice motion despite a similar reduction of effective pressure upglacier of Ice Stream B, albeit on a much smaller scale. Furthermore, the linear combination friction law is able to choose Coulomb and Weertman friction parameters (
$C_\ell$ and 
$C_\text{max}$), while simultaneously having a sliding regime that is spatially consistent with the degree of cavitation at the glacier bed. However, it still requires the modeler to choose the value of Iken’s bound, which introduces additional uncertainties to the model.
We have identified the importance of including evolving effective pressure in basal friction parameterizations, having a sliding regime that is spatially consistent with the degree of cavitation at the glacier bed and having the ability to choose both Coulomb and Weertman friction parameters. However, additional work is necessary to describe the evolution of basal friction in ice flow modeling in a way that is consistent with our understanding of sliding over cavities. In particular, a procedure must be developed to choose Weertman and Coulomb friction parameters simultaneously, the ideal functional form of the friction law must be investigated, and improved effective pressure parameterizations must be developed and implemented in uncoupled ice flow models.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10216.
Data availability statement
We use version 4.24 of the open-source ISSM software; ISSM is available for download from https://github.com/ISSMteam/ISSM. The datasets used in the initialization procedure are available at the references cited in the main text of this manuscript. Model outputs and scripts for running the models can be found at McArthur and others (Reference McArthur, Dow and Ehrenfeucht2026b).
Acknowledgements
K.R.K.M. was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), Canada Graduate Scholarship and the Ontario Graduate Scholarship. C.F.D. was supported by NSERC (NSERC RGPIN-03761-2017 and RGPIN-2025-03982) and by the Canada Research Chairs Program. S.E. was supported by NSERC (NSERC RGPIN-03761-2017). The authors would like to thank the Digital Research Alliance of Canada for access to supercomputing resources. The authors would like to thank the scientific editor Alan Rempel and three anonymous reviewers for their suggestions to improve this manuscript.
Open access
