BrAHMs V1.0: a fast, physically based subglacial hydrology model for continental-scale application

We present BrAHMs (BAsal Hydrology Model): a physically based basal hydrology model which represents water flow using Darcian flow in the distributed drainage regime and a fast down-gradient solver in the channelized regime. Switching from distributed to channelized drainage occurs when appropriate flow conditions are met. The model is designed for long-term integrations of continental ice sheets. The Darcian flow is simulated with a robust combination of the Heun and leapfrog–trapezoidal predictor–corrector schemes. These numerical schemes are applied to a set of flux-conserving equations cast over a staggered grid with water thickness at the centres and fluxes defined at the interface. Basal conditions (e.g., till thickness, hydraulic conductivity) are parameterized so the model is adaptable to a variety of ice sheets. Given the intended scales, basal water pressure is limited to ice overburden pressure, and dynamic time stepping is used to ensure that the Courant–Friedrichs–Lewy (CFL) condition is met for numerical stability.

The model is validated with a synthetic ice sheet geometry and different bed topographies to test basic water flow properties and mass conservation. Synthetic ice sheet tests show that the model behaves as expected with water flowing down gradient, forming lakes in a potential well or reaching a terminus and exiting the ice sheet. Channel formation occurs periodically over different sections of the ice sheet and, when extensive, displays the arborescent configuration expected of Röthlisberger channels. The model is also shown to be stable under high-frequency oscillatory meltwater inputs.

1 Introduction

Subglacial basal hydrology is a potentially critical control on basal drag and ice streaming. Furthermore, it is a clear control for subglacial sediment production/transport/deposition processes (Benn and Evans, 2010; Melanson, 2012). Subglacial water flows can also leave clear geological imprints. For instance, eskers are a geological footprint of past channelized subglacial drainage (Benn and Evans, 2010) that can in turn be used to better constrain past ice sheet evolution.

Many models relating to basal hydrology are either meant for short timescales (e.g., on the order of weeks to centuries), or are missing a key component of basal water flow (channelized flow). We present a computationally fast physics-based subglacial hydrology model for continental-scale ice sheet systems modelling over glacial cycles, which is meant to capture the relevant features of basal water flow for the above three contexts (including both distributed and channelized flow components).

This large spatio-temporal scale context places a high requirement on computational speed and justifies certain simplifications compared to glacier-scale models (Bartholomaus et al., 2011; Werder et al., 2013; de Fleurian et al., 2016). Glacial cycle models do not resolve daily or even weekly mean changes in basal drag and spatial scales are relatively coarse (10–50 km). As such, the detailed physics of cavity evolution and tunnel formation cannot be resolved (given their dependence on basal sliding velocities) nor, we posit, need they be resolved. The latter is justified on the large space–time scale difference between cavities and model grid. Furthermore, the lack of adequate constraint data for this scale dictates a more simplified approach to minimize the number of tunable parameters.

Only a few subglacial hydrology models have been described in the literature for continental-scale ice sheets. Of these models, some of the more advanced include the models developed by Flowers (2000, 2008), Johnson (2002), Arnold and Sharp (2002), Goeller et al. (2013), Bueler and van Pelt (2015), and Gudlaugsson et al. (2017). These models take various approaches to simulate the flow of basal water using physically based equations.

The original work of Flowers (2000) developed a physics-based, multi-component model that included englacial, subglacial, and groundwater (aquifer) hydrology. The subglacial component of this model simulated the flow of water as a distributed system via Darcian flux. The equations were cast in a finite-volume discretization (Patankar, 1980) and advanced in time using an iterative Newton–Krylov method on a 40×40 m grid. Later work on that model included a channelized flow that coexisted alongside the distributed system and allowed exchange between the two systems (Flowers, 2008).

Johnson (2002) developed a continental-scale model with a 5 km grid resolution. In this model, the water is transported underneath the ice sheet via a tunnel (channelized) system solved using the turbulent Manning pipe flow equation. The aquifer in this model was simply a parameter that drained a percentage of the available water in the grid cell. The equations of Johnson (2002) were solved using the Galerkin method for finite-element discretization.

The work of Arnold and Sharp (2002) attempts to model the flow of water with both distributed and channelized systems. The model determines the type of system operating in each grid cell based on the water flux in the grid cell. When the flux allows the grid cell to exceed the “orifice stability parameter” (Kamb, 1987), then the grid cell has a channelized system, otherwise, it is a distributed system. The model integrates the basal water fluxes down the hydraulic potential. From the fluxes, the model determines the drainage system present, which is a different method employed from those used in the previous two models and the one developed herein.

Goeller et al. (2013) considers a distributed system that covers the base of the ice sheet. As a simplification, the basal water pressure is assumed to be approximately equal to the ice overburden pressure. This simplifies the hydraulic gradient to follow the bed and ice geometries. Water flux out of a grid cell is limited in the case that it would lead to negative water by applying a multiplier to the out-fluxes that lower their values to the desired limit. Their model does not consider any channelized system. A similar model was used in Gudlaugsson et al. (2017) applied to the Eurasian ice sheets that covered northern Europe and parts of Asia during the last ice age.

Carter et al. (2017) created a 1-D model for the simulation of lake drainage beneath Antarctica. Their model did a detailed comparison to the rate and frequency of water drainage from a subglacial lake via R-channels and canals cut into the underlying sediment. Their model showed that the canal drainage system provided better estimations for cold ice, whereas R-channels would be more common in warm, temperate ice, such as near the terminus of Greenland that is also fed by surface run-off. Similar conclusions were drawn from Dow et al. (2015).

The basal hydrology model described here combines features from the above models to create a relatively fast subglacial hydrology model for continental-scale contexts. Following the work of Arnold and Sharp (2002), the basal drainage system is allowed to have both distributed and channelized drainage systems with a condition for determining which basal system is present. While conceptually similar, the implementation is rather different. In this model, the drainage system is initially assumed to be distributed, as in Flowers (2000), with basal fluxes computed under the same Darcy flow approximation. The distributed system in a grid cell is switched to a channelized drainage system when the flux exceeds a critical value developed in Schoof (2010). The switching condition explained in Schoof (2010) is for an R-Channel, but the model can allow for other conditions to be used that better suit other channelized drainage types. Starting at the grid cells that meet the switching condition, channelized systems are created by following the path of the steepest hydraulic gradient until a potential well or exit is reached. R-channel drainage is imposed instantaneously. For developmental expediency, the aquifer physics of Flowers (2000) are replaced with the drainage parameter from Johnson (2002).

Bueler and van Pelt (2015) developed a subglacial hydrology model for the Parallel Ice Sheet Model (PISM). Similar to BraHMs, their hydrology model simulates the subglacial water flow using a Darcian flux and limits the basal pressure to the ice overburden pressure (due to long timescales). Their model consists of several basal components, including a water-filled till layer and an effective cavity-based water storage. The model presented by Bueler and van Pelt (2015) does not have any channelized flow mechanisms, which is a major source of water flow/drainage beneath ice sheets, as discussed in Sect. 2. It is unclear as to how well this model compares to BrAHMs in terms of speed due to the vast difference in model grid and computer usage. The model of Bueler and van Pelt (2015) incorporates the opening and closure of cavities which is necessary for high-resolution modelling of present-day glaciers and ice sheets, but can be replaced by low-resolution physics for longer timescales and larger spatial scales, where data is sparse. The incorporation of cavity opening and closure would require computation resources that may be prohibitive for long-term, continental-scale models.

Calov et al. (2018) used the SICOPOLIS ice sheet model in their study on the future sea level contributions of the Greenland ice sheet. The basal hydrology model used in SICOPOLIS is for large-scale grid cells like BrAHMS. The model assumes a thin film of water, resulting in zero effective pressure (meaning the hydraulic potential is simply related to bed elevation and ice sheet thickness). The model first determines a down gradient path from each grid cell to the ocean/boundary. Any depressions are filled with water (akin to lakes) and are given a small gradient so that the down gradient solver can continue. From this the water level can be calculated based on the input of meltwater and the hydraulic gradient. This is a rather different approach to BrAHMs as BraHMs attempts to model the physical evolution of the water, allowing for varied flow of water, non-zero effective pressures, and the time evolution of lake growth under Darcy flow.

Hoffman and Price (2014) also developed a physically based model to be used as a part of the Community Ice Sheet Model (CISM). This model is rather detailed in combining cavity formation (providing water storage) and a method to form Röthlisberger channelized flow. Analysis of this model looked at fine (100 m×100 m) grids and for shorter periods of time (on the order of days). While they do not mention the speed of their model, given that they model the growth and decay of channels, it is unlikely that the model would be suitable for longer timescales as the time stepping must be small to capture the transient nature of channelized flow. Likewise, for larger grid sizes, the effects of distributed systems (cavities, thin films, etc.) can be averaged out, saving on computation with minimal loss of generality in the results.

A distinguishing feature of this model is the numerical time stepping scheme. The model uses a combination of Heun's method and the leapfrog–trapezoidal schemes, which are iterative predictor–corrector schemes. The latter scheme has been used in the more demanding case of ocean modelling (Shchepetkin and McWilliams, 2005). The combination of these two methods (see Appendix A2) proves to be robust and stable with quick convergence to the final solution.

The hydrological model has been incorporated into the glacial systems model (Tarasov and Peltier, 1999; Tarasov et al., 2012). Below, we further detail and validate the subglacial model. We document water pressure and thickness sensitivity to hydrological parameters. Example results for the past North American ice complex are presented. Conclusions are summarized in Sect. 7.

2 Subglacial drainage systems

Subglacial drainage systems can be characterized as belonging to one of two categories: distributed drainage systems or channelized drainage systems.

2.1 Distributed drainage systems

There are several ways that water can be distributed underneath the ice: it can be stored via a thin film (Weertman, 1972) between the bed and the ice; it can be stored on the lee side of bed protrusions to form a linked-cavity system (Kamb et al., 1985), where braided canals (Clark and Walder, 1994) are formed as water cuts into underlying sediment; and it can flow through a porous medium via Darcian flow (Flowers, 2000). Distributed systems are inefficient at draining water. Therefore, these types of systems lead to a build-up of basal water pressure under the ice sheet.

2.2 Channelized drainage systems

The channelized drainage system is, to a certain degree, the obverse of the distributed drainage system. This system has a lot of water concentrated in a small area of the glacial bed and transports water quickly. Since channelized systems are efficient at draining water, they tend to decrease the water pressure, which increases basal friction between the ice and the bed. Thus, channelized systems are associated with slow flowing ice regimes and are often seasonal. Bartholomew et al. (2011) provides evidence that sliding velocities near the margins of the Greenland ice sheet are lower in the late summer than earlier in the summer, likely as an indication of a switch from a distributed to a channelized drainage system. There are two types of channelized drainage systems: Nye-channels that are incised down into the substrate (Walder and Hallet, 1979), and R-channels that are incised up into the ice (Röthlisberger, 1972).

3 Glacial systems model

For the analyses presented herein, the subglacial hydrology model is passively coupled to the glacial systems model (GSM). Full two-way coupling was turned off to isolate the dynamical response of the basal hydrology model. The GSM is composed of a thermo-mechanically coupled ice sheet model (using the shallow ice approximation), locally 1-D diffusive permafrost resolving bed thermal model (Tarasov and Peltier, 2007), fully coupled diagnostic surface drainage and lake storage module (Tarasov and Peltier, 2006), visco-elastic bedrock response, positive degree-day surface mass-balance with refreezing, and both marine and lacustrine calving parameterizations (Tarasov and Peltier, 1999, 2002, 2004; Tarasov et al., 2012).

The evolving temperature field (T) of the ice sheet is determined from conservation of energy:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}\left(T\right){\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left\{K_{\mathrm{T}}\left(T\right){\displaystyle \frac{\partial T}{\partial z}}\right\}-{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}\left(T\right)\mathit{u}\cdot \mathbf{\nabla }T+E_{\mathrm{d}},\end{array}$$

where c_i represents the specific heat of ice, ρ_i is the density of ice (910 kg m⁻³), K_T is the thermal conductivity of ice, u is the sliding velocity of the ice, and E_d is the heat created from the deformation of ice. As is standard, the horizontal diffusion component is ignored given the scales involved. The fully coupled ice and bed thermodynamics are solved via an implicit finite-volume discretization in the vertical and explicitly for the horizontal advection component of the ice thermodynamics. Basal temperature is limited to a maximum of the pressure melting point, with excess heat used to melt basal ice.

Basal sliding uses Weertman type sliding relations (i.e., function of driving stress) with power law 3 for hard bed and power law 1 for soft bed with sliding onset linearly ramped up starting from 0.2 ^∘C below the pressure melting point.

4 Subglacial hydrology model

4.1 Model description

For brevity and clarity, this section discusses the physical and numerical concepts of the hydrology model developed in Kavanagh (2012). The appendix provides details on the spatial discretization of the equations and the time stepping using the Heun/Leapfrog–trapezoidal scheme. This scheme is second order accurate. The model dynamically adjusts its internal time step to ensure the CFL criteria is satisfied (with the time step set to F_CFL × minimum CFL time step). Both of these features contribute to the stability of the model.

The dynamical evolution of distributed drainage is extracted from the mass continuity equation. Written in conservative form, the equation is

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial w}{\partial t}}+\mathrm{\nabla }\cdot \mathit{Q}=\dot{b}+d_{\mathrm{s}:\mathrm{a}},\end{array}$$

where w is the water thickness, and $\dot{b}$ is the meltwater source from the ice (negative if water refreezes to the ice). For the scope of this initial investigation, we assume no transmission of ice surface melt to the base. Observationally this is known to be false (Zwally et al., 2002), but the dependence on ice thickness, ice temperature profile, and ice strain profile makes this an issue deserving of its own focused study. The other source term, d_s : a, represents the drainage into the underlying aquifer.

The water flux, Q, is given by Darcy's law:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathit{Q}=-{\displaystyle \frac{Kw}{{\mathit{\rho }}_{\mathrm{w}}g}}\mathbf{\nabla }\left\{P+{\mathit{\rho }}_{\mathrm{w}}g{z}_{\mathrm{b}}\right\},\end{array}$$

where K is the hydraulic conductivity of the underlying till, ρ_w is the density of water (1000 kg m⁻³), g is the acceleration due to gravity (9.81 m s⁻²), z_b is the topographical bed elevation, and P is the water pressure beneath the ice.

The distributed flow of water beneath the ice sheet can come in many forms (such as cavities, Nye-channels, thin film, and flow through porous sediment). However, the extent to which the details of these flow mechanisms matter under large spatial scale separation and mechanistic heterogeneity is unclear. Therefore, we use the large difference in scale between glacial cycle ice sheet model grid cells (𝒪(10 km))) and distributed subglacial flow structures (𝒪(10 m) or less) to justify the choice of the diffusive Darcy flow equation for BrAHMs.

We use an empirical relation for water pressure from Flowers (2000, p. 68) given by

$$\begin{array}{}\text{(4)}& {\displaystyle }P=P_{\mathrm{I}}min\left[{\left({\displaystyle \frac{w}{h_{\mathrm{c}}}}\right)}^{\mathrm{7}/\mathrm{2}},\mathrm{1}\right],\end{array}$$

where P_I is the ice overburden pressure. P is limited to ice overburden pressure when w≥h_c. Saturated sediment water thickness, h_c¹, equals till thickness times porosity and is effectively the water thickness that the till can hold before becoming oversaturated (at which point the excess water is stored between the till and the ice). Flowers (2000) derived this equation by considering sub-grid variation in bed elevation and associated sediment thickness (and water thickness, all for the context of 40 m×40 m grid-cell modelling of Trapridge Glacier). A further consideration (without the overburden limit) was that the non-linearity would address dynamic adjustments in porosity and prevent stiffness in the dynamic equations caused by unrealistic heterogeneity in the modelled water distribution. Though derived for glacier-scale flow through a heterogenous macroporous sediment layer, our working hypothesis is that this empirical relation approximately captures large-scale pressure responses for subglacial distributed flow through any heterogenous structure (be it a mix of cavities of different size, patchy sediment, Nye-channels, and so on). The limiting of water pressure to overburden is justified by the low likelihood of water pressure exceeding the overburden pressure for any significant amount of time on glacial cycle time step scales.

The channelized system is likened to a system of R-channels (tunnels incised upward into the ice). Numerically, the model first calculates the water flux from the Darcian flow (Eq. 3). Channelized flow is invoked when that flux exceeds a critical value for the stability of the distributed regime given in Schoof (2010) as

$$\begin{array}{}\text{(5)}& {\displaystyle }\left|\mathit{Q}\right|<{\displaystyle \frac{\left|{\mathit{u}}_{\mathrm{b}}\right|Z_{\mathrm{h}}}{{\left({\mathit{\rho }}_{\mathrm{i}}L\right)}^{-\mathrm{1}}\left(\mathit{\alpha }-\mathrm{1}\right)\left|\mathbf{\nabla }\left(P+{\mathit{\rho }}_{\mathrm{w}}gz\right)\right|}},\end{array}$$

where u_b is the basal sliding velocity of ice, Z_h is the bedrock protrusion height, L is the latent heat of fusion of ice, and $\mathit{\alpha }=\mathrm{5}/\mathrm{4}$.

To simulate the change between different drainage systems, at regular user-defined intervals, the channel flow subroutine is called. Grid cells for which the water flux exceeds the distributed flow stability criterion (Eq. 5) are marked as tunnel grid cells. A down hydraulic gradient solver (Tarasov and Peltier, 2006) instantaneously² moves water down the path of steepest potential gradient (channelizing grid cells along that path) until there is no grid cell with a lower hydraulic potential (so forms subglacial lake) or the water exits the ice sheet. The solver considers all adjacent grid cells (including corner adjacency) when searching for the lowest hydraulic potential. Once the tunnel water transport is complete, the tunnels are assumed closed and the distributed flow algorithm continues.

4.2 Model coupling

BrAHMs is highly modular and designed for asynchronous coupling at user specified time steps. Aside from basic grid information, for each call, the hydrology model requires the following input fields: ice thickness, basal elevation, sea level, basal ice temperature, basal melt rate, and basal sliding velocity of the ice. For two-way coupling, the relevant outputs from BrAHMS are basal water pressure and thickness.

Given that there is no lower limit to coupling time steps, synchronous coupling can also be implemented. For two-way coupling, sensitivity tests are recommended to determine the appropriate coupling time step for the relevant context.

5 Model validation

5.1 Set-up

The basal hydrology model was subject to several validation tests with synthetic axisymmetric ice sheets. The three set-ups are: symmetrical ice sheet on a flat bed, symmetrical ice sheet on a dilating (sinusoidally wavy) bed, and a symmetrical ice sheet on an inclined plane.

5.1.1 Ice sheet profile

The continental-scale ice sheet model used in these tests has a profile that follows a normal distribution from the centre of the ice sheet to the terminus and is symmetric around the centre (i.e., bell-shaped ice sheet), according to the following equation:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }H\left(d\right)=\left(H_{\max }-H_{\min }\right)\exp \left[-{\left({\displaystyle \frac{d}{\sqrt{\mathrm{2}}{H}_{\mathrm{d}}r}}\right)}^{\mathrm{2}}\right]+H_{\min }\\ \text{(6)}& {\displaystyle }& {\displaystyle }\text{for}\left(d<r_{\mathrm{t}}\right).\end{array}$$

In Eq.6 H_max is the ice thickness at the ice divide (the centre), H_min is the thickness at the terminus, H_d is a normalized (by radius) standard deviation that defines how the profile spreads out, r_t is the distance to the terminus, and d is the distance from the ice divide, given by the following equation:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }d\left(\mathit{\theta },\mathit{\varphi }\right)={\cos }^{-\mathrm{1}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }\left[\sin \left({\mathit{\theta }}_{\mathrm{c}}\right)\sin \left(\mathit{\theta }\right)+\cos \left({\mathit{\theta }}_{\mathrm{c}}\right)\cos \left(\mathit{\theta }\right)\cos \left(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{c}}\right)\right]R_{\mathrm{Earth}},\end{array}$$

where θ and ϕ are the latitude and longitude, and R_Earth is the radius of the Earth (all parameter values are listed in Table 1).

Table 1 Model parameters used in validation studies.

5.1.2 Incline and dilating bed profiles

The bed for the inclined plane (as a function of latitude) is given by

$$\begin{array}{ll}{\displaystyle }{d}_{\mathit{\theta }}& {\displaystyle }={\cos }^{-\mathrm{1}}\left[\sin \left({\mathit{\theta }}_{\mathrm{N}}\right)\sin \left(\mathit{\theta }\right)+\cos \left({\mathit{\theta }}_{\mathrm{N}}\right)\cos \left(\mathit{\theta }\right)\right]R_{\mathrm{Earth}}\\ {\displaystyle }{d}_{\mathrm{S}}& {\displaystyle }={\cos }^{-\mathrm{1}}\left[\sin \left({\mathit{\theta }}_{\mathrm{N}}\right)\sin \left({\mathit{\theta }}_{\mathrm{c}}\right)+\cos \left({\mathit{\theta }}_{\mathrm{N}}\right)\cos \left({\mathit{\theta }}_{\mathrm{c}}\right)\right]R_{\mathrm{Earth}}\\ \text{(8)}& {\displaystyle }z& {\displaystyle }=Z_{\max }\left(\mathrm{1}-{\displaystyle \frac{d_{\mathit{\theta }}}{d_{\mathrm{S}}}}\right)\end{array}$$

to test the flow of water, where Z_max is the maximum elevation, d_S is the distance of the southern-most grid cell from the centre of the ice sheet (likewise, d_θ refers to any grid cell under consideration), and the slope is calculated such that z=0 at the south (d_θ=d_S).

The dilating bed is given by

$$\begin{array}{ll}{\displaystyle }z\left(\mathit{\theta },\mathit{\varphi }\right)=\min & {\displaystyle }\left[\left(Z_{\min }\right)\cos \left({\displaystyle \frac{\left(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{C}}\right)R_{\mathrm{Earth}}}{\mathrm{5}}}\right)\right.\\ \text{(9)}& {\displaystyle }& {\displaystyle }\left.\cos \left({\displaystyle \frac{\left(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{C}}\right)R_{\mathrm{Earth}}}{\mathrm{5}}}\right),\mathrm{0}\right],\end{array}$$

where Z_min is the maximum depth of the bed (all parameter values are listed in Table 1).

Figure 1 Simple synthetic ice sheet testing scenarios. Transects are plotted as the absolute, normalized distance from the centre. (a) A dome-shaped ice sheet placed on a flat bed resulting in symmetric results. (b) The ice dome on a dilating bed. Lake formation occurs where the ice sheet is relatively flat and the topography dips.

5.2 Model runs

It should be noted that the model is based on spherical polar coordinates (as it is designed for modelling large sections of the Earth's surface), and so the figures presented here are akin to the Mercator projection³ (mapping polar coordinates to Cartesian).

Tables 1 and 2 list the parameters used in the validation analysis. For Table 2, the values used for the validation tests are listed in the “Value” column.

In the model runs, the ice sheets starts from the ground (at t_now=0) and grows until 50 % of the model runtime (t_half). The ice thickness grows according to Eq. (6) multiplied by the ratio t_now / t_half. When t_now is greater than t_half, the ice sheet is at its maximum size (as shown in Fig. 1a).

To facilitate the growth of the subglacial hydraulic system, a constant melting at the base of the ice is applied in a “ring” of uniform thickness near the terminus, with 0.6 m yr⁻¹ of melting at the terminus that then linearly decreases to 0.4 m yr⁻¹ inside the “ring”. However, if there is no ice where the ring of meltwater is defined, then the value of melt, M_d, is set to zero until there is ice, in which case it would take the value defined by the equation

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{M}_{\mathrm{d}}\left(d\right)=M_{\mathrm{t}}-\left(M_{\mathrm{t}}-M_{\mathrm{i}}\right)\left({\displaystyle \frac{r_{\mathrm{t}}-d}{c_{\mathrm{r}}}}\right)\\ \text{(10)}& {\displaystyle }& {\displaystyle }\text{for}\left(r_{\mathrm{t}}-c_{\mathrm{r}}<d<r_{\mathrm{t}}\right),\end{array}$$

where M_t is the melt rate at the terminus, M_i is the melt rate inside the melt ring, and c_r is the thickness of the ring from the terminus into the innermost melting point.

Figure 2 Maximum mass balance discrepancy over time.

5.3 Validation results

5.3.1 Symmetric ice sheet on flat bed

The first set-up tested was for the ice sheet on a flat bed. For the bell-shaped ice sheet on a flat surface (see Fig. 1a), the model is mass conserving on the order of 10⁻¹² m of water thickness within a grid cell (Fig. 2). For this case, the water drains radially away from the ice sheet under the influence of the basal water pressure. The average water thickness was 0.0801±0.1357 m for the grid cells that contained water.

Johnson (2002)Person et al. (2012)Flowers et al. (2005)Flowers et al. (2005)Flowers et al. (2005)Flowers et al. (2005)Kamb (1987)

Table 2 Chosen values for the baseline model run for synthetic and North American test runs.

The transects in Fig. 1a show how the water thickness/pressure profiles change along two cross sections of the sheet through the centre (transects are plotted along the absolute-value, normalized distance from the centre of the ice sheet. For example, the water profile north of the centre is plotted on top of the south water profile for comparison.) For the case of the flat bed, both transects appear to be single lines, showing that there is a N–S symmetry and an E–W symmetry. Both transects are slightly offset from each other due to the grid discretization, but otherwise look similar.

Figure 3 Plots of the bell-shaped ice sheet on an inclined plane. This series also tests the convergence of the model. “Normal” refers to the model using the same resolution as the previous tests and “double” refers to the tests that has twice the resolution (finer) grid than the other tests.

5.3.2 Symmetric ice sheet on dilating bed

The next test was designed to show the formation of lakes by the model. The ice sheet for this study, as seen in Fig. 1b, is similar to the other ice sheets, but has its centre height (H_max) lowered from 3500 to 2000 m to smooth out the ice sheet toward the terminus; this was due to the fact that the previous profile was too steep, which led to no lake formations.

Figure 1b shows that (in areas where the ice is relatively flat and there is a dip in the bed) the hydrology model does allow for the build-up of water into subglacial lakes, reaching up to 100 m of water thickness in places.

The transect plots show that there is a slight asymmetry that arises in the results (there appears to be two red curves). Under perfect symmetries, the tunnel solver will break symmetry in its down-slope search algorithm. While the results are not shown here for brevity, when the tunnel solver is turned off, the results do not show any discernible asymmetry. The asymmetry due to the inclusion of the tunnel solver is unlikely to be an issue in more realistic cases where the ice sheet would lack such symmetry.

5.3.3 Symmetric ice sheet on an inclined plane

Figure 3 indicates that there is an asymmetry of N–S transects (in the interior of the ice around 0.8 normalized distance); this does not occur in the E–W transect, which still maintains its symmetry as in the flat bed case. This is due to the slope of the bed that results in a build-up of water in the northern section of the ice sheet. The build-up of water then reduces the hydraulic gradient in comparison to the hydraulic gradient in the southern section. The average water thickness in this scenario is 0.0802±0.1357 m.

Figure 3 shows the results for both the normal and double resolution grids. These two results are nearly identical, except the basal water is slightly thicker in the higher resolution plot (0.0805±0.1365 m), and are slightly offset from each other due to the changes in grid resolution. This suggests that the model is convergent at finer grid resolutions.

5.3.4 Model stability test

Lastly, the model was tested for stability by shocking the system with sudden changes in the meltwater production. For this test, the base case of the ice dome lying on a flat bed was used (the same scenario as Fig. 1a), only the meltwater production from Eq. (10) was modified by a time-dependent multiplier:

$$\begin{array}{ll}{\displaystyle }{M}_{\mathrm{d}}^{*}\left(d,t^{\prime }\right)& {\displaystyle }=\mathit{\xi }\left(t^{\prime }\right)\times M_{\mathrm{d}}\left(d\right)\\ {\displaystyle }\mathit{\xi }\left(t^{\prime }\right)& {\displaystyle }=\mathrm{1}+\sin \left[\mathrm{2}\mathit{\pi }\left({\displaystyle \frac{t^{\prime }-\mathrm{1}/\mathrm{8}}{\mathrm{80}}}\right)\right]\\ {\displaystyle }& {\displaystyle }\left(\mathrm{\Delta }\left[t^{\prime }-{\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}\right]-\mathrm{\Delta }\left[t^{\prime }-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right]\right)\\ {\displaystyle }& {\displaystyle }+\sin \left[\mathrm{2}\mathit{\pi }\left({\displaystyle \frac{t^{\prime }-\mathrm{3}/\mathrm{8}}{\mathrm{800}}}\right)\right]\left(\mathrm{\Delta }\left[t^{\prime }-{\displaystyle \frac{\mathrm{3}}{\mathrm{8}}}\right]-\mathrm{\Delta }\left[t^{\prime }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right]\right)\\ {\displaystyle }& {\displaystyle }+\sin \left[\mathrm{2}\mathit{\pi }\left({\displaystyle \frac{t^{\prime }-\mathrm{5}/\mathrm{8}}{\mathrm{80}}}\right)\right]\left(\mathrm{\Delta }\left[t^{\prime }-{\displaystyle \frac{\mathrm{5}}{\mathrm{8}}}\right]-\mathrm{\Delta }\left[t^{\prime }-{\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}\right]\right)\\ \text{(11)}& {\displaystyle }& {\displaystyle }+\sin \left[\mathrm{2}\mathit{\pi }\left({\displaystyle \frac{t^{\prime }-\mathrm{7}/\mathrm{8}}{\mathrm{800}}}\right)\right]\left(\mathrm{\Delta }\left[t^{\prime }-{\displaystyle \frac{\mathrm{7}}{\mathrm{8}}}\right]\right),\end{array}$$

where $t^{\prime }=t_{\mathrm{now}}/t_{\mathrm{final}}$ is the normalized time (t_final=800 years), and Δ(t^′) is the Heaviside step function. ξ(t^′) is plotted in Fig. 4a. Each impulse lasts for 1∕8th (100 years) of the model run. The lower frequency impulses contain 10 wavelengths each (period of 10 years), and the higher frequency impulses contain 100 complete wavelengths (period of 1 year). Each impulse will vary the meltwater production from a range of 0–2 times the base value.

Figure 4b shows how BrAHMs is affected by the varying meltwater production for changes in the CFL condition affecting time step size, and how it is affected by changes to the grid resolution. The analysis showed that there was no discernible difference from halving the CFL condition, so the results of changing the CFL condition were omitted. While there are some differences (along this particular transect) related to changes in the grid resolution, none of the plots exhibit any unrealistic behaviour due to the cycling of the meltwater production.

Figure 4 Time series analysis of basal water level response to time-varying meltwater production, showing that the dynamic time stepping of the model is capable of providing stable, realistic results under sudden changes.

6 Model results coupled to the GSM

The North American ice sheet model used herein is from a large ensemble Bayesian calibration as detailed in Tarasov et al. (2012). Model runs start from 122 ka under ice free conditions.

6.1 The model parameter set

Due to the complex nature of basal hydrology and the spatial and temporal scales for the current context, there are many processes that are approximated through parameterizations. As such, there are a number of poorly constrained parameters in the hydraulic model (these are listed in Table 2).

The first parameter, F_CFL, is used to control the time stepping of the model. As the model runs an explicit time scheme, it is subject to the CFL condition for stability. To help prevent the model from breaking the CFL condition, the model time step is dynamically altered to prevent the maximum basal water velocity from exceeding the CFL velocity. F_CFL determines the maximum allowable basal water velocity as a fraction of the CFL velocity. Should a time step exceed the CFL condition (potentially leading to instabilities), the last time step is redone with a smaller Δt such that the CFL condition is not broken.

The simplified aquifer drainage of Johnson (2002), uses an aquifer that simply drains a percentage of the present water in a grid cell. The percentage of water drained in this model is represented by the D_r parameter ($d_{\mathrm{s}:\mathrm{a}}=D_{\mathrm{r}}w$).

Due to the small time steps (relative to glacial modelling) involved in the basal hydrology model, it would become computationally expensive to check for tunnels at each time step. As such, dt_tun determines the frequency at which the model checks for the formation of channelized flow.

For clarity of this initial analysis, results presented herein are with a uniform basal sediment cover over the whole bed for the duration of the run. The sediment cover was specified by h_c.

The water flux between grid cells is directly proportional to the hydraulic conductivity of the sediment. For each run, the conductivity was allowed to vary between a minimum and maximum value defined in the range of K_m. The hydraulic conductivity is allowed to vary to account for the fact that the till will expand as it is filled with water, allowing for the water to flow with less resistance. The transition between low and high conductivity is controlled by the parameters k_a and k_b given by the equation (Flowers, 2000)

$$\begin{array}{ll}{\displaystyle }\log \left(K\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\left(\log \left[{\displaystyle \frac{K_{max}}{K_{min}}}\right]\right){\tan }^{-\mathrm{1}}\left[k_{\mathrm{a}}\left({\displaystyle \frac{w}{h_{\mathrm{c}}}}-k_{\mathrm{b}}\right)\right]\\ \text{(12)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\log \left[K_{max}{K}_{min}\right]\right),\end{array}$$

which is a constitutive equation of the logarithmic form of the upstream area of a grid cell (Flowers, 2000). Data from Trapridge Glacier shows a similar relation between upstream area and water pressure as the hydraulic conductivity equation (a high and low regime with a transition zone). Flowers (2000) assumes that the upstream area is related to the connectivity in the grid cell (the more connected a grid cell is, the more upstream area it should have). This would suggest that the hydraulic conductivity (its connectivity) is dependent on the water level, and is of the form of Eq. (12).

In Eq. (12), k_b affects (as a fraction of h_c) when the transition between begins (k_b∕h_c is the halfway point of the transition). k_a affects the slope of the transition curve. For larger values of k_a the transition becomes sharper, leading to quicker transitions. Lower values of k_a lead to slower transitions with more intermediate values for the conductivity between the two extremes (see Fig. A1). Equation (12) is meant to capture the dependence of hydraulic conductivity on the pore size of the till, which is related to the amount of water in the till.

As the base of the ice sheet becomes colder, the ice should begin to freeze to the bed, preventing water from flowing there. Due to the 40 km resolution of the grid, it is unlikely that the entire bed in a grid cell would be frozen completely when the grid cell basal temperature crosses the pressure melting point. Therefore, water could potentially flow through a frozen grid cell (in the unfrozen places), but the water should have a harder time as it has fewer pathways to flow across. In the hydrology model, this is represented by parameter T_c, which acts to reduce the conductivity as a function of temperature. When the basal temperature is close to the pressure melting point (PMP), there is little change in the hydraulic conductivity. Conductivity decreases to an extremal low value as the temperature approaches the value of T_c. In the model simulations, the value of T_c, relative to PMP, is tested from −0.5 to −3.0 ^∘C . As a simplifying assumption, the hydraulic conductivity of a frozen grid is set to 10⁻¹⁴ m s⁻¹, but can be easily modified to follow a temperature-dependent profile to capture sub-grid variation in basal temperatures.

Tunnel formation has a direct impact on basal water pressure. To further test this impact, an enhancement factor, Q_sc, was introduced to Eq. (5) as a multiplier to the condition for tunnel flow. Higher values of Q_sc will increase the switching condition, leading to less tunnel formation, whereas lower Q_sc will increase the amount of tunnel formation.

6.2 The baseline model

Our choice of baseline model for the sensitivity analysis was based solely on mid-range values for parameter uncertainty ranges and not on any sort of tuning. As such, results presented here have an exploratory instead of predictive focus.

Figure 5 Basal water profiles for (a) 22 ka when the total water volume is high (mean thickness: 1.0644±2.8317 m; max thickness: 86.40 m), and (b) 18 ka, after a large reduction in total basal water volume (mean thickness: 0.6918±1.3160 m; max thickness: 24.91 m). Five hundred metre contour intervals for surface elevation are also shown.

Basal hydrology fields for the baseline model near the last glacial maximum (LGM) are shown in Fig. 5. There is a greater extent and generally thicker basal water at 22 ka than at 18 ka. Regions of low basal effective pressure (defined as ice overburden pressure minus basal water pressure) in the model are generally associated with ice streaming. As the results are for a model configured with no basal drag dependence on basal water, this correlation is potentially due to two factors. First, high basal velocities will increase basal heating and subsequent basal meltwater production. Second, the basal flux trigger threshold for the initiation of tunnel drainage (Eq. 5) is proportional to basal velocity.

As the water is removed from 22 to 18 ka, some of the areas experience a large increase in basal effective pressure.

Figure 6 Sensitivity plot at (a) 22 ka and (b) 18 ka. Water storage for lowest D_r value is off the scale (221×10¹² m³ and 60×10¹² m³ for 22 and 18 ka, respectively).

To account for dependence on baseline amounts of basal water, our sensitivity tests consider both the 22 and 18 ka time slices. Figure 6a shows model sensitivity at 22 ka when the baseline model total water volume is higher. The most important parameter is the aquifer drainage parameter, which is the proportion of water drained locally out of the system. This simplified parameterization of the aquifer can quickly drain a lot of water as it does not have to flow to the terminus to escape and does not return it to the ice–bed interface. In Fig. 6b, at 18 ka, the aquifer drainage is still the most important parameter, but its impact is less noticeable since there is less water to drain away from the bed.

The sediment thickness parameter (h_c) shows a 28 % drop in water volume over the range of values at 22 ka. At 18 ka the impacts of h_c are greatly reduced and have no effect on water volume when raised above the baseline value. This is due to the non-linear relation between water pressure and sediment thickness from Eq. (4). In areas where the water level is only a small fraction of the sediment thickness, the basal water pressure will be practically zero. At 18 ka, when the water level is low, an increase in sediment will have little effect on the results.

The runs with the basal freezing value closer to the PMP have about a 12 (%) increase in basal water volume, as expected due to the increased likelihood of ice frozen to the bed hindering the flow of water. In comparison to other parameter results, varying the value of the basal freezing parameter, T_c, does not significantly alter the water storage. This is expected as regions where the basal temperature is below the PMP have no subglacial meltwater production.

The tunnel criterion scaling factor, Q_sc, shows almost no impact in times of high water storage, but shows a drop of up to 80 % in water volume at 18 ka. During this time, the model is sensitive to Q_sc because the lower water levels are less likely to form tunnels than the thicker values at 22 ka. Lowering Q_sc allows more tunnels to form, which drains the water, keeping the water volume down.

The bedrock bump height, Z_h, has a similar effect to Q_sc since it also affects tunnel formation. Larger values of Z_h allow the cavity system to retain more water before filling up and becoming unstable. This allows the runs with higher Z_h to have thicker basal water (Schoof, 2010).

The results of changing the range of hydraulic conductivity (K_m), show little difference at higher water volumes for the different runs. However, at 18 ka there is a big difference in the results. It is found that as hydraulic conductivity increases, the total water volume decreases. This is expected since increasing the conductivity increases the water flow and tunnel formation, allowing the water to evacuate from the ice sheet. The variation of k_a and k_b has little impact on the model results. This is rather fortuitous since they are not physical parameters that can be easily measured, whereas the range of hydraulic conductivity values can be constrained based on the type of sediment from field studies.

Figure 7 Sensitivity plot at (a) 22 ka and (b) 18 ka.

The plot of the average basal effective pressure, in Fig.  7, shares similar properties to the water storage sensitivity in Fig. 6. During periods of high water volume, the two most important parameters are the aquifer drainage and saturated sediment thickness. Their effects are much closer in terms of effective pressure due to the limiting of the pressure to ice overburden; thus, the effects of the aquifer drainage are limited.

During the low water storage times, the other parameters become important to the basal effective pressure. The impact of saturated sediment thickness on basal effective pressure appears to be relatively insensitive to the amount of water storage, as the two plots in Fig. 7 show similar results for both cases. Otherwise, the parameter values that lead to higher basal effective pressure are the same values that prevent the water from flowing out of the ice sheet in Fig. 6.

Figures 6 and 7 both show the lack of importance of the frequency of tunnel formation checks (dt_tun) which stems from the timescale of grid cell water refill (typically greater than the largest value for dt_tun). The effect of lowering dt_tun may have a minor effect on when the tunnels form, but not on how often.

One important test result is the low sensitivity of the average basal water thickness and effective pressure to the maximum allowable time step (dt_max). There is only a 10 % water volume drop in the range of values chosen. Also, as the time steps become smaller (the lowest value was 1∕3 of the baseline value), they begin to converge to an answer somewhere in the vicinity of the baseline values. This shows model stability and convergence for decreasing time steps.

As a caveat, these initial sensitivity tests likely hide spatially localized parametric sensitivities. More critically, feedbacks in a two-way coupled ice sheet and basal hydrology model configuration may strongly change relative sensitivities to basal hydrology parameters. These analyses will be better placed in a future study examining fully coupled dynamics.

7 Conclusions

This paper presents a physically based hydrology model for numerical simulations over a glacial cycle at continental scales. The model considers two types of drainage systems: a distributed system that slowly drains basal water, and a fast draining channelized system. The distributed hydrology system is modelled with Darcy's law (Flowers, 2000) while the channelized system is likened to R-channels and solved using a fast down-gradient routing and lake solver (Tarasov and Peltier, 2006).

The model was tested over a set of synthetic ice profiles and topography. The results of these tests show that the model is mass conserving and that the water flows down the hydraulic potential gradient where it can exit the ice sheet or form subglacial lakes.

With the model validated using the synthetic ice sheets, the model was then one-way coupled to the GSM for testing on the North American ice sheet complex at LGM. The sensitivity results in Figs. 6 and 7 show that the significance of each parameter varies in time as the amount of basal water in the system changes. In times of high water input, the only significantly influential parameters are sediment pore space and aquifer drainage parameters. During times of lower water levels, other parameters begin to impact the basal water thickness and pressure as well. These parameters are related to tunnel formation, such as the bedrock bump height, tunnel criterion scaling factor, and the hydraulic conductivity.

The hydrology model also identified areas of low effective pressure, indicating areas of potentially fast flowing ice. These results were self-consistent with the GSM's parameterized areas of the fast-flowing ice.

The hydrology model presented here has been shown to be stable and robust for the range of parameters used in this study. The coupled model generally takes 5–8 h to run for a North American glacial cycle (0.5^∘ longitude by 1.0^∘ latitude resolution). This time includes the full GSM, suggesting that the hydrology model only contributes an hour or two of extra runtime over a full glacial cycle. The longest runs are those with the smallest time steps (1∕120 year) or frequent calls to the tunnel solver, both of which show insignificant changes to the model results. This shows that the combined Heun's method and Leapfrog–trapezoidal scheme can be a viable numerical method for subglacial hydrology modelling.

As an initial implementation of a 2-D basal hydrology solver, there were several simplifications made to facilitate the initial study of the basic properties of the subglacial water dynamics. One simplification was that the aquifer drainage parameter was used instead of a real aquifer drainage system (Flowers, 2000; Lemieux et al., 2008), which would provide a more realistic drainage and allow water to flow back into the subglacial system. The sediment thickness was simplified as a constant over the entire bed. Realistically, the sediment thickness would vary over different parts of the bed (e.g., thinly covered Canadian Shield bedrock as opposed to the thick cover of the prairies), as well as varying in time as the sediment cover changes due to sediment deformation (Melanson, 2012).

Appendix A: Model numerics

A1 Discretization of the mass balance equation

The model uses the mass continuity equation (Eq. 2) for subglacial water. Expanding the divergence of the flux terms from the mass balance equation gives

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial w}{\partial t}}={\displaystyle \frac{\mathrm{1}}{r\cos \mathit{\theta }}}\left[{\displaystyle \frac{\partial \left(Q_{\mathit{\varphi }}\right)}{\partial \mathit{\varphi }}}+{\displaystyle \frac{\partial \left(Q_{\mathit{\theta }}\cos \mathit{\theta }\right)}{\partial \mathit{\theta }}}\right]+\dot{b}+d_{\mathrm{s}:\mathrm{a}},\end{array}$$

with θ representing the latitudinal direction and ϕ representing the longitudinal direction.

Equation (A1) is integrated over a finite-control volume

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\iint {\displaystyle \frac{\partial w}{\partial t}}\mathrm{dV}\\ \text{(A2)}& {\displaystyle }& {\displaystyle }=\iint \left\{{\displaystyle \frac{\mathrm{1}}{r\cos \mathit{\theta }}}\left[{\displaystyle \frac{\partial \left(Q_{\mathit{\varphi }}\right)}{\partial \mathit{\varphi }}}+{\displaystyle \frac{\partial \left(Q_{\mathit{\theta }}\cos \mathit{\theta }\right)}{\partial \mathit{\theta }}}\right]+\dot{b}+d_{\mathrm{s}:\mathrm{a}}\right\}\mathrm{dV}\end{array}$$

using dV=r²cosθdϕdθ, Eq. (A2) becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\iint {\displaystyle \frac{\partial w}{\partial t}}\mathrm{dV}=\underset{\mathrm{n}}{\overset{\mathrm{s}}{\int }}\left\{\underset{\mathrm{e}}{\overset{\mathrm{w}}{\int }}{\displaystyle \frac{\partial \left(Q_{\mathit{\varphi }}\right)}{\partial \mathit{\varphi }}}d\mathit{\varphi }\right\}rd\mathit{\theta }\\ \text{(A3)}& {\displaystyle }& {\displaystyle }+\underset{\mathrm{e}}{\overset{\mathrm{w}}{\int }}\left\{\underset{\mathrm{n}}{\overset{\mathrm{s}}{\int }}{\displaystyle \frac{\partial \left(Q_{\mathit{\theta }}\cos \mathit{\theta }\right)}{\partial \mathit{\theta }}}d\mathit{\theta }\right\}rd\mathit{\varphi }+\iint \left\{\dot{b}+d_{\mathrm{s}:\mathrm{a}}\right\}\mathrm{dV}.\end{array}$$

Figure A1 Variations of subglacial hydraulic conductivity, K, with respect to changes in hydraulic parameters k_a and k_b for values of K ranging between $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{7}}$ and $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$.

This then simplifies to

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\iint {\displaystyle \frac{\partial w_{\mathrm{P}}}{\partial t}}d{V}_{\mathrm{P}}=\underset{\mathrm{n}}{\overset{\mathrm{s}}{\int }}\left\{Q_{\mathrm{w}}-Q_{\mathrm{e}}\right\}rd\mathit{\theta }\\ {\displaystyle }& {\displaystyle }+\underset{\mathrm{e}}{\overset{\mathrm{w}}{\int }}\left\{Q_{\mathrm{s}}\cos {\mathit{\theta }}_{\mathrm{s}}-Q_{\mathrm{n}}\cos {\mathit{\theta }}_{\mathrm{n}}\right\}rd\mathit{\varphi }\\ \text{(A4)}& {\displaystyle }& {\displaystyle }+\iint \left\{\dot{b}+d_{\mathrm{s}:\mathrm{a}}\right\}d{V}_{\mathrm{P}},\end{array}$$

where the subscripts n, e, s, w stand for north, east, south, and west interfaces, respectively, and P represents the central grid point.

Using the approximation V_P=r²cosθ_PΔϕΔθ, Eq. (A4) can be approximated as follows:

$$\begin{array}{ll}{\displaystyle \frac{\partial w_{\mathrm{P}}}{\partial t}}=& {\displaystyle \frac{\mathrm{1}}{r\cos {\mathit{\theta }}_{\mathrm{P}}\mathrm{\Delta }\mathit{\theta }}}\left\{Q_{\mathrm{s}}\cos {\mathit{\theta }}_{\mathrm{s}}-Q_{\mathrm{n}}\cos {\mathit{\theta }}_{\mathrm{n}}\right\}\\ \text{(A5)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{r\cos {\mathit{\theta }}_{\mathrm{P}}\mathrm{\Delta }\mathit{\varphi }}}\left\{Q_{\mathrm{w}}-Q_{\mathrm{e}}\right\}+\dot{b}+d_{\mathrm{s}:\mathrm{a}}.\end{array}$$

A2 Model time stepping

A2.1 Predictor time steps

The model presented in this paper uses two predictor–corrector methods with different predictors but identical corrector. The first predictor method, based on Heun's method (Mathews and Fink, 2004), is used when only the current water thickness field is self-consistently available. It generates the first time step during each call of the basal hydrology subroutine, and for any grid cell that has been activated as a tunnel just prior to the current Darcy flow time step. The first of these two conditions could be easily amended to just the first time step of the whole model run if desired (and would be required for the case of synchronous model coupling). We chose the current formulation to simplify module coupling with configurations that involved multiple ice sheets on separate grids.

The first step in Heun's method is to take some initial conditions ($w_{\mathrm{P}}^m$), and do an Euler forward time step:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{w}_{\mathrm{P}}^{(m+\mathrm{1}{)}^{*}}=w_{\mathrm{P}}^m\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{\Delta }t}{r\cos {\mathit{\theta }}_{\mathrm{P}}}}\left\{{\displaystyle \frac{Q_{\mathrm{s}}^m\cos {\mathit{\theta }}_{\mathrm{s}}-Q_{\mathrm{n}}^m\cos {\mathit{\theta }}_{\mathrm{n}}}{\mathrm{\Delta }\mathit{\theta }}}+{\displaystyle \frac{Q_{\mathrm{w}}^m-Q_{\mathrm{e}}^m}{\mathrm{\Delta }\mathit{\varphi }}}\right\}\\ \text{(A6)}& {\displaystyle }& {\displaystyle }+\left({\dot{b}}^{\mathrm{0}}+d_{\mathrm{s}:\mathrm{a}}^{\mathrm{0}}\right)\mathrm{\Delta }t,\end{array}$$

where $w_{\mathrm{P}}^{(m+\mathrm{1}{)}^{*}}$ is the tentative (predicted) value for the next time step. The source terms (${\dot{b}}^{\mathrm{0}}+d_{\mathrm{s}:\mathrm{a}}^{\mathrm{0}}$) do not change within a call to BrAHMs and retain the 0 time index (a positivity constraint ensures that w_P≥0).

When the previous time step value of w_P is self-consistently available, we use the second order accurate Leapfrog–Trapezoidal scheme (the Heun scheme is only first order accurate). The Leapfrog predictor for this scheme calculates the intermediate predicted value for the next time step, $w^{(m+\mathrm{1}{)}^{*}}$, as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{w}_{\mathrm{P}}^{(m+\mathrm{1}{)}^{*}}=w_{\mathrm{P}}^{m-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{2}\mathrm{\Delta }t}{r\cos {\mathit{\theta }}_{\mathrm{P}}}}\left\{{\displaystyle \frac{Q_{\mathrm{s}}^m\cos {\mathit{\theta }}_{\mathrm{s}}-Q_{\mathrm{n}}^m\cos {\mathit{\theta }}_{\mathrm{n}}}{\mathrm{\Delta }\mathit{\theta }}}+{\displaystyle \frac{Q_{\mathrm{w}}^m-Q_{\mathrm{e}}^m}{\mathrm{\Delta }\mathit{\varphi }}}\right\}\\ \text{(A7)}& {\displaystyle }& {\displaystyle }+\mathrm{2}\left({\dot{b}}^{\mathrm{0}}+d_{\mathrm{s}:\mathrm{a}}^{\mathrm{0}}\right)\mathrm{\Delta }t.\end{array}$$

A2.2 Trapezoidal corrector

Regardless of which predictor equation is active, the trapezoidal scheme is applied to give the corrected value, w^m+1, as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{w}_{\mathrm{P}}^{m+\mathrm{1}}=w_{\mathrm{P}}^m+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}r\cos {\mathit{\theta }}_{\mathrm{P}}}}\\ {\displaystyle }& {\displaystyle }\left\{{\displaystyle \frac{Q_{\mathrm{w}}^{(m+\mathrm{1}{)}^{*}}-Q_{\mathrm{e}}^{(m+\mathrm{1}{)}^{*}}}{\mathrm{\Delta }\mathit{\varphi }}}+{\displaystyle \frac{Q_{\mathrm{s}}^{(m+\mathrm{1}{)}^{*}}\cos {\mathit{\theta }}_{\mathrm{s}}-Q_{\mathrm{n}}^{(m+\mathrm{1}{)}^{*}}\cos {\mathit{\theta }}_{\mathrm{n}}}{\mathrm{\Delta }\mathit{\theta }}}\right.\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{Q_{\mathrm{w}}^m-Q_{\mathrm{e}}^m}{\mathrm{\Delta }\mathit{\varphi }}}+{\displaystyle \frac{Q_{\mathrm{s}}^m\cos {\mathit{\theta }}_{\mathrm{s}}-Q_{\mathrm{n}}^m\cos {\mathit{\theta }}_{\mathrm{n}}}{\mathrm{\Delta }\mathit{\theta }}}\right\}\\ \text{(A8)}& {\displaystyle }& {\displaystyle }+\left({\dot{b}}_{\mathrm{s}}^{\mathrm{0}}+d_{\mathrm{s}:\mathrm{a}}^{\mathrm{0}}\right)\mathrm{\Delta }t.\end{array}$$

A3 Discretization of the Darcian flux

The Darcian flux, Q, is given in Eq. (3). The values for hydraulic conductivity, basal water thickness, and pressure, along with bed topography are calculated at the grid cell centres. To calculate Q at the grid cell interfaces, the aforementioned values must be assigned values at the interfaces.

If we consider the case of the flux on the westward edge of the grid cell, Q_w, then the pressure gradient (∇{P+ρ_wgz_b} from Eq. 3) is simply the difference between the pressure values at the grid cell centres

$$\begin{array}{}\text{(A9)}& {\displaystyle }{Q}_{\mathrm{w}}={\displaystyle \frac{Kw}{{\mathit{\rho }}_{\mathrm{w}}g}}{\displaystyle \frac{P_{\mathrm{W}}-P_{\mathrm{P}}+{\mathit{\rho }}_{\mathrm{w}}g\left(z_{{\mathrm{b}}_{\mathrm{W}}}-z_{{\mathrm{b}}_{\mathrm{P}}}\right)}{r\cos \left(\mathit{\theta }{-}_{\mathrm{P}}\right)\mathrm{\Delta }\mathit{\varphi }}},\end{array}$$

where the W subscript indicates the value of the grid point to the west of the central point, and the P subscript represents the grid cell of interest.

Following the rules of Patankar (1980), the hydraulic conductivity at the grid cell interface is set to the geometric mean of the values at the adjacent grid cell centres:

$$\begin{array}{}\text{(A10)}& {\displaystyle }{Q}_{\mathrm{w}}=\left({\displaystyle \frac{\mathrm{2}{K}_{\mathrm{W}}{K}_{\mathrm{P}}}{K_{\mathrm{W}}+K_{\mathrm{P}}}}\right){\displaystyle \frac{w}{{\mathit{\rho }}_{\mathrm{w}}g}}{\displaystyle \frac{P_{\mathrm{W}}-P_{\mathrm{P}}+{\mathit{\rho }}_{\mathrm{w}}g\left(z_{{\mathrm{b}}_{\mathrm{W}}}-z_{{\mathrm{b}}_{\mathrm{P}}}\right)}{r\cos \left({\mathit{\theta }}_{\mathrm{P}}\right)\mathrm{\Delta }\mathit{\varphi }}}.\end{array}$$

To simplify the flux equation, the upwind scheme (Patankar, 1980) was used to give the value of the water at the interface (i.e., the value of w is equal to the water thickness of the grid cell with the highest pressure). This gives the final equation of the flux as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{Q}_{\mathrm{w}}=\left({\displaystyle \frac{\mathrm{2}{K}_{\mathrm{W}}{K}_{\mathrm{P}}}{K_{\mathrm{W}}+K_{\mathrm{P}}}}\right)\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}gr\cos \left({\mathit{\theta }}_{\mathrm{P}}\right)\mathrm{\Delta }\mathit{\varphi }}}\right)\\ {\displaystyle }& {\displaystyle }[max\left\{w_{\mathrm{W}}\left[P_{\mathrm{W}}-P_{\mathrm{P}}+{\mathit{\rho }}_{\mathrm{w}}g\left(z_{{\mathrm{b}}_{\mathrm{W}}}-z_{{\mathrm{b}}_{\mathrm{P}}}\right)\right],\mathrm{0}\right\}\\ \text{(A11)}& {\displaystyle }& {\displaystyle }-max\left\{-w_{\mathrm{P}}\left[P_{\mathrm{W}}-P_{\mathrm{P}}+{\mathit{\rho }}_{\mathrm{w}}g\left(z_{{\mathrm{b}}_{\mathrm{W}}}-z_{{\mathrm{b}}_{\mathrm{P}}}\right)\right],\mathrm{0}\right\}],\end{array}$$

where Q_w is positively defined if water flows eastward into the centre grid cell. Likewise all the other fluxes can be defined in a similar fashion. Outgoing fluxes are limited to ensure positive basal water thickness.

A4 Flowchart of model procedure

Figure A2 Hydrology model flow chart highlighting the processes involved in simulating basal water flow.