2.1 Hydrology model

The hydrology model is based on the concept that a certain amount of water gets stored in the sediments underlying the ice sheets, and, once saturated, the excess is transported in the direction of the hydrological gradient to the ice margin. Some components of our model derive from the routing scheme described by Bueler and van Pelt (Reference Bueler and van Pelt2015), but we have simplified the implementation to emphasize computation speed. Our model does not conserve mass, and transports water to the edge of the ice sheet without any time delay. When the routing of the excess water is computed, the water from upstream grid cells is added to each grid cell downstream. This is not entirely realistic, since the hydrological system can react at a time scale on the order of hours (Bartholomew and others, Reference Bartholomew2012). The time stepping in the model is usually on the order of days to months, so this simplification may be considered to be representative of average conditions. Ultimately, the output of the hydrology model is the effective pressure at the base of the ice sheet, which is then transferred into the basal sliding model. A schematic of the components of our model is shown on Figure 1, while Figure 2 shows the workflow of the model.

Fig. 1. Schematic of the components of the new basal conditions model. (a) Overview of ice sheet hydrology. (b) Overview of impact on sliding.

Fig. 2. Diagram showing the workflow of the model.

2.1.1 Water routing

The first component of the model is that it captures the surface melt. We are using the semi-analytical positive degree day (PDD) method module (Calov and Greve, Reference Calov and Greve2005). As implemented in PISM, it computes the amount of ice that melts at the surface as a diagnostic parameter. Our modification stores this value and passes it to our hydrology model. Within our model, there is an option to set the fraction of the meltwater that gets transferred to the base of the ice sheet (see Table 1 for a full list of command line options available for the model). The water transferred from the surface is added to the meltwater generated from heating at the base (Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012).

Table 1. Command line options available for the described models

The next step is a modification of the undrained plastic bed model (Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000; Bueler and Brown, Reference Bueler and Brown2009). In this model, a layer of sediment of a specified thickness and porosity fills with water until it is saturated, which is set within PISM as a ‘water thickness’ parameter, W _sed. The saturation, s, is:

(1)$$s = {W_{sed}\over W_{sed}^{\max}}.$$

W _sed is the amount of water in the sediments, represented as a layer below of the ice sheet that fills when there is water input into the subglacial hydrology system, while $W_{sed}^{\max }$ is the maximum thickness of that layer. In our simulations, $W_{sed}^{\max } = 1$ m, the value used in Niu and others (Reference Niu, Lohmann, Hinck, Gowan and Krebs-Kanzow2019). If the porosity of a deforming till is 40% (Blankenship and others, Reference Blankenship, Bentley, Rooney and Alley1987), this value implies that 2.5 m layer of subglacial sediment is active in the hydrology system of the ice sheet.

A certain amount of accumulated water within the sediment is removed at every time step in order to simulate drainage. This is determined by reducing the basal water layer thickness at a rate of 1  mm yr⁻¹, which is the default value in PISM. At every grid cell where s < 1, any subglacial water will be added to the sediments. Our modification from the default model is that the amount of water that can enter the sediments depends on the fraction of the subglacial surface that is covered in sediment. If the sediment cover is incomplete, then the sediments can fill with water to the maximum level faster than if there is complete cover since there is less sediment to accommodate the water. Note that our model does not take into account the possibility that the underlying sediments are impenetrable due to being frozen. The consequence of this is that the water flux is underestimated (since water would not be able to enter the sediments) and the area where sediment deformation happens would be overestimated (since frozen sediments cannot deform).

Any excess water in the grid cell after filling the sediments is transported to the edge of the ice sheet. We use a simple subglacial water routing routine, where the water is transported in the direction opposite to the hydrological potential gradient, $\nabla \phi _h$. Note that the routing of water happens after the sediment filling step, so none of the water added to a grid cell from upstream contributes to the water in the sediments. The equation for calculating the potential gradient at the base of the ice sheet is as in Cuffey and Paterson (Reference Cuffey and Paterson2010):

(2)$$-\nabla \phi_h = -\rho_i g \left[\,f_w \nabla S + \left[{\rho_w\over \rho_i} -f_w \right]\nabla B \right].$$

In this equation, ρ _i is the ice density, ρ _w is the water density, g is the gravitational acceleration, $\nabla S$ is the ice surface gradient, $\nabla B$ is the bed gradient, and f _w is the flotation fraction, which is the ratio of the water pressure and overburden pressure. The flotation fraction governs the relative influence of the bed and ice surface slopes on the direction of water flow. We have set it to be a constant, f _w = 0.8, which gives the surface slope a 2.7 times greater influence on the routing (Cuffey and Paterson, Reference Cuffey and Paterson2010). This ensures that the water will generally move toward the edge of the ice sheet. We calculate the gradient either using a third order finite difference method described in Skidmore (Reference Skidmore1989) or using a least squares method on a 5 × 5 grid (i.e. all the grid cells within 2 cells of cell where the gradient is calculated), the later which is the default.

The routing of water is accomplished by first sorting the hydrological potential, ϕ _h, values over the entire grid from highest to lowest. In order to avoid singularities, S and B are smoothed using a 5 × 5 average filter. The potential values are sorted from highest to lowest and the water is routed in the direction opposite to the gradient. This results in increasing amounts of water toward the edge of the ice sheet, where the potential will be the lowest as the ice sheet is thinnest. If the gradient of a cell is below a certain threshold (which we have set to be 1.0 N m⁻³), then no water is distributed as it is assumed that the water would be flowing too slowly to be distributed. The amount of water determined to go through each grid cell, which we define as T _w, is used to determine the effective pressure, which is described in the next step. The amount of water, T _w is expressed as the thickness of water over the entire grid cell per unit time (i.e. units of m s⁻¹).

2.1.2 Effective pressure

To calculate the effective pressure, we use a parameterization described by Schoof (Reference Schoof2010). This parameterization is based on the concept of water drainage at the bottom of the ice sheet being routed through efficient Röthlisberger channels (Röthlisberger, Reference Röthlisberger1972) or less efficient linked cavities (Kamb, Reference Kamb1987). This is a modification of other subglacial drainage models that have been proposed in the past (Fowler, Reference Fowler1987; Hewitt and Fowler, Reference Hewitt and Fowler2008), but allows for better switching between drainage styles. The style of drainage system is dependent on the amount of water available and the velocity of the ice. In this formulation, the effective pressure decreases up to a certain point, after which drainage becomes efficient enough that it causes the effective pressure to increase again.

The main component of this model is the switch between channel and cavity drainage systems. The type of drainage system is dependent on the total water flux, Q. The threshold water flux, Q _c is calculated by the following equation:

(3)$$Q_c = {u_b k\over c_1 ( \alpha - 1 ) \nabla \phi_h}.$$

The velocity of the ice at the base is u _b. In this model, the bed is assumed to be rough, with a protrusion height of k, which we have set to be 0.1 m. The constant c ₁ is related to the latent heat of fusion of ice, L, and is calculated by c ₁ = 1/(ρ _i L). The constant α = 5/4 is related to the Darcy–Weisbach equation friction factor for water flow in a conduit (Schoof, Reference Schoof2010). This quantity is calculated purely as a diagnostic value, and does not influence the modeled drainage system.

For the parameterization from Schoof (Reference Schoof2010), the effective pressure is calculated using the assumption that the water discharge is in steady state. This assumption reduces the complexity of the effective pressure calculation, since it does not have dependence on the size of the conduits or the style of drainage system. We use the total amount of water going through each cell, T _w, as calculated in the previous section, to determine the water flux. The water is assumed to be directed through a single channel. The total flux of water through a channel, Q, considering a grid cell of width dx is calculated as follows:

(4)$$Q = {T_w {\rm d}x^2\over {\rm d}x / r}.$$

The value of r is the spacing between channels. This value is set to a constant of 12 km, which is the average distance between eskers on the Canadian Shield (Storrar and others, Reference Storrar, Stokes and Evans2014). This formulation allows for the proper parameterization of water flux through the channel regardless of the actual width of the grid cell. Based on Eq. (3), if Q > Q _c, then the routing is via the tunnel system (efficient drainage), while if Q < Q _c, the drainage is via a cavity system (inefficient drainage). As a result of this formulation, if the ice velocity increases, the threshold amount of discharge to switch to efficient tunnel drainage also increases.

The effective pressure in the drainage system, N _hyd is calculated by the following equation (Schoof, Reference Schoof2010):

(5)$$N_{hyd}^n = {c_1 Q \nabla \phi_h + u_b h \over c_2 c_3^{-1/\alpha} Q^{1/\alpha} {\nabla \phi_h}^{-1 / ( 2 \alpha) } }.$$

The exponent, n is the Glen exponent, which by default is 3. The thickness of the ice is h. The velocity of the ice at the base is u _b. The constant c ₂ = 2 A n ⁻ⁿ includes parameters in Glen's law, where A is the ice softness. The default value is A = 3.1689 × 10²⁴ Pa⁻³  s⁻¹ (Huybrechts and Payne, Reference Huybrechts and Payne1996). The constant c ₃ is related to the relation for turbulent flow of water in the Darcy–Weisbach equation, where $c_3 = 2^{1/4} \sqrt {\pi + 2} / [ \pi ^{1/4} \sqrt {\rho _w f}]$, and f is a friction factor. We use the value f = 0.1 (Schoof, Reference Schoof2010).

There is a check so that the calculated effective pressure is not greater than the overburden pressure:

(6)$$N_{hyd} \leq {\rho}_i \,g h.$$

If the effective pressure is greater than the threshold, it is set to be equal to the overburden pressure. This can happen if the water flux, Q, becomes much smaller than u _b h. There is also a check to ensure that the effective pressure is greater than a minimum threshold, which we have set to be 0.01 times the overburden pressure. This should rarely happen except if the equation is solved where there is essentially no ice, or where there is no surface gradient and velocity. In reality, negative effective pressures can exist under glaciers when there is a rapid influx of water, which can cause the ice to temporarily float (Roberts, Reference Roberts2005). Here we have removed the possibility of negative effective pressure only to ensure the stability of the ice sheet model.

2.2 Basal sliding model

The sliding model that we use is basically a modification of the existing Mohr–Coulomb yield stress relationship that is generally used as the sliding law in PISM (Bueler and van Pelt, Reference Bueler and van Pelt2015). The general definition for the Mohr–Coulomb yield stress, τ _c, is a function of the effective pressure, N, the angle of internal friction, ϕ, and a cohesion parameter, c.

(7)$${\tau}_c = N \tan( \phi) + c.$$

The value of ϕ determines the angle that the material will fail if a normal stress is applied. In the default PISM sliding law, the entire base of the ice sheet is assumed to be covered in a layer of deformable sediments (i.e. soft bedded sliding), and ϕ is the shear friction angle of the sediments. For sediments, this value will depend on the dominant grain size, with clay materials having a lower value than sand and gravel. When a sediment under the ice sheet becomes water saturated, the effective pressure decreases, which increases the chance of failure. In general, the cohesion is regarded as being negligible in a deforming till (Cuffey and Paterson, Reference Cuffey and Paterson2010), so it is set to c = 0.

In PISM, the basal shear stress, τ _b that balances the driving stress is related to the yield stress τ _c by (Bueler and Brown, Reference Bueler and Brown2009):

(8)$${\tau}_{b, i, j} = -{\tau}_c {v_{i, j}\over ( v_1^2 + v_2^2) ^{{1}/{2}}}.$$

In this equation, v is the basal ice velocity, and the indices i, j refer to the directional components of the velocity. In PISM, the SSA is used to compute the stress balance only when v > 0, otherwise the non-sliding shallow ice approximation (SIA) is used. The value of τ _c used in the modified model is described below.

The modified sliding law has two components, sliding due to the deformation of saturated sediments, and sliding due to the interactions between the water in the drainage system and the ice-bed interface. The sliding between the ice and the substrate when the effective pressure is low due to high water pressures is considered to be analogous to a landslide, and can also be described using the Mohr–Coulomb relationship (Cuffey and Paterson, Reference Cuffey and Paterson2010). The PISM module we have created solves for both of these sliding mechanisms, and will chose the one that has a lower yield stress.

Our modified sliding law allows for spatially variable sediment cover, as places such as the Canadian Shield in North America did not have complete sediment cover (i.e. hard bedded sliding) (Fulton, Reference Fulton1995). This sliding law still allows for sediment deformation as utilized in the default PISM sliding law, and for sliding at the ice-bed interface. In this sliding law, the strength of the bed is calculated for both sediment deformation and sliding along the bed-ice interface, and the lower value is used.

First, we describe the yield stress due only to sediment deformation. The fraction of the area that is covered in sediment, S _f, can be spatially variable. This affects both components of the basal sliding model. For areas that have incomplete sediment cover, sediment deformation only happens for the fraction of the surface that has sediment, while the rest of the area is set to have a yield stress that is equal a user adjustable value, τ _bare (by default it is set to 100 kPa, which is a typical value of the yield stress at the base of a glacier (Cuffey and Paterson, Reference Cuffey and Paterson2010)). In reality this value will depend on factors such as the roughness of the bed, the debris content of the ice, and temperature, all factors that we do not estimate. This value is also set to be the maximum yield stress in sediment covered areas. The overall yield stress, τ _def, is:

(9)$$\tau_{def} = S_f N_{sed} \tan{\phi_{sed}} + ( 1 - S_f) \tau_{bare},\; $$

where τ _sed = N _sedtan (ϕ _sed) is the yield stress of the sediments. The result of this is that areas with incomplete sediment cover will be less likely to be influenced by sediment deformation as the primary mode of sliding. We sometimes refer to S _f as the percentage cover, although in these equations and in the code, it is expressed as a fraction. The effective pressure in the sediments, N _sed is the same as described by Bueler and van Pelt (Reference Bueler and van Pelt2015):

(10)$$N_{sed} = N_o \left({\delta P_o \over N_o} \right)^s 10^{ \left({e_0}/{C_c}\right)\left(1-s \right)}.$$

This equation has several constants, which in PISM are derived from Tulaczyk and others (Reference Tulaczyk, Kamb and Engelhardt2000). N _o = 1000 Pa is the reference effective pressure. e ₀ = 0.69 is the void ratio at the reference pressure. C _c = 0.12 is the compressibility of the sediments, which for this value refers to glacial till. P _o is the overburden pressure. The value s is the water saturation of the sediments, which is taken from the hydrology model described above.

For the second component of the sliding law with sliding along the ice-bed interface, the Mohr–Coulomb relationship is also used. In this case the ϕ value is related to the roughness of the interface between the ice and the bed (Iken, Reference Iken1981; Cuffey and Paterson, Reference Cuffey and Paterson2010). For clarity, we define the angle in this component as γ. A Coulomb-style law has been found to be sufficient to describe hard bedded sliding (Helanow and others, Reference Helanow, Iverson, Woodard and Zoet2021). In this model, the base of the ice sheet is covered by bumps, with an upslope angle that is equal to γ. There is a separate value for sediment covered areas (γ _sc) and areas where the bed is rock (γ _rc), as it is assumed that sediment covered areas will be smoother. First the model checks if the yield stress over sediment covered areas, N _hydtanγ _sc, is lower than the value from sediment deformation, N _sedtanϕ _sed. This lower value is taken as the yield stress over sediment covered areas and denoted as τ _sedfrac.

(11)$$\tau_{sedfrac} = {\text {min}}( N_{hyd} \tan{\gamma_{sc}},\; N_{sed} \tan{\phi_{sed}} ) .$$

As a result, if sediment cover is almost complete and the sediments are saturated with water, the effective yield stress will be similar to the default sliding law of PISM. The overall yield stress, τ _slide is:

(12)$$\tau_{slide} = S_f \tau_{sedfrac} + ( 1 - S_f) N_{hyd} \tan{\gamma_{rc}}.$$

The values of γ for sediment covered and bare areas can be set by the user. The effective pressure, N _hyd, is taken from the hydrology submodel described in the previous section.

After the yield stress for both sediment deformation and sliding at the base has been calculated, the lower of the two values is chosen as the yield stress for calculating sliding.

(13)$$\tau_c = {\text {min}}( {\tau}_{slide},\; {\tau}_{def}) .$$

In addition to sediment deformation and the combined hydrology and sediment deformation methods to find the yield stress, there is also an optional method to artificially impose a low value at the grounding line of the ice sheet when it is beside an ice shelf (the option is called ‘slippery grounding lines’ – sgl). The slippery grounding line option will ensure that the sediments are completely saturated. There is also an additional value of the potential yield stress, τ _sgl calculated, which scales the overburden pressure using the following relationship:

(14)$$\eqalignno{&\tau_{sgl} = ( 10^{-5}b + 0.2) gh\rho_i & ( b> -1000 \ {\text m}) ,\; \cr &\tau_{sgl} = ( 10^{-6}b + 0.019) gh\rho_i & ( -1000 \ {\text m} > b> -2000 \ {\text m}) ,\; \cr &\tau_{sgl} = 0.001gh\rho_i & ( b< -2000 \ {\text m}) ,\; \cr}( 14 ) $$

where b is the elevation, g is the gravitational acceleration and ρ _i is the density of ice. These equations scale the yield stress to be 0.1 to 0.2 times the overburden pressure above −1000 m, between 0.1 and 0.001 times the overburden pressure between −1000 and −2000 m, and 0.001 times the overburden pressure lower than −2000 m. This scaling will allow a reduction in the yield stress at the grounding line even where there is incomplete sediment coverage. The value of τ _sgl is used if it is lower than τ _slide and τ _def.

2.3 Limitations

In reality, if there was enough water under the ice sheet, it would cause the ice sheet to float (i.e. the water pressure would exceed the overburden pressure and N _hyd would be negative) (Schoof and others, Reference Schoof, Hewitt and Werder2012). The model does not take into account this possibility, and as a result limits the seasonal acceleration of the ice sheet. Another issue is the lack of water storage underneath the ice sheet. If there is a localized hydrological potential low point within the ice sheet, water would be routed toward this point. If enough water were to collect at such a point, it is likely that a subglacial lake would form. Since the influence of bed topography is reduced using the variable f _w, this problem is reduced, as the ice surface generally decreases toward the edge of the ice sheet. A future addition to this model that would make it more realistic would be to incorporate water conservation between time steps. This could be used to determine if a subglacial lake would form. The consequence of these limitations is that the modeled velocity of the ice sheet will be slower than reality, as a subglacial lake would essentially remove the resistance to flow at the base (Thoma and others, Reference Thoma, Grosfeld, Mayer and Pattyn2012). As an example, the ice velocity over a well studied subglacial lake in the Whillans and Mercer ice streams in Antarctica increased by up to 4% when it filled (Siegfried and others, Reference Siegfried, Fricker, Carter and Tulaczyk2016).

Our model does not take into account the possibility of spatially variable sediment thickness beyond having the possibility of having sediment free areas. This is not seen as being a major limitation, because sediment deformation mostly happens in the uppermost one meter of sediment (Boulton and others, Reference Boulton, Dobbie and Zatsepin2001). A larger possible consequence would be on the volume of water that could be stored subglacially in the sediments. Given the time scales of glacial cycle models, we regard this as a minor issue, as we expect that the aquifers would remain close to being full if water was consistently reaching the bed.

The effective pressure calculation uses an assumption that the water flux is in steady state. The style of drainage is implicit in the equation, and does not evolve if the flux is no longer in steady state. In reality, the drainage system does not necessarily switch back to an earlier state if there is a reduction of flux (Schoof, Reference Schoof2010). A more accurate drainage model would require explicit determination of the evolution of the geometry of the channels or tunnels.

Another limitation is the spatial resolution of the model simulation. The way the model is set up, it is assumed that the water is distributed to the adjacent cells. In a higher resolution model run, the pathway the water takes may become more focused than in the coarser tests that we have run. This would result in some pathways having a much higher water flux, while some adjacent cells would be much lower. The increased influence of the hydrological component of the model would be competing against adjacent cells that might not undergo a seasonal reduction in yield stress.