the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
SUHMO: an adaptive mesh refinement SUbglacial Hydrology MOdel v1.0
Anne M. Felden
Daniel F. Martin
Esmond G. Ng
Water flowing under ice sheets and glaciers can have a strong influence on ice dynamics, particularly through pressure changes, suggesting that a comprehensive ice sheet model should include the effect of basal hydrology. Modeling subglacial hydrology remains a challenge, however, mainly due to the range of spatial and temporal scales involved – from subglacial channels to vast subglacial lakes. Additionally, networks of subglacial drainage channels dynamically evolve over time. To address some of these challenges, we have developed an adaptive mesh refinement (AMR) model based on the Chombo software framework. We extend the model proposed by Sommers et al. (2018) with a small but significant change to accommodate the transition from unresolved to resolved flow features. We handle the strong nonlinearities present in the equations by resorting to an efficient nonlinear full approximation scheme multigrid (FASMG) algorithm. We outline the details of the algorithm and present convergence analysis results demonstrating its good performance. Additionally, we present results validating our approach, using test cases from the Subglacial Hydrology Model Intercomparison Project (SHMIP) (de Fleurian et al., 2018). We finish by presenting a more complex, 100 kmby100 km synthetic test case with peaks and valleys that we use to investigate the effective pressure distribution as the number of AMR levels increases. These preliminary results suggest that a minimum spatial resolution is needed to properly capture channel features, but additional work is required to precisely quantify this and its impact on accurately modeling the coupled ice sheet–hydrology system. The efficiency of our approach, relying on localized refinement, is also demonstrated. Future work will include coupling the SUbglacial Hydrology MOdel (SUHMO) with the BISICLES AMR ice sheet model (Cornford et al., 2013), both built on the same numerical framework.
 Article
(7560 KB)  Fulltext XML

Supplement
(1580 KB)  BibTeX
 EndNote
