the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### Mark Kavanagh

### Lev Tarasov

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.

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.

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).

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:

where *c*_{i} represents the specific heat of ice, *ρ*_{i}
is the density of ice (910 kg m^{−3}), *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.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

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:

where *K* is the hydraulic conductivity of the underlying till,
*ρ*_{w} is the density of water (1000 kg m^{−3}), *g* is the
acceleration due to gravity (9.81 m s^{−2}), *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

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}^{1}, 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

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^{2} 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.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:

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:

where *θ* and *ϕ* are the latitude and longitude, and *R*_{Earth}
is the radius of the Earth (all parameter values are listed in
Table 1).

### 5.1.2 Incline and dilating bed profiles

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

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

where *Z*_{min} is the maximum depth of the bed (all parameter values
are listed in Table 1).

## 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^{3} (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^{−1} of melting at the terminus that
then linearly decreases to 0.4 m yr^{−1} 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

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.

## 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^{−12} 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.

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.

### 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:

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.

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, *d**t*_{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)

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^{−14} m s^{−1}, 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.

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.

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.

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 (*d**t*_{tun}) which
stems from the timescale of grid cell water refill (typically greater
than the largest value for *d**t*_{tun}). The effect of lowering
*d**t*_{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 (*d**t*_{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.

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).

Basal hydrology code with validation drivers is freely available at http://doi.org/10.5281/zenodo.1230046 (Kavanagh and Tarasov, 2018).

## 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

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

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

using dV=*r*^{2}cos*θ**d**ϕ**d**θ*, Eq. (A2)
becomes

This then simplifies to

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*^{2}cos*θ*_{P}Δ*ϕ*Δ*θ*,
Eq. (A4) can be approximated as follows:

## 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:

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}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\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

### 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

## 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*+*ρ*_{w}*g**z*_{b}} from Eq. 3) is simply the
difference between the pressure values at the grid cell centres

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:

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

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

JH, BN, and TB performed the majority of the source code reconfiguration. BN provided overall oversight. JH designed the experiments, and carried them out with assistance from YM and DG. JH and BN prepared the manuscript with contributions from all co-authors.

The authors declare that they have no conflict of interest.

We thank Jesse Johnson and Julien Seguinot for their thoughtful and helpful
reviews. This work also significantly benefited from a review of the
associated Thesis by Jesse Johnson. Entcho Demirov offered some helpful
numerical suggestions. This work was supported by a NSERC Discovery Grant
(Lev Tarasov), the Canadian Foundation for Innovation (Lev Tarasov), and the
Atlantic Computational Excellence Network (ACEnet).

Edited by: Philippe Huybrechts

Reviewed by: Jesse Johnson and Julien Seguinot

Arnold, N. and Sharp, M.: Flow variability in the Scandinavian ice sheet: modelling the coupling between ice sheet flow and hydrology, Quaternary Sci. Rev., 21, 485–502, 2002. a, b, c

Bartholomaus, T. C., Anderson, R. S., and Anderson, S. P.: Growth and collapse of the distributed subglacial hydrologic system of Kennicott Glacier, Alaska, USA, and its effects on basal motion, J. Glaciol., 57, 985–1002, 2011. a

Bartholomew, I. D., Nienow, P., Sole, A., Mair, D., Cowton, T., King, M. A., and Palmer, S.: Seasonal variations in Greenland Ice Sheet motion: Inland extent and behaviour at higher elevations, Earth Planet. Sc. Lett., 307, 271–278, 2011. a

Benn, D. I. and Evans, D. J. A.: Glaciers and Glaciation, 2nd edn., Hodder Education, London, 2010. a, b

Bueler, E. and van Pelt, W.: Mass-conserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635, https://doi.org/10.5194/gmd-8-1613-2015, 2015. a, b, c, d

Calov, R., Beyer, S., Greve, R., Beckmann, J., Willeit, M., Kleiner, T., Rückamp, M., Humbert, A., and Ganopolski, A.: Simulation of the future sea level contribution of Greenland with a new glacial system model, The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-23, in review, 2018. a

Carter, S. P., Fricker, H. A., and Siegfried, M. R.: Antarctic subglacial lakes drain through sediment-floored canals: theory and model testing on real and idealized domains, The Cryosphere, 11, 381–405, https://doi.org/10.5194/tc-11-381-2017, 2017. a

Clark, P. U. and Walder, J. S.: Subglacial drainage, eskers, and deforming beds beneath the Laurentide and Eurasian ice sheets, Geol. Soc. Am. Bull., 106, 304–314, 1994. a

de Fleurian, B., Morlighem, M., Seroussi, H., Rignot, E., van den Broeke, M. R., Munneke, P. K., Mouginot, J., Smeets, P. C. J. P., and Tedstone, A. J.: A modeling study of the effect of variability on the effective pressure beneath Russell Glacier, West Greenland, J. Geophys. Res.-Earth, 121, 1834–1848, https://doi.org/10.1002/2016JF003842, 2016. a

Dow, C. F., Kulessa, B., Rutt, I. C., Tsai, V. C., Pimentel, S., Doyle, S. H., van As, D., and Lindbäck, K.: Modeling of subglacial hydrological development following rapid supraglacial lake drainage, J. Geophys. Res.-Earth, 120, 1127–1147, 2015. a

Flowers, G. E.: A multicomponent coupled model of glacier hydrology, PhD thesis, The University of British Columbia, BC, Canada, 2000. a, b, c, d, e, f, g, h, i, j, k, l

Flowers, G. E.: Subglacial modulation of the hydrograph from glacierized basins, Hydrol. Process., 22, 3903–3918, 2008. a, b

Flowers, G. E., Marshall, S., Björnsson, H., and Clarke, G. K. C.: Sensitivity of Vatnajökull ice cap hydrology and dynamics to climate warming over the next 2 centuries, J. Geophys. Res., 110, F02011, https://doi.org/10.1029/2004JF000200, 2005. a, b, c, d

Goeller, S., Thoma, M., Grosfeld, K., and Miller, H.: A balanced water layer concept for subglacial hydrology in large-scale ice sheet models, The Cryosphere, 7, 1095–1106, https://doi.org/10.5194/tc-7-1095-2013, 2013. a, b

Gudlaugsson, E., Humbert, A., Andreassen, K., Clason, C. C., Kleiner, T., and Beyer, S.: Eurasian ice-sheet dynamics and sensitivity to subglacial hydrology, J. Glaciol., 63, 556–564, 2017. a, b

Hoffman, M. J. and Price, S.: Feedbacks between coupled subglacial hydrology and glacier dynamics, J. Geophys. Res.-Earth, 119, 414–436, 2014. a

Johnson, J. V.: A Basal Water Model for Ice Sheets, Doctorate Thesis Abstract, PhD thesis, The University of Maine, USA, 2002. a, b, c, d, e, f

Kamb, B.: Glacier surge mechanism based on linked-cavity configuration of basal water conduit system, J. Geophys. Res., 92, 9083–9100, 1987. a, b

Kamb, B., Raymond, C. F., Harrison, W. D., engelhardt, H., Echelmeyer, K. A., Humphrey, N., Brugman, M. M., and Pfeffer, T.: Glacier Surge Mechanism: 1982–1983 Surge of Variegated Glacier, Alaska, Science, 227, 469–479, 1985. a

Kavanagh, M.: A fast, physically-based subglacial hydrology model applied to the North American ice complex over the last glacial cycle, Master's thesis, Memorial University of Newfoundland, Canada, 2012. a

Kavanagh, M. and Tarasov, L.: BrAHMs V1.0: A fast, physically-based subglacial hydrology model for continental-scale application, Zenodo, https://doi.org/10.5281/zenodo.1230046, last access: 24 April 2018. a

Lemieux, J. M., Sudicky, E. A., Peltier, W. R., and Tarasov, L.: Dynamics of groundwater recharge and seepage over the Canadian landscape during the Wisconsinian glaciation, J. Geophys. Res., 113, F01011, https://doi.org/10.1029/2007JF000838, 2008. a

Mathews, J. H. and Fink, K. K.: Numerical Methods Using Matlab, chap. 9, 4th edn., Prentice-Hall Inc., 2004. a

Melanson, A.: Numerical modelling of subglacial erosion and sediment transport and its application to the North American ice sheets over the last glacial cycle, unpublished masters thesis, Memorial University of Newfoundland, St. John's, Canada, 2012. a, b

Patankar, S. V.: Numerical Heat Transfer and Fluid Flow, Series in computational methods in mechanics and thermal sciences, Hemisphere Publishing Corporation, Washington, 1980. a, b, c

Person, M., Bense, V., Cohen, D., and Banerjee, A.: Models of ice-sheet hydrogeologic interactions: a review, Geofluids, 12, 58–78, 2012. a

Röthlisberger, H.: Water Pressure in Intra- and Subglacial Channels, J. Glaciology, 11, 177–203, 1972. a

Schoof, C.: Ice-sheet acceleration driven by melt supply variability, Nature, 468, 803–806, 2010. a, b, c, d

Shchepetkin, A. F. and McWilliams, J. C.: The regional oceanic modeling system (ROMS):a split-explicit, free-surface, topography-following-coordinate oceanic model, Ocean Model., 9, 347–404, 2005. a

Tarasov, L. and Peltier, W. R.: Impact of Thermomechanical ice sheet coupleing on a model of the 100 kyr ice age cycle, J. Geophys. Res., 104, 9517–9545, 1999. a, b

Tarasov, L. and Peltier, W. R.: Greenland glacial history and local geodynamic consequences, Geophys. J. Int., 150, 198–229, 2002. a

Tarasov, L. and Peltier, W. R.: A Geophysically constrained large ensemble analysis of the deglacial history of the North Americal ice-sheet complex, Quaternary Sci. Rev., 23, 359–388, 2004. a

Tarasov, L. and Peltier, W. R.: A calibrated deglacial drainage chronology for the North American continent: evidence of an Artic trigger for the Younger Dryas, Quaternary Sci. Rev., 25, 659–688, 2006. a, b, c

Tarasov, L. and Peltier, W. R.: The co-evolution of continental ice cover and permafrost extent over the last glacial-interglacial cycle in North America, J. Geophys. Res., 112, F02S08, https://doi.org/10.1029/2006JF000661, 2007. a

Tarasov, L., Dyke, A. S., Neal, R. M., and Peltier, W. R.: A data-calibrated distribution of deglacial chronologies for the North American ice complex from glaciological modeling, Earth Planet. Sc. Lett., 315–316, 30–40, 2012. a, b

Walder, J. S. and Hallet, B.: Geometry of former subglacial water channels and cavities, J. Glaciol., 23, 335–346, 1979. a

Weertman, J.: General theory of water flow at the base of a glacier or ice sheet, Rev. Geophys., 10, 287–333, 1972. a

Werder, M. A., Hewitt, I. J., Schoof, C. G., and Flowers, G. E.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res.-Earth, 118, 2140–2158, 2013. a

Zwally, H. J., Abdalati, W., Herring, T., Larson, K., Saba, J., and Steffen, K.: Surface Melt-Induced Acceleration of Greenland Ice-Sheet Flow, Science, 297, 218–222, https://doi.org/10.1126/science.1072708, 2002. a

^{1}

*h*_{c} can be thought of as the
critical water thickness before the basal water pressure reaches overburden
and quickly rises as more water is added.

^{2}

During the tunnel flow, no model time is stepped as tunnel drainage is computed diagnostically.

^{3}

Set
between 17 and 40^{∘} N for the simplicity of avoiding any issues such
as having negative values for generating the synthetic geometry.