The extensive and accelerating retreat of glaciers observed over the last 150 years has fueled interest in the behavior of the cryosphere. It is common knowledge that ice sheets have grown and retreated many times over the past 2.6 million years, but some studies suggest that the current interglacial period is straying from its expected course and could last much longer than originally anticipated (e.g., Berger and Loutre, 2002; Ganopolski et al., 2016). The amount of water stored in the Greenland Ice Sheet alone has the potential to raise the globalmean sea level (GMSL) by about 7 m (Aschwanden et al., 2019) at a rate being characterized by deep uncertainties to external factors (Edwards et al., 2021). The contribution of glaciers and ice sheets to GMSL has been increasing over the past decades, combining to more than 50 % of the total change over the period 2006–2018 (Edwards et al., 2021; MassonDelmotte et al., 2021). Recent studies predict a likely GMSL rise of about 0.6 m by 2100 under intermediate greenhouse gas emissions scenarios, but ice sheet processes could drive GMSL to rise up to about 2 m by 2100 and 5 m by 2150 under a very high emissions scenario (MassonDelmotte et al., 2021). Since more than 30 % of the world's population today lives in what can be considered coastal areas, predicting and understanding the mechanics governing the melting of the cryosphere can be considered a problem of pressing scientific and societal importance.
The evolution of the cryosphere can hardly be decoupled from that of its environment, and numerical ice sheet models should ideally include simplified but accurate representations of the interactions with the other components of the Earth system (see, e.g., Fig. 1 of Goelzer et al., 2017). In particular, it is now widely accepted – through observational evidence and theoretical considerations – that basal hydrology strongly influences the dynamics of glaciers and ice sheets, mainly through changes in the basal pressure (e.g., Iken, 1981; Bindschadler, 1983; Fricker et al., 2007; Stearns et al., 2008; Doyle et al., 2018). Concurrent ice uplift and acceleration have been observed in both Greenland and the Antarctic Ice Sheet (Das et al., 2008; Nienow et al., 2017; Tuckett et al., 2019), suggesting significant routing of surface lake water to the base and an active and dynamic subglacial hydrology network. Persistent subglacial water structures like conduits – from one draining event to the next – and lakes have also been postulated to exist in Antarctica, which could be of critical importance in the estimation of the ice sheet mass balance and, eventually, GMSL projections (Siegfried and Fricker, 2018; Kirkham et al., 2019; Malczyk et al., 2020).
Modeling subglacial hydrology remains a challenge, mainly due to the large discrepancy in spatial and temporal scales between the physics driving subglacial phenomena versus those driving ice sheet dynamics. Relevant timescales for the latter are on the order of tens to thousands of years, while subglacial water cycles occur over days or months. Englacial structures can organize into channels a few meters wide, while Antarctica covers an area of over 10 million km^{2}. Additionally, the variety of subglacial and englacial water structures, including sheets (Weertman, 1962) and cavities (Walder, 1986; Kamb, 1987), icewalled conduits (Röthlisberger, 1972), bedrock channels, lakes and everything in between, makes it difficult to construct a comprehensive subglacial hydrology model (see, for example, Fig. 2 of Flowers, 2015). Complicating matters further, existing observations are insufficient to fully validate any modeling effort and, to date, no complete, physically based theory has been developed. Most models follow more or less the same “recipe” made up from drainagesystem elements that are combined depending on what is thought to be physically relevant and on the numerical framework at hand. These models typically resort to a differentiation between inefficient flow configurations (linked cavities) and efficient ones (channels). Both are described by a balance between opening and closing terms: cavities typically open up due to sliding over the bedrock, while channels open up due to dissipative heating and melting of their walls, and closing in both configurations is due to ice creep. While a review of all subglacial modeling efforts is out of the scope of the present paper, we briefly discuss recent relevant twodimensional, multielement efforts which share the ultimate goal of coupling with an ice sheet model, to provide context to our study. The reader is referred to the thorough review of Flowers (2015) as well as to the references and participating models in the Subglacial Hydrology Model Intercomparison Project (SHMIP) (de Fleurian et al., 2018) for additional information.
“Nextgeneration” efforts (terminology from Flowers (2015)) started emerging in the early 2000s, with pioneering work to simulate the Weichselian Scandinavian Ice Sheet by Arnold and Sharp (2002). In their study, ice velocities relied on calculations of subglacial water pressures and the use of a waterpressuredependent sliding law. The twodimensional basal hydrology is made up of either inefficient (linked cavities) or efficient (Röthlisberger channels) flow configurations – both cannot coexist in the same cell, and findings highlighted that the spatial distribution of ice flow is greatly impacted by the presence of the subglacial water. More recently, mainly motivated by observations made in Greenland linking surface meltwater and basal sliding, the last decade has seen a sustained interest in comprehensive subglacial hydrology models. Schoof (2010) used a model with a single equation for the crosssectional area of discrete conduits that could behave as either cavities or channels to investigate the link between ice velocities and subglacial water channelization patterns generated by seasonal and shortterm water supply variations. Hewitt et al. (2012) and Hewitt (2013) introduced a twodimensional subglacial hydrology model coupling a continuum sheet and discrete channel elements – each requiring a different set of equations. This model is linked to an ice flow model, enabling a parametric study on a synthetic sheetlike geometry emulating Greenland Ice Sheet margins. Werder et al. (2013) extended this effort to an unstructured finiteelement mesh where the distributed continuum sheet is solved using finite elements on a set of triangular cells while channels are located along the edges of the cells. The resulting model, the Glacier Drainage System (GlaDS) model, is coupled to both the Icesheet and Sealevel System Model (ISSM) (Larour et al., 2012) and the Elmer/Ice ice sheet model (Gagliardini et al., 2013) and has been successfully employed to study both the Greenland and Antarctic ice sheets (e.g., Dow et al., 2016; Gagliardini and Werder, 2018). SHAKTI (Sommers et al., 2018), also coupled to ISSM, removes the distinction between channels and cavities/sheets by using a single set of equations to evolve the water gap height (unlike the model of Schoof (2010) that evolves a drainage crosssectional area). Such an approach is attractive because water structures are free to evolve and merge anywhere in the domain (without being restricted to cell faces, for example), without the need for an explicit coupling between them. It does, however, raise the question of how fine the mesh must be to properly resolve the various subglacial features. The cost of finer resolution to accommodate the formation of channels could very quickly become prohibitive. Fortunately, fine structures are expected to be localized both spatially and temporally, making this a promising target for adaptive mesh refinement (AMR).
We follow the approach of Sommers et al. (2018) with adaptations to implement a subglacial hydrology model based on the Chombo AMR framework (Adams et al., 2001–2021). Chombo provides a set of tools for implementing finitevolume methods for the solution of partial differential equations on blockstructured adaptively refined rectangular grids. We propose a small but significant modification to the set of equations presented in Sommers et al. (2018) in order to seamlessly transition from underresolved to resolved channels, alleviating an unphysical asymptotic behavior as the mesh size begins to allow resolution of typical channel scales. Building on the AMR framework allows development of a robust numerical approach to solve the resulting nonlinear system of PDEs which achieves secondorder convergence in space.
The paper is structured as follows. In Sect. 2, we summarize the set of equations used in the SUbglacial Hydrology MOdel (SUHMO). In Sect. 3 we provide the details of the nonlinear full approximation scheme (FAS) algorithm employed to solve the governing equations of Sect. 2. Convergence analyses demonstrating the efficiency and accuracy of our implementation are presented in Sect. 4. In Sect. 5, we present additional validating results, choosing three representative test cases extracted from SHMIP (de Fleurian et al., 2018). Results from a transient, largerscale AMR simulation with random bed roughness and interesting topographic features are discussed in the final Sect. 6. We finish with concluding remarks.
We start with a set of equations similar to that used in SHAKTI (Sommers et al., 2018), which is a parallelized, finiteelement subglacial hydrology model currently implemented as part of the opensource Icesheet and Sealevel System Model (Larour et al., 2012). We provide a brief overview before introducing a novel diffusion component that we believe represents an improvement to the existing model. For additional details concerning the derivation of the original equations, the reader is referred to Sommers et al. (2018) or to earlier work by, e.g., Hewitt (2011) or Hewitt et al. (2012).
2.1 Equations
The governing equation set starts with a twodimensional expression for the conservation of mass – assuming we are dealing with an incompressible fluid:
where b is the subglacial waterfilled gap height (m), b_{e} is the volume of water stored englacially per unit area of bed (m), q is the gapintegrated basal water flux (m^{2} s^{−1}), $\dot{m}$ is the melt rate (kg m^{−2} s^{−1}), ρ_{w} is the density of water (kg m^{−3}), and e_{s} encompasses all external sources of meltwater (produced englacially or surface meltwater, for example) (m s^{−1}).
An approximate momentum equation for water velocity integrated over the gap height gives rise to an expression for the water flux, based on equations developed for flow in rock fractures (e.g., Zimmerman et al., 2004):
where g is the gravitational acceleration (m s^{−2}), ν is the kinematic viscosity of water (m^{2} s^{−1}), ω is a dimensionless parameter controlling the nonlinear transition from laminar to turbulent flow and Re is the Reynolds number. The hydraulic head h (m) is defined as
where P_{w} is the water pressure (Pa) and z_{b} is the bed elevation (m). The Reynolds number follows a classical definition:
where v is the average flow velocity across the gap height. Equation (2) is an important piece in the SHAKTI model (Sommers et al., 2018). It allows for a spatially and temporally variable hydraulic transmissivity in the system and facilitates the representation of the simultaneous coexistence of laminar, transitional, and turbulent flow in subregions of the domain. The reader is referred to the aforementioned publication for additional details.
The melt rate $\dot{m}$ includes heat produced at the bed (geothermal flux and frictional heat due to sliding over the bed) along with heat generated through internal dissipation (mechanical energy converted to thermal energy by the flow), which is effectively melting the drainage system's walls and ceiling:
where L is the latent heat of fusion of water (J kg^{−1}), G is the geothermal flux (W m^{−2}), u_{b} is the ice basal velocity vector (m s^{−1}), τ_{b} is the stress exerted by the bed onto the ice (Pa), c_{t} is the change in pressuremelting point with temperature (K Pa^{−1}) and c_{w} is the heat capacity of water (J kg^{−1} K^{−1}). Note that the last term takes into account changes in sensible heat due to pressuremelting point variations. It can become important when considering realistic topographies with sharp changes in bed elevation (z_{b}) but is often neglected in similar models (de Fleurian et al., 2018). We also neglect this term in SUHMO, unless explicitly stated otherwise (such as for the AMR test case in Sect. 6).
Since we do not allow for the drainage space to be partially filled, b also evolves according to opening and closing terms that are typically modelspecific. As in Sommers et al. (2018), opening can be due to melt and sliding over bumps on the bed, while closing is solely due to ice creep:
where ρ_{i} is the ice density (kg m^{−3}), u_{b} is the magnitude of the sliding velocity (m s^{−1}), A is the ice flow parameter (Pa^{−3} s^{−1}), n is the flow law exponent (typically, n=3) and P_{i} is the ice overburden pressure (Pa). The parameter $\mathit{\beta}=max\left(\right({b}_{\mathrm{r}}b)/{l}_{\mathrm{r}},\mathrm{0})$ is dimensionless; it governs opening by sliding and is a function of the typical bed bump height (b_{r}) and bump spacing (l_{r}) in such a way that opening by sliding only occurs where the gap height is less than the typical local bump height (Werder et al., 2013). The quantity l_{c} is the creep length scale, which is defined as follows:
with b_{c} a critical gap height controlling the creep. Note that previous work, including SHAKTI (Sommers et al., 2018), typically sets l_{c}=b. By using Eq. (7) to define the ice creep length scale, we ensure that ice creep can be more rapidly cut off as soon as b≤b_{c}, thus enabling sheetlike systems to survive. This can be useful, for instance, when the bed topography is such that big lubricated areas are known to persist.
Equations (1) and (6) can then be combined to produce an equation for the evolution of the hydraulic head:
2.2 Model parameters
Constants and parameters presented in the previous section are summarized in Table 1, along with typical values. Note that the englacial storage volume, b_{e}, present in the original set of equations (Sommers et al., 2018) is not currently used in SUHMO. Under this assumption, Eq. (8) becomes a standard elliptic partial differential equation (PDE).
2.3 Introduction of a diffusion term
As described in the introduction, the subglacial drainage system is made of various coexisting and dynamically evolving subglacial drainage structures. These are typically broken into two categories: distributed (or inefficient) flow structures and channelized (or efficient) flow structures. The advantage of this set of governing equations is that only one equation (i.e., Eq. 6) is necessary to model the drainage space, and this equation accommodates both inefficient and efficient elements. One drawback from Eq. (6), however, is that it was originally derived by considering a sheet drainage system, central to which is the growth of the water sheet thickness, while equations governing channel growth are intrinsically twodimensional. In practice, adding a simple meltopening term in Eq. (6) to accommodate channels causes the water sheet thickness to grow locally (e.g., in each computational cell), because the second dimension, the channel width, has been dropped from the problem formulation.
This behavior is not an issue if relatively coarse meshes are used in configurations where the primary interest is effective pressure fields, but it does prevent channel geometries from converging with mesh resolution. Convergence with mesh resolution is an important feature of consistent numerical methods; additionally, AMR requires consistent convergence with mesh resolution to be effective. As a first step towards addressing this issue, we modify Eq. (6) by adding a diffusionlike term as follows:
where the diffusion coefficient depends on the dissipation term of the melt rate:
With this formulation, we aim to represent melting of channel walls as heat dissipation is no longer limited to channel/cavity ceilings. Adding this diffusion term and neglecting englacial storage, Eq. (8) then becomes
We solve Eqs. (9) and (11) on a hierarchy of blockstructured, Cartesian meshes using a finitevolume discretization, facilitated by the Chombo framework (Adams et al., 2001–2021). We extend the Chombo toolbox to solve the nonlinear evolution equation for the hydraulic head implicitly using the full approximation scheme (Briggs et al., 2000). The resulting algorithm accuracy is second order in space and first order in time. For completeness, Appendix A gives a brief summary of the main features of our AMR framework. We adopt the notation used in several previous studies (Berger and Colella, 1989; Martin et al., 2008; Cornford et al., 2013; Parkinson et al., 2020), and the reader is referred to these prior publications for additional details. All that follows is specific to a twodimensional application with an isotropic Cartesian mesh.
3.1 Full approximation scheme for variablecoefficient, nonlinear elliptic PDE
3.1.1 Multigrid methods
Multigrid (MG) methods are commonly employed to solve linear or nonlinear problems which can be expressed as
If we denote by $\stackrel{\mathrm{\u0303}}{u}$ an approximation to the exact solution u of this problem, we can define the error e as $e=u\stackrel{\mathrm{\u0303}}{u}$ and the residual r as $r=\mathcal{F}\mathcal{A}\left(\stackrel{\mathrm{\u0303}}{u}\right)$. If 𝒜 is a linear operator, then we obtain the residual equation:
with which we can start a recursive process of “restricting” (averaging) the residual onto a coarser grid to solve for a coarse correction of the error that we then interpolate back (usually referred to as “prolongation”) to the fine grid (Briggs et al., 2000). For linear systems, multigrid is highly efficient and can be implemented using a matrixfree approach. It is also straightforward to extend to AMR mesh hierarchies (Martin and Cartwright, 1996).
3.1.2 FAS multigrid
If 𝒜 is a nonlinear operator, we cannot make use of Eq. (13). Instead we use the residual in a new problem formulation:
with which we can follow a similar recursive process to a traditional MG method. For completeness, the main steps of one FAS MG iteration are described here and illustrated in Fig. 1, focusing on the Vcycle scheduling of events; a more thorough description can be found elsewhere (Briggs et al., 2000; Henson, 2003). We have chosen to work with h instead of u to be consistent with our main variable – the hydraulic head. We will use the following general expression for the nonlinear operator 𝒜:
where the nonlinear piece is contained in 𝒢(h). The restriction process is illustrated on the left of Fig. 1, by the arrow pointing down, while the prolongation is illustrated on the right with the arrow going up. The original level is labeled the finest level, and it has a mesh spacing of Δx_{0}. This level is coarsened n times by a factor of 2, so that coarsened grid n has a grid size of 2^{n}Δx_{0}.
We start by finding an approximate solution $\stackrel{\mathrm{\u0303}}{h}$ to the original problem Eq. (12) on the finest grid by using a few iterations of an iterative method (the relaxation method is described in Appendix B). The residual r is then evaluated, and both $\stackrel{\mathrm{\u0303}}{h}$ and r are averaged down to the first coarsened grid via a restriction operator denoted as ${R}_{\mathrm{\Delta}{x}_{\mathrm{0}}}^{\mathrm{2}\mathrm{\Delta}{x}_{\mathrm{0}}}$. We rewrite the problem as in Eq. (14) (where the new righthand side – RHS – is denoted ${\mathcal{F}}^{\mathrm{2}\mathrm{\Delta}{x}_{\mathrm{0}}}$ in Fig. 1) and perform a number of relaxation iterations to find an approximate solution, $\stackrel{\mathrm{\u0303}}{{h}^{\mathrm{2}\mathrm{\Delta}{x}_{\mathrm{0}}}}$. This process is repeated until we reach the coarsest level. The coarsest error, ${e}^{{\mathrm{2}}^{n}\mathrm{\Delta}{x}_{\mathrm{0}}}$, can then be extracted from the solution and interpolated up to the next finer grid to be used to correct the local approximation $\stackrel{\mathrm{\u0303}}{{h}^{{\mathrm{2}}^{n\mathrm{1}}\mathrm{\Delta}{x}_{\mathrm{0}}}}$ (see Fig. 1). This “prolongation” process is repeated, going up the levels until we reach the finest level again. We usually perform between two and four relaxation iterations while correcting the solution at each level.
3.1.3 Special treatment for the variable coefficient
We rewrite Eq. (11) to make it consistent with the operator that has been described in Fig. 1:
where ${\mathit{\rho}}_{\mathrm{c}}=(\frac{\mathrm{1}}{{\mathit{\rho}}_{\mathrm{w}}}\frac{\mathrm{1}}{{\mathit{\rho}}_{\mathrm{i}}})$. We discuss our treatment of the remaining nonlinearity and h dependencies in the RHS in Sect. 3.2. Note that the coefficient B is both spatially variable and a function of the primary variable h, via the coupling with Re. In our final implementation of the algorithm, B is recomputed on the finest grid at the beginning of every V cycle and then averaged down on all coarser grids. We experimented with fully lagging B, estimating the Reynolds number before the first Vcycle iteration and freezing its value for the remaining of the solve, but this resulted in poor overall algorithmic efficiency, as discussed in Appendix C.
3.2 Map of the algorithm
In order to solve the coupled set of Eqs. (9) and (16), we use a combination of FASMG iterations as described in Sect. 3.1 to solve for the hydraulic head and traditional MG to solve for the gap height. We note that the lagged term in Eq. (16), as well as any term involving the water flux, will depend on h in a nontrivial way. To ensure these are treated properly, we embed the FASMG solve for h in external Picard iterations (in which the value of the gap height is frozen). We have found that the required number of Picard iterations is, on average, typically one or two. This number can occasionally increase, particularly during transient phases.
We use a backward Euler method to discretize the temporal term in Eq. (9), reorganizing to be consistent with Eq. (12):
where we have replaced the effective pressure (P_{i}−P_{w}) with N, and n=3 has been assumed.
The main steps required to advance our main variables from t=t^{n} to $t={t}^{n+\mathrm{1}}$ are summarized in Algorithm 1. The superscripts n and k refer to the current time step and the current Picard iteration, respectively. For better readability, n is omitted when discussing Picard iterations. When relevant, “CC” will refer to cellcentered variables, while “FC” will refer to facecentered variables.
Following common usage (for example, in Martin and Colella, 2000), we enforce boundary conditions at domain boundaries and between AMR grid patches (both at the same refinement level and between levels) through the use of a ring of ghost cells (“GC”) around each logically rectangular patch.
4.1 Distributed test case
We start by analyzing a simple distributedflow test case. The topography considered is the synthetic representation of a landterminating ice sheet margin from SHMIP (de Fleurian et al., 2018), illustrated in Fig. 2a. The ice sheet domain measures 100 km in the x direction and 20 km in the y direction, the bed is flat, and a parabolic ice surface z_{s} is prescribed by
In order to evaluate the convergence properties of the complete algorithm, we evolve the system to a fixed time with increasing resolution, halving Δx with each refinement. The solution error is then computed by comparing Φ^{c}, the solution at resolution Δx, with the finer solution ϕ^{f} computed using $\mathrm{\Delta}x/\mathrm{2}$ and averaged onto the coarser grid (Φ^{f→c}). The L^{2} norm of the error for a simulation with n_{tot} cells is
Figure 3a shows the error convergence using six grid resolutions for two variables of interest: h and N. The slope of the error demonstrates secondorder convergence for both variables.
4.2 Channelized test case
We now turn our attention to a channelizing test case. The domain is a rectangle of 64 m in the x direction by 16 m in the y direction. The bed is sloped in the x direction (with a +2 % slope) and is topped with a slab of ice of constant thickness everywhere (500 m). A moulin delivering 30 m^{3} s^{−1} of water is located 16 m from the outlet of the domain in the x direction. The geometry is shown in Fig. 2c. Periodicity is assumed in the y direction, a homogeneous Dirichlet boundary condition (outlet) is prescribed at x=0 m, and a homogeneous Neumann boundary condition is prescribed at x=64 m. Note that the moulin source term follows a spatial Gaussian profile, as shown in the top plot of Fig. 4a – the convergence of which is exactly second order. Additionally, the moulin input is gradually increased in time, from 0 at time t=0 s to the maximum value after about a month, and the simulation proceeds until steady state is reached.
We use Eq. (19) to evaluate the convergence properties of SUHMO in this slightly more challenging test case. Results using seven grid resolutions are presented in Fig. 3b and c. The original formulation (Eqs. 6 and 8) fails to converge: as the resolution increases, the channel width becomes smaller and smaller, with the limit being the cell size Δx. This phenomenon is illustrated in the middle plot of Fig. 4a and in Fig. 4b for a simulation with Δx=0.125 m. In this case, the gap height is seen to exceed 12 m, a situation deemed unphysical in most cases, and we clearly observe that all the flow is routed through a single cell. As can be seen in Fig. 3b, adding the diffusion term to the formulation (Eqs. 9 and 16) enables secondorder convergence of all the variables of interest. In this case, as is evidenced by the bottom plot of Fig. 4a and by Fig. 4b, the channel width can extend over several Δx and the overall shape and aspect ratio of the conduit better fit the idea of what a channel should be.
To demonstrate our AMR implementation, we also perform a numerical convergence analysis using several levels of refinement. The AMR scheme should ideally produce solutions with comparable accuracy to a uniform mesh solution with the same (finest) resolution. Using the previously described test case, a refinement criterion based on the local melting rate enables refinement of the region where channelization occurs. Starting from a baseline, singlelevel simulation with a cell size of Δx, we enable the computation to continue and allow either one or two extra level(s) of refinement, where the finest level will have a cell size of either $\mathrm{\Delta}x/\mathrm{2}$ or $\mathrm{\Delta}x/\mathrm{4}$, respectively. This two or threelevel simulation is then compared to results at a finer resolution, $\mathrm{\Delta}x/\mathrm{4}$ or $\mathrm{\Delta}x/\mathrm{8}$, respectively. The results shown in Fig. 3b indicate that the error using two or three AMR levels is comparable to that of the singlelevel solution with the same effective resolution. The entire plot can be read horizontally, where an imaginary line drawn from a point on the onelevel line should intersect the corresponding equivalent resolution on the twolevel or threelevel line, as is the case here.
Having demonstrated the convergence and efficiency of SUHMO, we turn our attention to a set of simple test cases from SHMIP (de Fleurian et al., 2018) to demonstrate the validity of our implementation. SHMIP is built around six synthetic suites of experiments (labeled from A to F), each consisting of a set of four to six numerical experiments designed to show the formation and evolution of the different drainage elements (sheets and channels) in the context of different input scenarios. Two geometries are considered, a landterminating ice sheet margin and a synthetic valley glacier geometry, as shown in Fig. 2a and b, respectively. In the paper, the results from 13 different models are presented, including the SHAKTI model from Sommers et al. (2018).
We show steadystate cases (Suites A, B and E) in Sect. 5.1 and use Suite F in Sect. 5.2 to explore seasonal forcing. We compare all of our results to those obtained with SHAKTI (all results from the SHMIP project are open source and freely available online; see Werder et al., 2018). We believe that such comparisons are appropriate since SHAKTI is the closest in formulation to SUHMO, and so, by tuning to the same set of parameters and showing that we exhibit similar behavior, we hope to validate our approach. Note that the diffusion term is included in what follows, but its impact is negligible in these nonchannelizing test cases. Additionally, the geothermal flux is not used in Suites A and B, in order to reach a steady state and retrieve the prescribed discharge, and the creep cutoff length scale b_{c} is always set to 0 so that the creep length scale l_{c} reverts to b (see Eq. 7), consistent with the SHAKTI contribution in de Fleurian et al. (2018). Finally, while Suite E has been designed in part to investigate the effect of bed slope on model predictions, we elected to turn off the pressuremelting contribution and focus on validating our implementation on this challenging geometry (our numerical implementation uses a structured, unskewed mesh). Results for Suite E with an added pressuremelting contribution are available in the Supplement.
5.1 Steadystate test cases
Results for the SHMIP test cases A, B and E are presented in Figs. 5, 6 and 7, respectively. The longitudinal evolution of the widthaveraged effective pressure N is displayed on the lefthand side of all the figures. The righthand side shows the total discharge and its various contributions to be compared with the total recharge. The discharge is the evolution along the x axis of the y integral of the facecentered water flux q^{x}, while the volumetric recharge contains contributions from both moulin input (when applicable) and melting. Using the notation introduced in Sect. A2, the discharge at each x_{i} location can be computed on the coarsest grid as
while the recharge is the cellintegrated righthand side of Eq. 1):
At steady state, the recharge and discharge at the domain outlet should be exactly the same. Some figures also display efficient and inefficient contributions to the total discharge, which are computed based on a degree of channelization (DoC) variable. The DoC is a cellcentered variable used to quantify the relative contribution of the two iceopening terms in the RHS of Eq. (9):
with values between 0 and 1 in each cell. A DoC close to 1 indicates a high degree of channelization, while a value close to 0 is indicative of a sheetlike drainage system. The efficient and inefficient contributions to the total discharge at x_{i} are expressed as dis(x_{i})DoC(x_{i}) and dis(x_{i})(1−DoC(x_{i})), respectively.
Suites A and B use the landterminating ice sheet margin geometry (Fig. 2a). In Suite A, a steady and spatially uniform water input is prescribed, with total recharge increasing as we progress from A1 to A6. In Suite B, the same amount of water as in case A5 is fed into an increasing number of moulins (from 1 in B1 to a 100 in B5). For additional details on the parameterization of the different suites, the reader is referred to the online instructions (https://shmip.bitbucket.io/instructions.html, last access: 1 January 2023). Singlelevel (no AMR) simulations are performed, with a fixed cell size Δx^{0}=312.5 m and a fixed time step dt=1 h, which is consistent with values reported from other twodimensional models in de Fleurian et al. (2018). We apply a Dirichlet boundary condition h=z_{b} at the left edge (outlet) of the domain and Neumann boundary conditions with 0 prescribed flux on all other domain boundaries. The steady state is quickly reached in all cases, and simulations are run for approximately 400 d. Results obtained with SUHMO compare well to those obtained with SHAKTI, as expected since both models are built on the same set of equations. We note from Figs. 5b and 6b that SUHMO exhibits slightly smaller contributions from the efficient system in cases that display a hybrid flow configuration (see A6 in Fig. 5b, for example). These discrepancies are attributed to discretization differences between both models. SHAKTI uses an unstructured mesh, while SUHMO uses regular Cartesian meshes. We note that discretization features are probably the cause of the overshoot observed in Fig. 6b for case B1 with SHAKTI. The recharge provided by each moulin in Suite B is another potential source of inconsistency. In this study, a Gaussian distribution is assumed.
Suite E uses the synthetic valley glacier geometry (Fig. 2b). In this experiment, water input is uniformly distributed at the bed of the glacier. Here also, singlelevel simulations are performed, with Δx^{0}=23.4375 m and dt=1 h. The boundary conditions are similar to those in Suites A and B. As already specified, the value of c_{t} is set to 0 for this suite of experiments (removing the pressuremelting term in Eq. 5). The steady state is quickly reached in all cases, and simulations are run for approximately 400 d. While the effective pressure distributions obtained using SUHMO compare well to those obtained using SHAKTI, we note bigger differences in the spatial distribution of the hybrid flow configuration in Fig. 7b. These are again attributed to differences in mesh and cell size between SUHMO and SHAKTI. We note the same tendencies as the overdeepening of the valley increases (from E1 to E5), however, with an increasingly sheetlike distribution throughout.
5.2 Suite F: seasonal cycle with valley topography
We turn our attention to SHMIP Suite F. The results presented in Sect. 5.1 focused on the effects of the geometry and water input type on otherwise steadystate configurations. Suite F prescribes a seasonal water forcing in the synthetic valley glacier geometry of case E1 (Bench Glacier reference geometry). The water input increases from run F1 to run F5, as can be seen in Fig. 8b. The setup follows that of Suite E. A total of 6 years are simulated, allowing sufficient time to settle into a periodic state. Year 6 results are presented in Fig. 8a. Time evolutions of the averaged effective pressure N are extracted at three locations of interest, labeled low, middle, and high bands, depicted in color in Fig. 2b. As before, results obtained with SUHMO closely follow those of SHAKTI.
Overall, while comparisons with SHMIP do not enable a true validation of our results, they do help validate our algorithm and provide an idea of how SUHMO compares to other subglacial hydrology models available in the literature.
6.1 Case description
The test cases in Sect. 5 were all singlelevel experiments. In the present section, we now consider a synthetic square topography of 100 km by 100 km, generated with the intent of emulating catchment areas found at ice sheet margins. Both the bed geometry and ice thickness are shown in Fig. 9a. The bed height varies from 0 m to just under 1000 m, while the ice thickness increases from 100 m in the bottomleft corner (red dot location) to 700 m in the topright corner. Zero flow via a homogeneous Neumann boundary condition is imposed on the two interior boundaries (y=0 and x=100 km), while Dirichlet boundary conditions are prescribed on the other two boundaries (with h=z_{b}, so that P_{w}=0).
Five runs will be presented hereafter – labeled from R_{0} to R_{4} (see Table D1), all using a base mesh with 256 cells in both the x and y directions. This value is chosen so Δx^{0}=390.625 m, which is typical of many ice sheet simulations. The runs are forced by 63 randomly placed moulins, delivering a total water input of 5180 m^{3} s^{−1}. The location of the moulins is shown in Fig. 9b. We emphasize that the topography, moulin location and amount of water input used here have not been designed to reproduce an existing glacier area. The water delivered by the moulins is constant, no seasonal cycle is considered, and it is also quite high: this experiment should be taken as a demonstration of the robust behavior of the system even under prolonged high melting scenarios, when a high degree of channelization is expected. Our purpose is to demonstrate the importance of spatial resolution when looking at subglacial water patterns and ultimately to examine the impact of resolution on the effective pressure distribution.
All runs start from an established, steadystate, singlelevel simulation with no moulins. The moulins activate at t=0 s and the influx ramps up over a period of 2 months until the maximum is reached. The simulations are run for another 22 months after that, bringing the total simulated time to 2 years. The time step is dt=2 h in the first year before increasing to 5 h for the remainder of the simulation.
6.2 Overview of computational requirements
An example of mesh configuration for runs R_{1} to R_{3} is shown in Fig. 10b. In every case, a regridding operation is performed each simulated week, so that the dynamic meshing can follow the water patterns and add or remove refinement around the channels as they develop or retract. We use both the gap heights and melting rates to tag cells for refinement, for this example, cells in which the melt rate exceeds $\mathrm{2}\times {\mathrm{10}}^{\mathrm{5}}$ kg m^{−2} s^{−1} are tagged for refinement, along with cells in which the gap height exceeds 0.1 m (this is effectively only active for R4, since the gap height is smaller for R3 and under). As can be seen, this criterion is very efficient, and only a small area of our entire computational domain ends up requiring up to four levels of refinement (due to the very small cell size reached, R_{4} does not provide any additional visual insight and is therefore omitted from the figures).
A singlelevel run with the same resolution as R_{4} (4096×4096) would require evolving over 16 million cells, when 20 times fewer cells are required in the AMR run to capture what we believe to be most of the important features, as can be seen in Fig. 10a (the cell count in the AMR runs is obtained by averaging the total number of cells evolved over the course of the simulation). The total time to solution is also shown in Fig. 10a. For completeness, results from singlelevel simulations with resolution matching each finest level of refinement k of R_{k} are also reported. The total time to solution in this case does not scale exactly with the increase in resolution: the simulation contains a transient ramp up where moulins activate, and refining means creating more channels and spatial stiffness, as evidenced by the increasing number of FAS MG iterations (not reported here). The timetosolution ratio of single level versus AMR simulations, however, increases with refinement and is well over an order of magnitude for the most refined case presented here. Additionally, note that despite being performed very frequently, the cost of regridding never accounts for more than 0.7 % of the total computational time.
6.3 Results and analysis of the effective pressure distribution
The first three rows of Fig. 11 show fields of gap height (b) and effective pressure (N) for runs R_{0} and R_{2} to R_{4}. A closeup on the valley 1 area (see Fig. 9) illustrates the extent and shape of the central channel. One interesting observation is that no channel forms in R_{0}; no channelization ever occurs, even if the simulation is extended for another 10 years. R_{1} is very similar to R_{0} in that no real channel inception can be seen and is therefore omitted in Fig. 11. This appears to demonstrate a “minimum resolution requirement” to capture channelization behavior in this example. Similar minimum resolution requirements were demonstrated for marine ice sheets in Cornford et al. (2016). Such a minimum resolution is, however, difficult to quantify in practice and may be casedependent. More investigation (currently ongoing) is required to understand tendencies in resolution requirements, but a design criterion could be a much needed numerical tool.
The bottom row of Fig. 11 displays differences in effective pressure (ΔN) between R_{0} and, respectively, R_{2}, R_{3} and R_{4}, from left to right. As expected, N increases significantly in the valley 1 area and near the top boundary (y=100 km) for ${R}_{k,k>\mathrm{1}}$, where channels are seen to develop. For R_{3} and R_{4}, ΔN can reach up to 10 % of the global maximum of N and, more importantly, up to 25 % of the local N value (up to 0.6 MPa). These differences are deemed nonnegligible in the context of evaluating a locally varying friction law, such as the ones from Schoof (2005) or Tsai et al. (2015). Figure 1 in Brondex et al. (2017), for example, shows strong nonlinearities in certain lowpressure regimes, where a difference of this order could result in very different basal drag profiles. Such sensitive areas are more likely to be located near the grounding line, where ice is thinner, and the potential impact on ice velocity warrants further investigation.
Note that all flow features may not be properly resolved in the present configuration. Indeed, even the run with the most AMR levels (R_{4}) still has a resolution of about 25 m (see Table D1), when channels are commonly believed to have a diameter of about 1–10 m. The extent to which channels need to be resolved to sufficiently account for their impact in a coupled subglacial hydrology–ice sheet modeling system is still an open question. For example, are friction law results sufficiently accurate when, say, differences of less than 5 % are observed in the N profile as the resolution is doubled (or a new AMR level is added)? Regardless, this suite of simulations demonstrates that commonly used resolutions (on the order of 300–500 m) can sometimes fail to even reveal channelized features.
In this paper, we present and validate a novel AMR subglacial hydrology model, SUHMO, based on the Chombo framework (Adams et al., 2001–2021). We solve equations similar to those in Sommers et al. (2018), with the addition of a pseudodiffusion to recover the wall melting in channels that was discarded in the derivation of the original equations. We demonstrate the usefulness of this additional term in achieving consistent spatial convergence as finer resolution begins to resolve flow configurations. Our algorithm uses an efficient combination of nonlinear MG iterations embedded in external Picard iterations. We show that results with SUHMO closely follow those obtained with SHAKTI (Sommers et al., 2018) on a broad selection of the SHMIP suite of test cases (de Fleurian et al., 2018). A more complex, multilevel test case is also presented which hints at the need for sufficient spatial resolution to correctly resolve subglacial system dynamics; computational performance analysis demonstrates the efficiency of AMR on such largescale hydrologic problems when compared to a singlelevel run with the same spatial discretization. The AMR approach will eventually enable better icebed boundary conditions for ice sheet simulations at a reasonable computational cost.
With that in mind, future work will focus on the coupling of SUHMO with the BISICLES AMR ice sheet model (Cornford et al., 2013) in order to further investigate the sensitivity of model predictions to spatially accurate modeling of basal conditions. Indeed, while the precise topography of the subglacial network is generally deemed unimportant to the overall ice sheet dynamics, there is to our knowledge no real numerical proof of this assessment, and the results presented here allude to the potential impact of the finer resolution of these systems. A numerical tool capable of resolving the structure of channels, following them as they emerge and disappear, will be an asset in helping to determine whether this is indeed the case. The efficiency provided by AMR can also be leveraged if and when it is determined that sufficient resolution in specific regions is important.
A1 Proper nesting
Calculations are performed on a hierarchy of ℓ_{max} nested, cellcentered level domains. For each AMR level $\mathrm{\ell}=\mathrm{0},\mathrm{\dots},{\mathrm{\ell}}_{\mathrm{max}}$, the problem domain Ω^{ℓ} is discretized by a uniform Cartesian grid Γ^{ℓ} with grid spacing Δx^{ℓ}. Level 0 is the coarsest level, encompassing the entire geometry, while each subsequent finer level, ℓ+1, is a factor ${n}_{\mathrm{ref}}^{\mathrm{\ell}}=\frac{\mathrm{\Delta}{x}^{\mathrm{\ell}}}{\mathrm{\Delta}{x}^{\mathrm{\ell}+\mathrm{1}}}$ finer than level ℓ (${n}_{\mathrm{ref}}^{\mathrm{\ell}}$ is a power of 2, usually 2). Each Ω^{ℓ} is constructed from one or more rectangular subsets of Γ^{ℓ}, as can be seen in Fig. A1a: Ω^{1}, for example, is built from two separate rectangular blocks, each with their own subgrid ${\mathrm{\Gamma}}^{\mathrm{\ell},*}$. An important property is that each domain level is properly nested; that is, no interfaces exist between Ω^{ℓ} and Ω^{ℓ±2}, only between two subsequent levels or the domain boundary.
Certain derived quantities, such as fluxes, are located on two supplementary hierarchies of facecentered level domains that will be denoted Ω^{ℓ,x} and Ω^{ℓ,y} for x and ycentered faces, respectively, on level ℓ.
A2 Cell and facecentered data
Variables can be cellcentered or facecentered. We define a grid vector, p∈ℤ^{2}, choosing to number cells starting at (0, 0), and grid basis vectors ${\mathit{e}}_{x}=(\mathrm{1},\mathrm{0})$ and ${\mathit{e}}_{y}=(\mathrm{0},\mathrm{1})$. Cell centers within Ω^{ℓ} are then located at ${\mathit{x}}_{p}^{\mathrm{\ell}}=\mathrm{\Delta}{x}^{\mathrm{\ell}}(\mathit{p}+\frac{\mathrm{1}}{\mathrm{2}}({\mathit{e}}_{x}+{\mathit{e}}_{y}\left)\right)$ and the midpoints of cell faces within ${\mathrm{\Omega}}^{\mathrm{\ell},*}$ at ${\mathit{x}}_{p}^{\mathrm{\ell}}\pm \frac{\mathrm{\Delta}{x}^{\mathrm{\ell}}}{\mathrm{2}}{\mathit{e}}_{*}$. Cellcentered level variables ${\mathit{\varphi}}_{p}^{\mathrm{\ell}}=\mathit{\varphi}\left({\mathit{x}}_{p}^{\mathrm{\ell}}\right)$ and facecentered level variables ${\mathit{\varphi}}_{p}^{\mathrm{\ell},*}=\mathit{\varphi}({\mathit{x}}_{p}^{\mathrm{\ell}}\frac{\mathrm{\Delta}{x}^{\mathrm{\ell}}}{\mathrm{2}}{\mathit{e}}_{*})$ follow naturally. Notice that ${\mathit{\varphi}}_{p}^{\mathrm{\ell},x}$ is located on the “western” face of the cell p and that ${\mathit{\varphi}}_{p}^{\mathrm{\ell},y}$ is located on the “southern” face of the cell p.
A3 Coarsening operator
We identify cells at different levels which occupy the same geometric regions by means of the coarsening operator ${\mathcal{C}}_{r}\left(\mathit{p}\right)=\frac{\mathit{p}}{r}$ and its inverse, the refinement operator. In that sense, ${\mathcal{C}}_{r}^{\mathrm{1}}\left(\mathit{p}\right)$ is the set of all cells in a grid r times finer that represent the same geometric region (in a finitevolume sense) as the cell p (Martin and Colella, 2000).
A4 Composite variables
Discrete representations of continuous fields are then cell and facecentered composite variables ϕ^{comp}, made up of “valid” (or uncovered by a finer level) portions of the level variables. First, level domains Ω^{ℓ} are divided into valid (${\mathrm{\Omega}}_{\mathrm{valid}}^{\mathrm{\ell}}$) and invalid (${\mathrm{\Omega}}_{\mathrm{invalid}}^{\mathrm{\ell}}$) regions, such that ${\mathrm{\Omega}}_{\mathrm{valid}}^{\mathrm{\ell}}={\mathrm{\Omega}}^{\mathrm{\ell}}{\mathcal{C}}_{{n}_{\mathrm{ref}}^{\mathrm{\ell}}}\left({\mathrm{\Omega}}^{\mathrm{\ell}+\mathrm{1}}\right)$. Valid level domains for facecentered quantities are defined in the same way, ${\mathrm{\Omega}}_{\mathrm{valid}}^{\mathrm{\ell},*}={\mathrm{\Omega}}^{\mathrm{\ell},*}{\mathcal{C}}_{{n}_{\mathrm{ref}}^{\mathrm{\ell}}}\left({\mathrm{\Omega}}^{\mathrm{\ell}+\mathrm{1},*}\right)$. A composite variable is then defined on the union of all valid regions, $\mathbf{\Omega}={\bigcup}_{\mathrm{\ell}}{\mathrm{\Omega}}_{\mathrm{valid}}^{\mathrm{\ell}}$, where ${\mathit{\varphi}}^{\mathrm{comp}}\left(\mathbf{\Omega}\right)={\bigcup}_{\mathrm{\ell}}{\mathit{\varphi}}^{\mathrm{comp}}\left({\mathrm{\Omega}}_{\mathrm{valid}}^{\mathrm{\ell}}\right)$. Likewise, composite vector fields are valid on all faces not overlain by finer faces.
We also construct ghost regions ${\mathrm{\Omega}}_{\mathrm{ghost}}^{\mathrm{\ell}}$ that surround Ω^{ℓ}. These usually contain one or two extra cells and exist purely for numerical convenience – to compute gradients or other facecentered quantities. These buffer regions contain either boundaryspecified values or are used to store extrapolated data or data calculated by interpolation from valid regions of a coarser level. Details pertaining to the computation of composite operators such as gradients and Laplacians can be found in prior publications (Martin and Colella, 2000).
A5 Level variables – averaging down
It is sometimes necessary to transfer information from finer grids to coarser ones: ${\mathcal{C}}_{{n}_{\mathrm{ref}}^{\mathrm{\ell}}}\left({\mathrm{\Omega}}^{\mathrm{\ell}+\mathrm{1}}\right)$ is typically filled from appropriate cellcentered (or facecentered) arithmetic averaging of level ℓ+1 data. An example case is illustrated in Fig. A1b: where levels 0 and 1 meet, facecentered quantity ${\mathit{\varphi}}_{\mathit{p}+\mathrm{1}}^{\mathrm{0},x}$ would be replaced by $\frac{{\mathit{\varphi}}_{\mathbf{2}\mathbf{p}+\mathbf{2}{\mathbf{e}}_{\mathbf{x}}}^{\mathrm{1},x}+{\mathit{\varphi}}_{\mathbf{2}\mathbf{p}+\mathbf{2}{\mathbf{e}}_{\mathbf{x}}+{\mathbf{e}}_{\mathbf{y}}}^{\mathrm{1},x}}{\mathrm{2}}$.
A6 Regridding
We regrid every n_{regrid} time step. n_{regrid} is typically fixed at the beginning of a run. During this process, the solution at each grid cell and on each level, whether valid or currently covered, is tested against some specified criteria (or combination thereof) to determine whether refinement is required, in which case the cell is tagged for refinement. A new set of grids is then generated to ensure all tagged cells are covered by a finer level whilst still satisfying the rules introduced above regarding proper nesting. This procedure enables the refinement or coarsening of the grids over time, following regions of interest. The appropriate refinement criteria vary depending on the type of application. In the case of SUHMO, we typically refine based on high values of the melting rate and/or gap height to ensure we resolve the channelization process.
To relax Eq. (14) on each FAS level, we employ nonlinear Gauss–Seidel method with multicolor ordering. With the formulation of Eq. (15) for the nonlinear operator, and using the notation introduced in Sect. A, we obtain the following discretized equation for each p cell on a given level where the notation ℓ has been omitted:
which can be rewritten as ℋ(h_{p})=0, where
Equation (B2) can be solved by resorting to Newton's method (for a scalar):
where
As mentioned in Sect. 3.1.3, the variable coefficient in the PDE Eq. (16) requires special treatment due to coupling with the main variable h. We tested two different approaches for the treatment of B. In the first, which we will call B fixed, the value of B is evaluated once per Picard iteration and is frozen during the FAS solve. With the second approach, which we will call B on the fly, B is recomputed and averaged down on each MG grid at the beginning of each V cycle. We investigate the overall efficiency of the algorithm in terms of the number of Picard and total FAS iterations and CPU time as a function of the FAS solve tolerance. While keeping the (relative) tolerance of the outer Picard solver at a constant value of 10^{−8}, statistics are collected for about 20 time steps in a transient simulation (using the channelized test case described in Sect. 4.2), and the average per time step is presented in Fig. C1. Across the entire range of FAS tolerances considered, B on the fly is found to be more efficient, with a lower average number of FAS iterations per step and smaller computational time (albeit by a smaller margin due to the small computational overhead of recomputing B and updating the FAS multigrid hierarchy). When using the B fixed approach, tight FAS tolerances lead to a large increase in the average FAS iteration count without having a significant effect on the Picard iteration count. In contrast, B on the fly leads to a significant reduction of the Picard iteration count, because part of the nonlinearity is handled by the FAS.
Based on these results, we opt to use the B on the fly approach for all of our computations, and we fix the FAS tolerances to 10^{−10}, appearing to provide a good compromise in terms of computational time.
We used the released v1.1 version of the publicly available repository of the SUHMO subglacial hydrology model that can be found on GitHub https://github.com/EnnaDelfen/SUHMO (last access: 1 January 2023). SUHMO is written in a combination of C++ and FORTRAN and is built upon the Chombo AMR software framework (Adams et al., 2001–2021). We have registered both SUHMO v1.1 (https://doi.org/10.5281/zenodo.7487485; Felden and Martin, 2022) and the specific forked SUHMO version of Chombo v3.2 (https://doi.org/10.5281/zenodo.7487502; Adams et al., 2022) on Zenodo, including scripts (in Python), inputs and postprocessing tools enabling any user to reproduce all results presented in the paper. Results will be produced either in plain text (for all onedimensional plots presented in this manuscript) or in a HDF5 format (twodimensional plots). For convenience and easy testing of reproducibility, we have added plain text results for all onedimensional plots reported in Sects. 4 and 5.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd164072023supplement.
AMF and DFM formulated the modeling approach. AMF built the SUHMO model on a preexisting framework from DFM. AMF performed simulations and compiled the paper with contributions from DFM. EGN provided overall guidance.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We would like to thank Aleah Sommers and Colin Meyer for many fruitful discussions. Anne M. Felden would also like to thank Mathieu Morlighem and Basile de Fleurian for their help setting up test cases from the SHMIP intercomparison project.
Support for this work was provided through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy (DOE), Office 45 of Science, Biological and Environmental Research and Advanced Scientific Computing Research programs, as a part of the ProSPect SciDAC Partnership. Work at Berkeley Lab was supported by the Director, Office of Science, of the U.S. Department of Energy (contract no. DEAC0205CH11231).
This paper was edited by Ludovic Räss and reviewed by two anonymous referees.
Adams, M., Colella, P., Graves, D. T., Johnson, J. N., Keen, N. D., Ligocki, T. J., Martin, D. F., McCorquodale, P. W., Modiano, D., Schwartz, P. O., Sternberg, T. D., and Van Straalen, B.: Chombo Software Package for AMR Applications – Design Document, Tech. Rep. LBNL6616E, Lawrence Berkeley National Laboratory, https://commons.lbl.gov/display/chombo/Chombo++Software+for+Adaptive+Solutions+of+Partial+Differential+Equations (last access: 1 January 2023), 20012021. a, b, c, d
Adams, M., Colella, P., Graves, D. T., Johnson, J. N., Keen, N. D., Ligocki, T. J., Martin, D. F., McCorquodale, P. W., Modiano, D., Schwartz, P. O., Sternberg, T. D., and Van Straalen, B.: EnnaDelfen/Chombo_3.2: Chombo_SUHMO 1.1 (Chombo_SUHMO_1.1), Zenodo [code], https://doi.org/10.5281/zenodo.7487502, 2022. a
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, https://doi.org/10.1016/S02773791(01)000592, 2002. a
Aschwanden, A., Fahnestock, M. A., Truffer, M., Brinkerhoff, D. J., Hock, R., Khroulev, C., Mottram, R., and Khan, S. A.: Contribution of the Greenland Ice Sheet to sea level over the next millennium, Sci. Adv., 5, eaav9396, https://doi.org/10.1126/sciadv.aav9396, 2019. a
Berger, A. and Loutre, M.F.: An exceptionally long interglacial ahead?, Science, 297, 1287–1288, https://doi.org/10.1126/science.1076120, 2002. a
Berger, M. J. and Colella, P.: Local adaptive mesh refinement for shock hydrodynamics, J. Comput. Phys., 82, 64–84, https://doi.org/10.1016/00219991(89)900351, 1989. a
Bindschadler, R.: The importance of pressurized subglacial water in separation and sliding at the glacier bed, J. Glaciol., 29, 3–19, https://doi.org/10.3189/S0022143000005104, 1983. a
Briggs, W., Henson, V., and McCormick, S.: A Multigrid Tutorial, 2nd edn., SIAM, ISBN 9780898714623, 2000. a, b, c
Brondex, J., Gagliardini, O., GilletChaulet, F., and Durand, G.: Sensitivity of grounding line dynamics to the choice of the friction law, J. Glaciol., 63, 854–866, https://doi.org/10.1017/jog.2017.51, 2017. a
Cornford, S. L., Martin, F. D., Graves, D. T., Ranken, D. F., Le Brocq, A. M., Gladstone, R. M., Payne, A. J., Ng, E. G., and Lipscomb, W. H.: Adaptive mesh, finite volume modeling of marine ice sheets, J. Comput. Phys., 232, 529–549, https://doi.org/10.1016/j.jcp.2012.08.037, 2013. a, b, c
Cornford, S. L., Martin, D. F., Lee, V., Payne, A. J., and Ng, E. G.: Adaptive mesh refinement versus subgrid friction interpolation in simulations of Antarctic ice dynamics, Ann. Glaciol., 57, 1–9, https://doi.org/10.1017/aog.2016.13, 2016. a
Das, S. B., Joughin, I., Behn, M. D., Howat, I. M., King, M. A., Lizarralde, D., and Bhatia, M. P.: Fracture Propagation to the Base of the Greenland Ice Sheet During Supraglacial Lake Drainage, Science, 320, 778–781, https://doi.org/10.1126/science.1153360, 2008. a
de Fleurian, B., Werder, M. A., Beyer, S., Brinkerhoff, D. J., Delaney, I., Dow, C. F., Downs, J., Gagliardini, O., Hoffman, M. J., Hooke, R. L., Seguinot, J., and Sommers, A. N.: SHMIP The subglacial hydrology model intercomparison Project, J. Glaciol., 64, 897–916, https://doi.org/10.1017/jog.2018.78, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o
Felden, E. and Martin, D. F.: EnnaDelfen/SUHMO: v1.1 (v1.1), Zenodo [code and data], https://doi.org/10.5281/zenodo.7487485, 2022. a
Dow, C. F., Werder, M. A., Nowicki, S., and Walker, R. T.: Modeling Antarctic subglacial lake filling and drainage cycles, The Cryosphere, 10, 1381–1393, https://doi.org/10.5194/tc1013812016, 2016. a
Doyle, S. H., Hubbard, B., Christoffersen, P., Young, T. J., Hofstede, C., Bougamont, M., Box, J., and Hubbard, A.: Physical conditions of fast glacier flow: 1. Measurements from boreholes drilled to the bed of Store Glacier, West Greenland, J. Geophys. Res.Earth Surf., 123, 324–348, https://doi.org/10.1002/2017JF004529, 2018. a
Edwards, T. L., Nowicki, S., Marzeion, B., Hock, R., Goelzer, H., Seroussi, H., Jourdain, N. C., Slater, D. A., Turner, F. E., Smith, C. J., McKenna, C. M., Simon, E., AbeOuchi, A., Gregory, J. M., Larour, E., Lipscomb, W. H., Payne, A. J., Shepherd, A., Agosta, C., Alexander, P., Albrecht, T., Anderson, B., AsayDavis, X., Aschwanden, A., Barthel, A., Bliss, A., Calov, R., Chambers, C., Champollion, N., Choi, Y., Cullather, R., Cuzzone, J., Dumas, C., Felikson, D., Fettweis, X., Fujita, K., GaltonFenzi, B. K., Gladstone, R., Golledge, N. R., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huss, M., Huybrechts, P., Immerzeel, W., Kleiner, T., Kraaijenbrink, P., Le clec’h, S., Lee, V., Leguy, G. R., Little, C. M., Lowry, D. P., Malles, J.H., Martin, D. F., Maussion, F., Morlighem, M., O'Neill, J. F., Nias, I., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Radić, V., Reese, R., Rounce, D. R., Rückamp, M., Sakai, A., Shafer, C., Schlegel, N.J., Shannon, S., Smith, R. S., Straneo, F., Sun, S., Tarasov, L., Trusel, L. D., Van Breedam, J., van de Wal, R., van den Broeke, M., Winkelmann, R., Zekollari, H., Zhao, C., Zhang, T., and Zwinger, T.: Projected land ice contributions to twentyfirstcentury sea level rise, Nature, 593, 74–82, https://doi.org/10.1038/s4158602103302y, 2021. a, b
Flowers, G. E.: Modelling water flow under glaciers and ice sheets, P. Roy. Soc. A, 471, 20140907, https://doi.org/10.1098/rspa.2014.0907, 2015. a, b, c
Fricker, H. A., Scambos, T., Bindschadler, R., and Padman, L.: An active subglacial water system in West Antarctica mapped from space, Science, 315, 1544–1548, https://doi.org/10.1126/science.1136897, 2007. a
Gagliardini, O. and Werder, M. A.: Influence of increasing surface melt over decadal timescales on landterminating Greenlandtype outlet glaciers, J. Glaciol., 64, 700–710, https://doi.org/10.1017/jog.2018.59, 2018. a
Gagliardini, O., Zwinger, T., GilletChaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a newgeneration ice sheet model, Geosci. Model Dev., 6, 1299–1318, https://doi.org/10.5194/gmd612992013, 2013. a
Ganopolski, A., Winkelmann, R., and Schellnhuber, H. J.: Critical insolation–CO_{2} relation for diagnosing past and future glacial inception, Nature, 529, 200–203, https://doi.org/10.1038/nature16494, 2016. a
Goelzer, H., Robinson, A., Seroussi, H., and Van De Wal, R. S.: Recent progress in Greenland ice sheet modelling, Current Climate Change Reports, 3, 291–302, https://doi.org/10.1007/s406410170073y, 2017. a
Henson, V. E.: Multigrid methods nonlinear problems: an overview, in: Proc. SPIE, 5016, 36–48, https://doi.org/10.1117/12.499473, 2003. a
Hewitt, I.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sci. Lett., 371–372, 16–25, https://doi.org/10.1016/j.epsl.2013.04.022, 2013. a
Hewitt, I. J.: Modelling distributed and channelized subglacial drainage: the spacing of channels, J. Glaciol., 57, 302–314, https://doi.org/10.3189/002214311796405951, 2011. a
Hewitt, I. J., Schoof, C., and Werder, M. A.: Flotation and free surface flow in a model for subglacial drainage. Part 2. Channel flow, J. Fluid Mech., 702, 157–187, https://doi.org/10.1017/jfm.2012.166, 2012. a, b
Iken, A.: The effect of the subglacial water pressure on the sliding velocity of a glacier in an idealized numerical model, J. Glaciol., 27, 407–421, https://doi.org/10.3189/S0022143000011448, 1981. a
Kamb, B.: Glacier surge mechanism based on linked cavity configuration of the basal water conduit system, J. Geophys. Res., 92, 9083, https://doi.org/10.1029/jb092ib09p09083, 1987. a
Kirkham, J. D., Hogan, K. A., Larter, R. D., Arnold, N. S., Nitsche, F. O., Golledge, N. R., and Dowdeswell, J. A.: Past water flow beneath Pine Island and Thwaites glaciers, West Antarctica, The Cryosphere, 13, 1959–1981, https://doi.org/10.5194/tc1319592019, 2019. a
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res.Earth Surf., 117, F01022, https://doi.org/10.1029/2011JF002140, 2012. a, b
Malczyk, G., Gourmelen, N., Goldberg, D., Wuite, J., and Nagler, T.: Repeat Subglacial Lake Drainage and Filling Beneath Thwaites Glacier, Geophys. Res. Lett., 47, e2020GL089658, https://doi.org/10.1029/2020gl089658, 2020. a
Martin, D. F. and Cartwright, K. L.: Solving Poisson's Equation using Adaptive Mesh Refinement, Tech. Rep. UCB/ERL M96/66, U.C. Berkeley Electronics Research Laboratory, 1996. a
Martin, D. F. and Colella, P.: A CellCentered Adaptive Projection Method for the Incompressible Euler Equations, J. Comput. Phys., 163, 271–312, https://doi.org/10.1006/jcph.2000.6575, 2000. a, b, c
Martin, D. F., Colella, P., and Graves, D.: A cellcentered adaptive projection method for the incompressible Navier–Stokes equations in three dimensions, J. Comput. Phys., 227, 1863–1886, https://doi.org/10.1016/j.jcp.2007.09.032, 2008. a
MassonDelmotte, V., Zhai, P., Pirani, A., Conners, S., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J., Maycock, T., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B.: IPCC, 2021: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, Tech. rep., PM Cambridge University Press, https://www.ipcc.ch/report/ar6/wg1/ (last access: 1 January 2023), 2021. a, b
Nienow, P. W., Sole, A. J., Slater, D. A., and Cowton, T. R.: Recent Advances in Our Understanding of the Role of Meltwater in the Greenland Ice Sheet System, Current Climate Change Reports, 3, 330–344, https://doi.org/10.1007/s4064101700839, 2017. a
Parkinson, J. R., Martin, D. F., Wells, A. J., and Katz, R. F.: Modelling binary alloy solidification with adaptive mesh refinement, J. Comput. Phys., 5, 100043, https://doi.org/10.1016/j.jcpx.2019.100043, 2020. a
Röthlisberger, H.: Water pressure in intraand subglacial channels, J. Glaciol., 11, 177–203, https://doi.org/10.3189/S0022143000022188, 1972. a
Schoof, C.: The effect of cavitation on glacier sliding, P. Roy. Soc. A, 461, 609–627, https://doi.org/10.1098/rspa.2004.1350, 2005. a
Schoof, C.: Icesheet acceleration driven by melt supply variability, Nature, 468, 803–806, https://doi.org/10.1038/nature09618, 2010. a, b
Siegfried, M. R. and Fricker, H. A.: Thirteen years of subglacial lake activity in Antarctica from multimission satellite altimetry, Ann. Glaciol., 59, 42–55, https://doi.org/10.1017/aog.2017.36, 2018. a
Sommers, A., Rajaram, H., and Morlighem, M.: SHAKTI: Subglacial Hydrology and Kinetic, Transient Interactions v1.0, Geosci. Model Dev., 11, 2955–2974, https://doi.org/10.5194/gmd1129552018, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
Stearns, L. A., Smith, B. E., and Hamilton, G. S.: Increased flow speed on a large East Antarctic outlet glacier caused by subglacial floods, Nat. Geosci., 1, 827–831, https://doi.org/10.1038/ngeo356, 2008. a
Tsai, V. C., Stewart, A. L., and Thompson, A. F.: Marine icesheet profiles and stability under Coulomb basal conditions, J. Glaciol., 61, 205–215, https://doi.org/10.3189/2015jog14j221, 2015. a
Tuckett, P. A., Ely, J. C., Sole, A. J., Livingstone, S. J., Davison, B. J., van Wessem, J. M., and Howard, J.: Rapid accelerations of Antarctic Peninsula outlet glaciers driven by surface melt, Nat. Commun., 10, 4311, https://doi.org/10.1038/s41467019120392, 2019. a
Walder, J. S.: Hydraulics of Subglacial Cavities, J. Glaciol., 32, 439–445, https://doi.org/10.3189/s0022143000012156, 1986. a
Weertman, J.: On the sliding of glaciers, J. Glaciol., 3, 33–38, https://doi.org/10.3189/S0022143000024709, 1957. a
Weertman, J.: Catastrophic glacier advances, Tech. rep., Cold Regions Research and Engineering Laboratory (U.S.), http://hdl.handle.net/11681/5909 (last access: 1 January 2023), 1962. a
Werder, M., De Fleurian, B., Beyer, S., Brinkerhoff, D. J., Delaney, I. A., Dow, C. F., Downs, J., Gagliardini, O., Hoffman, M. J., Hooke, R. L., Seguinot, J., and Sommers, A. N.: Subglacial Hydrology Model Intercomparison Project (SHMIP) Data Submissions, ETH Zürich [data set], https://doi.org/10.3929/ETHZB000249168, 2018. 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 Surf., 118, 2140–2158, https://doi.org/10.1002/jgrf.20146, 2013. a, b
Zimmerman, R. W., AlYaarubi, A., Pain, C. C., and Grattoni, C. A.: Nonlinear regimes of fluid flow in rock fractures, Int. J. Rock Mech. Min., 41, 163–169, https://doi.org/10.1016/j.ijrmms.2003.12.045, 2004. a
 Abstract
 Introduction
 Conservation equations for the subglacial drainage system
 Algorithm details
 Analysis of the algorithm efficiency
 SHMIP suite of test cases
 AMR synthetic experiment
 Concluding remarks
 Appendix A: AMR structure and notation
 Appendix B: Relaxation method in the FASMG algorithm
 Appendix C: A case for the treatment of the B coefficient in the FASMG algorithm
 Appendix D: Details of runs R_{0} to R_{4}
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Conservation equations for the subglacial drainage system
 Algorithm details
 Analysis of the algorithm efficiency
 SHMIP suite of test cases
 AMR synthetic experiment
 Concluding remarks
 Appendix A: AMR structure and notation
 Appendix B: Relaxation method in the FASMG algorithm
 Appendix C: A case for the treatment of the B coefficient in the FASMG algorithm
 Appendix D: Details of runs R_{0} to R_{4}
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement