Improved representation of laminar and turbulent sheet flow in subglacial drainage models

Abstract

Subglacial hydrology models struggle to reproduce seasonal drainage patterns that are consistent with observed subglacial water pressures and surface velocities. We modify the standard sheet-flow parameterization within a coupled sheet–channel subglacial drainage model to smoothly transition between laminar and turbulent flow based on the locally computed Reynolds number in a physically consistent way (the ‘transition’ model). We compare the transition model to standard laminar and turbulent models to assess the role of the sheet-flow parameterization in reconciling observed and modelled water pressures under idealized and realistic forcing. Relative to the turbulent model, the laminar and transition models improve seasonal simulations by increasing winter water pressure and producing a more prominent late-summer water pressure minimum. In contrast to the laminar model, the transition model remains consistent with its own internal assumptions across all flow regimes. Based on the internal consistency of the transition model and its improved performance relative to the standard turbulent model, we recommend its use for transient simulations of subglacial drainage.

1. Introduction

The subglacial drainage system beneath the flanks of the Greenland Ice Sheet is subject to seasonal variations in surface melt input, resulting in strong seasonal cycles in subglacial water pressure and ice flow (e.g., Joughin and others, 2008; Moon and others, 2014; Davison and others, 2020; Vijay and others, 2021). The seasonal velocity patterns, and how they vary with increasing volumes of surface melt, are key to understanding ice-discharge-related sea-level contributions from Greenland (e.g., King and others, 2020). However, it remains difficult to model seasonal water pressure and corresponding ice-flow velocities (e.g., Koziol and Arnold, 2018; Cook and others, 2022; Ehrenfeucht and others, 2023) that are consistent with observations of water pressure and ice velocity (e.g., Andrews and others, 2014; Moon and others, 2014; Nienow and others, 2017), limiting the ability of existing models to explain ice-flow patterns and their seasonal variations.

Modern subglacial hydrology models represent water flow through a variety of flow elements, most commonly including efficient drainage through R-channels (Röthlisberger, 1972) and inefficient drainage through linked cavities (Kamb, 1987). Models take different forms (e.g., Flowers, 2015), including those with coupled distributed–channelized flow and spatially extensive channel networks (e.g., Hewitt, 2013; Werder and others, 2013; Hoffman and others, 2018), as well as those comprised of a single set of flow elements that transition between inefficient and efficient drainage (Schoof, 2010; Sommers and others, 2018; Felden and others, 2023). Some models represent physical processes in more detail, for example by including a weakly connected drainage system (e.g., Hoffman and others, 2016), while others trade process detail for computational efficiency (e.g., de Fleurian and others, 2014; Bueler and van Pelt, 2015).

Models that explicitly represent distributed and channelized flow elements (e.g., Hewitt, 2013; Werder and others, 2013; Hoffman and others, 2018) capture much of the presently understood physics of real subglacial drainage and have had success when applied to steady-state ice-sheet hydrology (e.g., Hager and others, 2022), with modelled drainage pathways resembling those inferred from radar data (e.g., Dow and others, 2020). However, these models have difficulty producing realistic water-pressure variations when applied to ice-sheet-scale domains and forced with seasonally varying surface melt inputs. Specifically, models tend to (1) underpredict winter water pressures (de Fleurian and others, 2018; Poinar and others, 2019; Ehrenfeucht and others, 2023) compared to winter water pressure inferred from seasonal velocity patterns (e.g., Moon and others, 2014; Vijay and others, 2021) and observed in borehole water pressures (e.g., Andrews and others, 2014; Wright and others, 2016) (c.f., Downs and others, 2018), (2) fail to capture the late-summer pressure minimum (e.g., Koziol and Arnold, 2018; Cook and others, 2020) that is inferred from typical Greenland outlet glacier velocity records (e.g., Davison and others, 2020), and (3) require a priori assumptions about distributed flow being fully laminar or turbulent (e.g., Hewitt, 2013; Werder and others, 2013). It is unclear whether the assumptions in (3) hold across the typical spatiotemporal domain of these models. Resolving the discrepancies enumerated above is important for capturing the complete relationship between surface melt, subglacial drainage, and ice flow in numerical models.

Most subglacial drainage models require specification of the relationship between water flux or discharge and the hydraulic potential gradient driving flow at the scale of drainage elements. Here we investigate the role of this relationship within distributed drainage components in controlling seasonal pressure variations as modelled with the Glacier Drainage System (GlaDS) model (Werder and others, 2013), a representative example of an explicitly channel-resolving model. We compare seasonal water-pressure variations modelled for different flux models to assess the influence on the shortcomings identified above. On the basis of our results, we make recommendations for the parameterization of distributed water flux in this popular class of channel-resolving drainage models.

2. Methods

2.1 Subglacial hydrology model

Subglacial drainage is modelled with the Glacier Drainage System (GlaDS) model (Werder and others, 2013) as implemented in MATLAB (commit 040032e). GlaDS conceptualizes subglacial water flow occurring through a distributed drainage system composed of linked cavities and through an efficient drainage system composed of R-channels (Hewitt and others, 2012; Schoof and others, 2012; Werder and others, 2013). GlaDS is a representative example of the broader class of multi-component models that share common physical processes (e.g., Hewitt, 2013; Hoffman and others, 2018), and primarily differs in the discrete nature of subglacial channels from models that represent individual elements as transitioning between distributed and channelized flow (e.g., Schoof, 2010; Sommers and others, 2018; Felden and others, 2023).

GlaDS requires specification of a number of parameters that control the formation of subglacial cavities, water flow within distributed and channelized drainage elements, basal sliding, englacial water storage, and the strength of sheet–channel coupling. Constraining drainage model parameters with direct measurements is difficult and has only rarely been done for a few model parameters (e.g., Werder and others, 2009; Pohle and others, 2022). Inferring parameter values via drainage model inversions (e.g., Brinkerhoff and others, 2021; Irarrazaval and others, 2021) is a promising direction given a variety of observational data sources (e.g., Derkacheva and others, 2021; Nanni and others, 2021; Rada Giacaman and Schoof, 2023), however, the limited availability of observational data continues to make parameter inference challenging. In this study, model parameter values (Table 1; Section S1.1) are chosen to obtain summer water pressures near overburden with channels extending ~30 km inland. These values are similar to existing model applications to Greenland-scale catchments with seasonal melt forcing (e.g., Downs and others, 2018; Gagliardini and Werder, 2018; Cook and others, 2022). The size of the controlling bed obstacle (including both the bump height $h_{\rm {b}}$ and the bump length $l_{\rm {b}}$), the width of sheet flow contributing to channel discharge ($l_{\rm {c}}$), and the channel conductivity ($k_{\rm {c}}$) in particular are larger here than typically used for alpine glaciers (e.g., Werder and others, 2013) or steady state Antarctic applications (e.g., Dow and others, 2022; Hager and others, 2022), potentially reflecting the physically larger scale compared to alpine glaciers and the increased size of drainage elements compared to Antarctic applications. Of these key parameters ($h_{\rm {b}}$, $l_{\rm {b}}$, and $k_{\rm {c}}$), the greatest sensitivity is to the bed bump height since it controls the rate of cavity opening and sets the typical sheet thickness.

Table 1. Constants (top group) and model parameters (bottom group) for GlaDS simulations

$^{a}\tilde A$ differs from the canonical rheology parameter A by a factor of ${2\over 27}$. The listed value for $\tilde A$ corresponds to the recommended value A = 2.4 × 10⁻²⁴ s⁻¹ Pa⁻³ for temperate ice (Cuffey and Paterson, 2010).

We intentionally disallow cavities from opening by ice creep when water pressure exceeds ice overburden by setting the ice creep constant $\tilde A_{\rm {s}} = 0$ when $p_{\rm {w}} > p_{\rm {i}}$. We expect that unrepresented physical mechanisms would take over when $p_{\rm {w}}$ exceeds $p_{\rm {i}}$ (e.g., Tsai and Rice, 2010; Schoof and others, 2012; Dow and others, 2015). Based on model experiments, allowing cavities to creep open as a rough approximation of these mechanisms leads to undesirable behaviour: cavities grow arbitrarily large within overpressurized regions, preventing channels from developing and leading to persistent and extensive pressure above overburden. Disabling creep opening is therefore a suitable modelling choice for the configuration presented here (Section S1.2).

While GlaDS is a representative example of a channel-resolving subglacial drainage model, there are physical processes that are missing in its formulation, especially the representation of hydraulically unconnected or weakly connected bed patches (e.g., Murray and Clarke, 1995; Andrews and others, 2014; Hoffman and others, 2016). Since GlaDS represents only hydraulically connected drainage, winter water pressures may be expected to be lower than observations of winter water pressure within disconnected patches. For example, Rada Giacaman and Schoof (2023) report mean winter water pressure ~90% of overburden within hydraulically connected boreholes and >100% of overburden for hydraulically unconnected boreholes for a small alpine glacier.

2.2 Sheet-flow model

The distributed water flux parameterization used by GlaDS and similar models (e.g., Hewitt, 2013; Hoffman and others, 2018) assumes that water flow is exclusively laminar or turbulent. We aim to test the validity of this assumption and develop a parameterization that can transition between laminar and turbulent flow depending on local drainage characteristics. Previously, progress has been made in addressing the shortcomings listed above through adjustments to the distributed drainage parameterization, including representing flow within the distributed drainage system as laminar (Hewitt, 2013; Banwell and others, 2016; Gagliardini and Werder, 2018; Cook and others, 2022), by explicitly parameterizing sheet conductivity as a function of surface melt rate (e.g., Downs and others, 2018) or by including a Reynolds-number-dependent transmissivity (Sommers and others, 2018, 2023).These models share the common feature that resistance to water flow in the distributed drainage system is sensitive to the volume of water supplied from surface and basal melt to the subglacial system. This sensitivity is obtained in different ways, but with similar impacts on the modelled winter water pressure. Therefore, adjusting the distributed water flux parameterization to include both laminar and turbulent flow (e.g., Sommers and others, 2018) is a promising direction to improve modelled seasonal water pressure variations.

2.2.1 Standard sheet-flow model

We consider two primary forms for the distributed water flux parameterization with GlaDS. The standard discharge-per-unit-width (q) parameterization for subglacial drainage models intends to represent the average flux through many sub-grid-scale linked cavities (e.g., Hewitt, 2013; Werder and others, 2013; Hoffman and others, 2018) and can be written

(1)$${\bf q} = -k_{\rm{s}} h^{\alpha_{\rm{s}}}\vert \nabla \phi\vert ^{\beta_{\rm{s}} - 2} \nabla \phi,\; $$

for conductivity $k_{\rm {s}}$, water thickness h, hydraulic potential ϕ, and exponents $\alpha _{\rm {s}}$ and $\beta _{\rm {s}}$.

Choosing values for $\alpha _{\rm {s}}$ and $\beta _{\rm {s}}$ requires an assumption about the relationship between water flux, cavity height, and the hydraulic potential gradient. Values of $\alpha _{\rm {s}} = 3$ and $\beta _{\rm {s}} = 2$ correspond to purely laminar flow (e.g., Creyts and Schoof, 2009; Hewitt, 2013; Cook and others, 2022), while $\alpha _{\rm {s}} = 5/4$ and $\beta _{\rm {s}} = 3/2$ are typically explained as representing fully turbulent flow according to the Darcy–Weisbach relationship (e.g., Schoof and others, 2012; Werder and others, 2013; Hoffman and others, 2018). It is worth noting that the parameterization for channel discharge is written in an analogous way, with the same interpretation of the exponents $\alpha _{\rm {c}}$ and $\beta _{\rm {c}}$.

The validity of the laminar or turbulent assumption can be assessed by inspecting the Reynolds number, Re. In the context of standard fluid dynamics, the Reynolds number predicts whether a specified flow is laminar or turbulent. For a general flow with representative velocity V, length scale D, and for a fluid with kinematic viscosity ν, the Reynolds number is the ratio of the inertial and viscous forces, Re = VD/ν. In the context of the discharge-per-unit-width parameterization (Eqn. (1)), the length scale D is set to the water sheet thickness h, so the Reynolds number becomes Re = q/ν, where q = |q| = VD.

The transition between laminar and turbulent flow is best understood for the simple case of flow through circular pipes. In this case, the empirical relationship between Re and the Darcy friction factor $f_{\rm {D}}$ is summarized by the Moody diagram (Moody, 1944), which demonstrates the clear differences in the behaviour of laminar and turbulent flows (Fig. 1). Laminar flow results in an inverse relationship between Re and $f_{\rm {D}}$ that is independent of roughness (straight line in Fig. 1). Fully turbulent flow is represented by the friction factor being independent of Re as Re → ∞. The transition from laminar to fully turbulent flow can be approximated by the Colebrook–White equation (Colebrook and White, 1937).

Figure 1. Moody diagram, representing the friction factor $f_{\rm {D}} = \frac{h_{\rm {l}}}{\left({L\over D}\right){V^2\over 2g}}$ (for head loss $h_{\rm {l}}$ over a pipe of length L, diameter D, and with flow velocity V), as a function of the Reynolds number Re = VD/ν for different relative roughness scales ($\varepsilon$). The transition region (shaded grey, 1000 ≤ Re ≤ 3000) separates regions of laminar flow and turbulent flow. The laminar friction factor is $f_{\rm {D}} = \frac{64}{\rm {Re}}$ (Moody, 1944), and the friction factor in the transition and turbulent regimes is computed using the Colebrook–White equation (Colebrook and White, 1937).

The fully turbulent behaviour from the Moody diagram can be carried over to the context of distributed subglacial water flow through a macroporous sheet by writing the Darcy–Weisbach equation (e.g., Moody, 1944) for flow between parallel plates and in terms of the flux q instead of the flow velocity (Section S1.3.1, Eqn. S.4). By doing this, fully turbulent flow would require a flow exponent $\alpha _{\rm {s}} = 3/2$, as in the SHAKTI model (Sommers and others, 2018) and in contrast to the assumed value of 5/4 for GlaDS and similar models; however, given the conceptual differences between flow through rough pipes, on which the Moody diagram is based, and the subglacial linked cavity system, we test the sensitivity of modelled water pressure to turbulent flow exponent values $\alpha _{\rm {s}} = 3/2$ and $\alpha _{\rm {s}} = 5/4$. We denote the model using Eqn. (1) with $\alpha _{\rm {s}} = 5/4$ ‘turbulent 5/4,’ with $\alpha _{\rm {s}} = 3/2$ ‘turbulent 3/2,’ and with $\alpha _{\rm {s}} = 3$ and $\beta _{\rm {s}} = 2$ as ‘laminar’ (Table 2). All models use $\beta _{\rm {s}} = 3/2$ to represent turbulent flow.

Table 2. Summary of sheet-flow parameterizations with parameter values substituted in the general forms (Eqns. (1) and (2))

2.2.2 Sheet-flow model with laminar–turbulent transitions

Equation (1) assumes that water flow everywhere and at all times is either purely laminar or purely turbulent. To remove this limitation and develop a model appropriate for the entire Re range, we replace Eqn. (1) with a model that represents both laminar and turbulent flow, with the partitioning governed by the local Reynolds number:

(2)$$-k_{\rm{s}}h^3\nabla\phi = {\bf q} + \omega{\rm{Re}}\left({h\over h_b} \right)^{3 - 2\alpha_{\rm{s}}}{\bf q},\; $$

for bed bump height $h_{\rm {b}}$. Substituting Re = q/ν yields a quadratic equation that can be solved exactly for q (Eqn. S6–S8, Table 2). The transition parameter ω governs partitioning between laminar and turbulent flow, with the transition occurring at approximately Re = 1/ω. The exponent $\alpha _{\rm {s}}$ controls the behaviour of the model in the fully turbulent limit (ωRe ≫ 1).

We call Eqn. (2), which transitions between laminar and turbulent flow based on the local Reynolds number, the ‘transition’ model. In the laminar regime (ωRe ≪ 1), the first term on the right hand side dominates and Eqn. (2) reduces to the laminar model (Eqn. (1) with $\alpha _{\rm {s}} = 3$ and $\beta _{\rm {s}} = 2$). In the turbulent regime (ωRe ≫ 1), the second term on the right hand side dominates and Eqn. (2) reduces to the turbulent model (Eqn. (1) with $\alpha _{\rm {s}}$ specified by the turbulent assumption and $\beta _{\rm {s}} = 3/2$) with an effective turbulent conductivity given by $k_{\rm {t}}^2 = k_{\rm {s}} \frac{\nu }{\omega } h_{\rm {b}}^{3 - 2\alpha _{\rm {s}}}$. In the intermediate regime (ωRe ~ 1), Eqn. (2) smoothly blends laminar and turbulent flow. Table 2 summarizes the five flux parameterizations obtained by applying Eqns. (1) and (2) with turbulent flow exponents $\alpha _{\rm {s}} = 5/4$ and $\alpha _{\rm {s}} = 3/2$.

Figure 2 compares the flux dependence on sheet thickness for the transition (Eqn. (2)), laminar and turbulent models (Eqn. (1)) for a fixed hydraulic potential gradient. The nondimensional sheet thickness, $\tilde h = {h}/{h_{\rm {crit}}}$, is scaled using the critical sheet thickness, defined as the sheet thickness that produces the critical Reynolds number (ωRe = 1). That is, $h_{\rm {crit}}$ is defined to satisfy

(3)$$1 = {\omega\over \nu} k_{\rm{s}} h_{\rm{crit}}^3 \overline{\nabla \phi},\; $$

where $\overline {\nabla \phi }$ is the mean hydraulic potential gradient assuming water pressure is equal to overburden for a given ice geometry. Equation (3) is derived from the laminar model, but with the sheet conductivity chosen for the turbulent model (Section 2.2.3), the critical sheet thickness is identical for laminar and turbulent models. Sheet flux is represented by ωRe = qω/ν, such that values <1 correspond to laminar flow and values >1 represent turbulent flow.

Figure 2. Scaled sheet thickness $\tilde h = {h}/{h_{\rm {crit}}}$ and scaled sheet discharge qω/ν for the five flux parameterizations in Table 2 and with a fixed hydraulic potential gradient. The sheet thickness is scaled by $h_{\rm {crit}}$, the sheet thickness that produces a Reynolds number equal to the transition threshold (ωRe = 1) for turbulent and laminar models.

Transitioning between laminar and turbulent flow in this way means that Eqn. (2) is consistent with the Moody diagram (Fig. 1). The flux is more sensitive to changes in cavity height h and potential gradient $\nabla \phi$ in the laminar regime than in the turbulent regime. By changing the sensitivity to h and $\nabla \phi$ as a function of Re, the transition model (Eqn. (2)) should allow for restricted flow during winter compared to a turbulent model. If the Reynolds number reaches or exceeds the transition point (set by 1/ω), the flux becomes less sensitive to h and $\nabla \phi$, such that the minimum flow resistance (measured by the friction factor $f_{\rm {D}}$) is set by the fully turbulent limit, in contrast to the laminar model where there is no lower bound on the friction factor (e.g., the ‘Turbulent’ region of the Moody diagram; Fig. 1).

The transition parameterization (Eqn. (2)) is similar in form to the Forchheimer equation used for non-Darcy flow through porous media, where the potential gradient is balanced by the sum of a linear term (with respect to flux, or equivalently velocity) representing laminar flow, and a quadratic term representing turbulent flow (e.g., Ward, 1964; Bear, 1972; Venkataraman and Rao, 1998). In the glaciological context, Stone and Clarke (1993) applied the Forchheimer equation to represent drainage within till beneath Trapridge Glacier. The result of Eqn. (2) has a similar effect as the Flowers and Clarke (2002) model, where sheet conductivity is a non-linearly increasing function of water thickness, such that the flux parameterization accommodates a large range in flux magnitudes and approximates both laminar and turbulent flows. Equation (2) is most closely related to the flux parameterization used by the SHAKTI (Sommers and others, 2018) and SUHMO (Felden and others, 2023) models and the rock fracture-flow models these parameterizations are based on (e.g., Zimmerman and others, 2004; Chaudhuri and others, 2013). However, compared to SHAKTI and SUHMO, we apply this parameterization to represent flow exclusively within the distributed drainage system, whereas Sommers and others (2018) and Felden and others (2023) apply a similar parameterization to represent flow within the drainage system as a whole. We have further introduced a free conductivity parameter $k_{\rm {s}}$ to the transition model (Eqn. (2)) in order to recover the standard GlaDS model in laminar and turbulent limits. We retain the standard turbulent flux parameterization for subglacial channels (Werder and others, 2013, Eqn. (12)) based on modelled Reynolds number within the turbulent regime ($\rm {Re} \gtrsim 2000$) for channels discharge above a minimum discharge Q = 10⁻² m³ s⁻¹.

2.2.3 Turbulent model sheet conductivity

The turbulent models in Table 2 prescribe the conductivity $k_{\rm {s}}$ in units that depend on the value of $\alpha _{\rm {s}}$, and differ from the units of $k_{\rm {s}}$ in the laminar and transition models. The conductivity for the turbulent models must therefore be scaled appropriately to obtain a fair comparison between models (Section S1.3.3). For the same reason, implementation of the transition model must be done with caution for models that choose to non-dimensionalize the governing equations (e.g., Werder and others, 2013). The conductivity for the turbulent models, $k_{\rm {t}}$, is computed by setting the turbulent and laminar flux models equal with $h = h_{\rm {crit}}$ (Eqn. (3)) and with the mean hydraulic potential gradient (allowing for $\alpha _{\rm {s}} = 3/2$ or 5/4 for the turbulent model),

(4)$$k_{\rm{t}} h_{\rm{crit}}^{\alpha_{\rm{s}}} \overline{\vert \nabla\phi\vert }^{1/2} = k_{\rm{s}} h_{\rm{crit}}^3 \overline{\vert \nabla \phi\vert }.$$

This scaling choice sets the laminar and turbulent models to intersect at $h = h_{\rm {crit}}$ and ωRe = 1 in Fig. 2. The turbulent models could, instead, be set to match the trajectory of the transition model in the fully turbulent limit. Matching the turbulent trajectories, however, would result in the turbulent models significantly overestimating sheet flux relative to the transition and laminar models for the entire range shown in Fig. 2, rendering the models incomparable. A similar scaling could be done to set the transition model to intersect the laminar and turbulent models at $h = h_{\rm {crit}}$ and ωRe = 1; however, we have chosen to match the laminar and transition models in the laminar regime (the slight offset in Fig. 2 for $\tilde h < 1$ represents the small contribution of the second term in Eqn. (2) and is a consequence of the log-scale).

2.3 Synthetic experiment design

We apply GlaDS with the flux parameterizations in Table 2 to a synthetic ice-sheet margin domain with both synthetic and realistic temperature forcings. The synthetic domain and temperature forcing isolates differences between the models by reducing external controls on the drainage configuration, while the realistic temperature forcing allows us to assess differences in seasonal pressure patterns given plausible variations in surface melt rate that impact the development of efficient drainage in summer.

2.3.1 Domain and geometry

The model is applied to a 100 km × 25 km domain with ice thickness similar to the SHMIP experiment (de Fleurian and others, 2018) (Fig. 3a). The domain is adapted to coarsely represent the K-transect in western Greenland to ensure the surface melt forcing (Section 2.3.2) and geometry are consistent. The bed is flat with an elevation of 350 m, which approximates the ice-margin elevation near the K-transect (Smeets and others, 2018). The minimum ice-surface elevation is 390 m at the terminus (approximately equal to the elevation of the lowest K-transect station; van de Wal and others, 2005). The surface elevation is computed as

(5)$$z_{\rm{s}} = 6 \left(\sqrt{x + 5000} - \sqrt{5000} \right) + 390$$

for x measured in metres from the terminus. The maximum surface elevation is 1909 m, which is near or above the modern-day ELA of >1700 m a.s.l. (Smeets and others, 2018).

Figure 3. Overview of synthetic model domain and moulin distribution. (a) Surface and bed elevation with moulins indicated by black circles. The bands at 15, 30, and 70 km indicate where model variables are aggregated in other figures. (b) Target moulin density (derived from Yang and Smith 2016) and density of randomly generated synthetic moulin design as a function of surface elevation.

2.3.2 Melt forcing

The subglacial model is forced with steady basal melt (0.05 m w.e. a⁻¹, Table 1) and seasonally varying surface melt. The basal melt rate, representing the total melt by external heat sources (i.e., geothermal flux and basal drag), is chosen to be in line with measured (e.g., 3 to 8 cm w.e. a⁻¹, Harper and others, 2021) and inferred basal melt rates (Karlsson and others, 2021) in west Greenland. Since our focus is on seasonal evolution of subglacial drainage, we neglect diurnal variations in surface melt rate. We have found that seasonal water pressure patterns and the relative performance of the flux models (Table 2) are not sensitive to diurnal variations (Fig. S7). Spatially distributed surface melt rates are computed from a prescribed sea-level temperature T ₀(t) using a temperature-index model,

(6)$$\dot m( z,\; t;\; \Gamma) = \max\ \left(0,\; f_{\rm{m}} \left(T_0( t) - \Gamma z\right)\right),\; $$

for melt factor $f_{\rm {m}}$, temperature lapse rate Γ, and elevation above sea level z. The melt factor $f_{\rm {m}} = 0.01\, {\rm m\, w.e.\, a}^{-1}\, ^\circ {\rm C}^{-1}$ is taken from the SHMIP experiment (de Fleurian and others, 2018), and the temperature lapse rate $\Gamma = 0.005\,^\circ {\rm {C}}\, \rm {m}^{-1}$ is chosen to be consistent with summer lapse rates observed in west Greenland (Fausto and others, 2009).

GlaDS is forced with two sea-level temperature timeseries:

1. 1. ‘Synthetic’ forcing using a sea-level temperature parameterization adapted from the SHMIP experiment case D3 (de Fleurian and others, 2018):

   (7)$$T_0( t) = -a\cos\left({2\pi t\over T_{\rm{year}}}\right) + b,\; $$
   where constants a and b control the intensity and duration of surface melt, and $T_{\rm {year}}$ is the number of seconds in a year.

2. 2. ‘KAN’ forcing using daily mean air temperatures recorded at the PROMICE KAN_L weather station (How and others, 2022). We use temperatures from 2014, a representative year in terms of total volume and duration of surface melt over the 2009–2022 period (Section S1.4, Fig. S1)

Prior to applying the above forcings, we forced the model with surface melt identical to that of the SHMIP experiment case D3. Modelled subglacial drainage for the turbulent 5/4 model (as used in the SHMIP experiment) recreates the key features of the published SHMIP outputs (Fig. S9) (de Fleurian and others, 2018).

The constants a and b for the synthetic forcing scenario presented here are computed to retain the same duration of positive sea-level temperatures as the SHMIP experiment and to result in the same total melt volume as the KAN scenario so that only the temporal variations in surface melt rate, and not the total melt volume, vary by scenario. We also tested the sensitivity to total melt volume by increasing the temperatures in the KAN timeseries to produce the same total melt volume as the original SHMIP case D3 (Fig. S10), but present the results for the observed melt volume since these results are expected to be more realistic.

2.3.3 Moulins

Surface meltwater drains into the subglacial system through discrete moulin locations. Supraglacial catchments are generated by randomly placing catchment centroids throughout the domain according to a space-filling maximin design (i.e., a design that maximizes the minimum distance between moulins) and with an elevation-dependent density derived from supraglacial mapping (Yang and Smith, 2016) (Fig. S2). The moulin density is greatest at 1138 m a.s.l., and we assign a total of 68 supraglacial catchment centroids, computed from the product of the observation-derived density and the hypsometry of our domain.

Supraglacial catchments are generated by drawing a Voronoi diagram from the catchment centroids (i.e., assigning each node in the mesh to the catchment of the nearest centroid), and moulins are placed as the node with the lowest surface elevation within each catchment subject to the constraints: (1) the minimum distance between neighbouring moulins is 2.5 km, and (2) moulins can not be placed on boundary nodes or within 5 km of the terminus. Fig. S2 illustrates the moulin and catchment generation scheme in more detail (Section S1.5)

Surface meltwater is accumulated within catchments and instantly routed into moulins. This scheme neglects the impact of supraglacial hydrology, which characteristically delays the diurnal peak and reduces the diurnal amplitude of surface inputs to moulins compared to the diurnal cycle of surface melt rate (e.g., Muthyala and others, 2022). This simplification is appropriate in our synthetic model setup considering the idealized nature of our experiments and since we are not attempting to resolve diurnal cycles in water pressures in response to diurnal variations in moulin inputs.

2.3.4 Boundary and initial conditions

The subglacial model is posed on an unstructured triangular mesh. We apply GlaDS on a mesh with 4156 nodes and a mean edge length of 883 m. This mesh resolution was chosen from mesh refinement tests as a suitable tradeoff between precision and computation time (Fig. S3). Boundary conditions consist of a zero-pressure boundary condition at the terminus (x = 0 km) and a zero-flux condition elsewhere.

GlaDS simulations involve a steady-state spin-up used as initial conditions for periodic runs. The spin-up is accomplished in three phases to ensure numerical stability: (1) 25 years with no surface inputs, starting with a uniform water depth equal to half the bed bump height and no subglacial channels (a sufficient duration for the model to evolve to an intermediate winter-like state that is independent of the uniform initial condition); (2) 25 years with a linear ramp-up of surface melt intensity; and (3) 50 years with constant melt rates to reach a final steady state (evaluated based on the rate-of-change of average water pressure). Each of these steps is longer than strictly necessary. For example, a steady state drainage configuration is typically reached well before the end of step (3), but with implicit and adaptive timestepping the extra spin-up time is associated with negligible increases in runtime.

Periodic simulations are run for two years, and only results from the second year are analysed. It would also be possible to begin seasonal simulations directly from the uniform initial condition. This, however, would require the transient simulations to be run for several melt seasons to reach a periodic state as remnant channels with areas up to S ≈ 12 m² persist through the winter near the terminus and require multiple melt seasons to reach their equilibrium size. Given this nonuniform winter condition, it is faster to approach the periodic state from a channelized system, i.e., from the steady simulation.

3. Results

3.1 Synthetic scenario

To illustrate the differences between modelled water pressure for the five flux parameterizations (Table 2), we first present modelled subglacial water pressure (normalized by ice overburden) and channel discharge for the synthetic forcing scenario (Fig. 4). The primary differences in modelled subglacial drainage are found during winter and above the maximum extent of surface melt (i.e., above ~70 km). The most significant differences are a result of the flux parameterization family (i.e., turbulent, laminar, and transition), with only minor differences related to $\alpha _{\rm {s}}$ (i.e., between turbulent 5/4 and turbulent 3/2, and transition 5/4 and transition 3/2).

Figure 4. Synthetic forcing scenario. Floatation fraction ${p_{\rm {w}}}/{p_{\rm {i}}}$ and channel discharge on 9 July (a-e) for turbulent 5/4 (a), turbulent 3/2 (b), laminar (c), transition 5/4 (d) and transition 3/2 (e) models, and width-averaged floatation fraction on 9 July (f). Width-averaged pressure in bands at x = 15 ± 2.5 km (g), x = 30 ± 2.5 km (h), and x = 70 ± 2.5 km (i) and imposed air temperature at 390 m a.s.l. used to force the temperature-index model (j). The centre of bands used for (g–i) are indicated by vertical lines in (a–f), and the time of (a–f) is shown by vertical lines in (g–i).

These model outputs confirm the well-known winter water pressure problem for the standard turbulent 5/4 model, which tends to produce unrealistically low winter and high summer water pressures (e.g., de Fleurian and others, 2018; Poinar and others, 2019; Ehrenfeucht and others, 2023). For this scenario, the turbulent models predict winter water pressures of 43% of overburden at 30 km with $\alpha _{\rm {s}} = 5/4$ and 51% with $\alpha _{\rm {s}} = 3/2$ (Table 3). These modelled winter water pressures are low compared to borehole observations close to overburden (e.g., winter water pressure higher than 95% of overburden 7 km from the ice margin (van de Wal and others, 2015); ~80–100% of overburden 27 to 33 km from the ice margin (Wright and others, 2016)), even after accounting for the difference in pressure between connected and disconnected bed patches (e.g., Rada Giacaman and Schoof, 2023). The winter water pressure is improved for the laminar, transition 5/4 and transition 3/2 models (each produce water pressure 67% of overburden at 30 km). Summer water pressure is between 81% to 85% of overburden for all models (Table 3). The relative performance of the five models in Fig. 4 is the same as that obtained with higher surface melt rates used in the SHMIP experiment D3 (Fig. S9). The increased melt volume in the SHMIP forcing scenario compared to our synthetic forcing experiment results in summer water pressure above overburden for all except the laminar model.

Table 3. Water pressure normalized by overburden (i.e., floatation fraction) for synthetic and KAN temperature-forcing scenarios. Winter floatation fraction is computed as the average value within x = 30 ± 2.5 km (Fig. 3) during the two months preceding the initial onset of surface melt. Summer floatation fraction is computed as the 95th-percentile width-averaged water pressure produced during the melt season within x = 30 ± 2.5 km. The bracketed number beside summer floatation fractions for the KAN scenario indicates the number of days water pressure exceeded overburden. Water pressure does not exceed overburden in the Synthetic scenario.

The differences in water pressure between the flux parameterizations can be understood by considering the spatial and seasonal pattern in modelled Reynolds number, transmissivity, water depth, hydraulic potential, and conductivity (Fig. 5). The turbulence index (ωRe) highlights regions and times where the turbulent and laminar assumptions are inconsistent (Figs. 5a, b). The turbulent model assumes ωRe ≫ 1 everywhere and for all times, so that the turbulent model is applied inappropriately above x = 20 km and outside of the peak summer season. On the other hand, the laminar model is inappropriate near x = 20 km, near the terminus, and during elevated summer water pressures.

Figure 5. Turbulence index ωRe (log scale; a, b), transmissivity T (log scale; c, d), water depth h (log scale; e, f), potential gradient $\vert \nabla \phi \vert$ (linear scale; g, h), and effective turbulent conductivity (log scale; i, j) on 14 June (left column), and averaged for the band x = 30 ± 2.5 km (right column).

The transmissivity, $T = \rho _w g{q}/{\vert \nabla \phi \vert }$, measures the discharge-per-unit-width associated with a specified potential gradient (Figs. 5c, d). It has similar spatial and seasonal patterns as the turbulence index ωRe. The four order-of-magnitude spatial variation produced by the laminar and transition models is the same as that produced by the SHAKTI model in winter for Helheim Glacier (Sommers and others, 2023).

The spatial and seasonal patterns in turbulence index ωRe can be decomposed into individual contributions from the water depth h (Figs. 5e, f) and potential gradient $\vert \nabla \phi \vert$ (Figs. 5g, h). Of the two components, the water depth h more strongly controls the turbulence index than the potential gradient. This is in line with the mathematically stronger dependence on h than the potential gradient, especially for the laminar and transition models.

The differences in seasonal water pressure variations between the turbulent, laminar, and transition models are largely explained by variations in the effective turbulent conductivity, defined as $k_{\rm {eff}} = q/{h^{5/4}\vert \nabla \phi \vert ^{1/2}}$ (Figs. 5i, j). By this definition, $k_{\rm {eff}} = k_{\rm {s}}$ for the turbulent 5/4 model, meaning that variations in the effective turbulent conductivity for other models allow them to be directly compared to the more commonly used turbulent 5/4 model. For the remaining models, $k_{\rm {eff}}$ is a function of the water thickness and potential gradient, with h again being the main driver based on its higher exponent. The $k_{\rm {eff}}$ for the laminar and transition models varies over two orders of magnitude in space (Fig. 5i) and more than one order of magnitude in time (Fig. 5j). The reduced effective conductivity for the laminar and transition models in winter explains the higher winter water pressure compared to the turbulent models, while the large seasonal changes in effective conductivity explain the reduced seasonal amplitude in water pressure compared to the turbulent models.

The modelled Reynolds number further highlights the conceptual inconsistencies with purely turbulent or laminar models (Fig. 6). Distributed water flow is mostly laminar for all models, with turbulent flow limited to a narrow band that migrates upglacier during the first half of the melt season. The laminar model is inconsistent around the distributed–channelized drainage transition, introducing uncertainty into the onset of channelization predicted with this model.

Figure 6. Reynolds number and channel discharge for synthetic scenario for turbulent 5/4 (a), turbulent 3/2 (b), laminar (c), transition 5/4 (d), and transition 3/2 (e) models.

The results in Figs. 4 and 5 align with what is expected based on the Moody diagram (Fig. 7). Here the spread in the curves for the turbulent 5/4 (lighter blue) and transition 5/4 (lighter yellow) models is a result of the Re–$f_{\rm {D}}$ relationship depending on the hydraulic potential gradient, which varies in space and time. As shown by the effective turbulent conductivity (Figs. 5i, j), the turbulent models have significantly less flow resistance compared to the laminar and transition models except when the Reynolds number approaches the laminar–turbulent transition point. The opposite slope of the turbulent 5/4 cluster of points further suggests a structural problem where flow resistance decreases with decreasing water supply (e.g., flow resistance decreases without bound during winter), regardless of the chosen model parameters. This behaviour is not supported by the other models or the empirical friction factor curves. Of all the models, the transition 3/2 model (darker yellow) is closest to the empirical friction factor curves.

Figure 7. Moody diagram computed from model outputs in the synthetic scenario for the five flux parameterizations (Table 2). The turbulent 3/2 model appears as a horizontal line since its friction factor is independent of Re and $\nabla \phi$. The results from the laminar model are displayed as a thick dashed line to distinguish the modelled results from the theoretical curves. The transition Reynolds number is shown as the solid black line at Re = 2000. For reference, the classical pipe-flow Moody diagram from Fig. 1 is shown in the background (thin black lines, right axis). Note that the scaling between the two axes is arbitrary.

3.2 KAN scenario

The evolution of summer water pressure is sensitive to the temporal pattern of surface melt (Fig. 8). Despite identical total melt volumes between the synthetic and KAN temperature forcing scenarios, peak summer water pressures are higher with KAN temperature forcing (Table 3, S1).

Figure 8. KAN forcing scenario. Floatation fraction ${p_{\rm {w}}}/{p_{\rm {i}}}$ and channel discharge on 9 July (a–e) for turbulent 5/4 (a), turbulent 3/2 (b), laminar (c), transition 5/4 (d) and transition 3/2 (e) models, and width-averaged floatation fraction on 9 July (f). Width-averaged pressure in bands at x = 15 ± 2.5 km (g), x = 30 ± 2.5 km (h), and x = 70 ± 2.5 km (i) and imposed air temperature at 390 m a.s.l. used to drive the temperature-index model (black curve, right axis g–i). The centre of bands used for (g–i) are indicated by vertical lines in (a–f), and the time of (a–f) is shown by vertical lines in (g–i).

The turbulent 5/4 and turbulent 3/2 models once again predict low winter water pressure (44% and 51% of overburden at x = 30 km) compared to the laminar and transition models (67%) (Table 3). Compared to the synthetic forcing case, there is a slightly more prominent late-summer water-pressure minimum observed with the laminar and transition models. At 15 km and 30 km from the terminus, the late-summer water pressure minima are below the winter baseline value, and from September to November the water pressure is generally increasing. This trend is not observed as strongly for the turbulent models given their lower winter baseline water pressure. Peak summer water pressures are broadly similar for all models, exceeding overburden by ~50% in the spring at 30 km. The laminar model has the lowest peak pressure (126% of overburden), perhaps as a result of the breakdown of the laminar assumption (Fig. S5). Following the spring pressure maximum, peak water pressure during melt events remains near overburden, with the 95th-percentile summer water pressure at 30 km between 100% and 104% of overburden (Table 3).

The controls on differences in water pressure between the flux parameterizations are the same as for the synthetic scenario. The opposing sensitivity of the friction factor (i.e., flow resistance) to the Reynolds number (i.e., flow intensity) for the turbulent models compared to the laminar and transition models (Fig. 7) results in significantly lower winter water pressure and a larger variation between winter and summer water pressure for the turbulent models.

To ensure the qualitative differences observed between the synthetic and KAN forcings are not a function of seasonal melt volume, we re-ran the KAN simulations with the original SHMIP D3 (larger) seasonal melt volume. To do this, we increased the temperatures in the KAN timeseries by $2.43^\circ {\rm {C}}$ and adjusted the lapse rate to $\Gamma = -0.0075^\circ \rm {C}\, \rm {m}^{-1}$ to produce the desired seasonal melt volume. The qualitative differences related to the flux parameterizations are robust with respect to this change in total melt volume, however the modelled water pressure is unrealistically high during melt events, reaching ~200% of overburden for all models (Fig. S10).

3.3 Parameter sensitivity

The results for the synthetic (Fig. 4) and KAN (Fig. 8) scenarios represent a single set of parameter values. To see how the performance of these models varies across parameter space, especially as we push the laminar and turbulent models further outside their regimes of applicability, we test parameter settings that represent high and low channel-extent end-members. To achieve this, we tune the sheet ($k_{\rm {s}}$) and channel ($k_{\rm {c}}$) conductivities and $l_{\rm {c}}$, which controls channel initialization (Table S3)

3.3.1 High channelization end-member

To increase the propensity for channelized drainage, the sheet conductivity $k_{\rm {s}}$ is decreased from 0.05 to 0.02 Pa s⁻¹, channel conductivity $k_{\rm {c}}$ is increased from 0.5 to 1.0 m^3/2 s⁻¹, and $l_{\rm {c}}$ is increased from 10 to 25 m. Compared to the reference case (Fig. 8), this sensitivity test results in minimum summer water pressures below the steady winter value for all models. The turbulent model once again produces the lowest winter water pressure, while all models exhibit similar overall behaviour in summer (Fig. S12).

3.3.2 Low channelization end-member

The low channelization end-member uses a sheet conductivity $k_{\rm {s}} = 0.1\, {\rm {Pa\, s}}^{-1}$, channel conductivity $k_{\rm {c}} = 0.2\, \rm {m}^{3/2}\, {\rm s}^{-1}$, and $l_{\rm {c}} = 10\, {\rm {m}}$. These parameters result in channels limited to the lowest 20 km for the turbulent and transition models, with only scattered channels below 15 km produced by the laminar model. All models produce higher peak summer water pressure than the reference case (Fig. 8), with the laminar model producing notably lower pressures closer to overburden (Fig. S14).

Across the sensitivity tests (Figs. S12, S14) and the KAN scenario (Fig. 8), the turbulent model has the behaviour closest to that of the laminar and transition models in the high channelization end-member. However, this scenario produces nearly entirely laminar flow (Fig. S13), rendering the turbulent model physically inconsistent. Similarly, the laminar model outperforms the transition models in the low channelization end-member, where the Reynolds number indicates nearly fully turbulent flow (Fig. S15).

4. Discussion

4.1 Distributed water flux parameterizations

We have presented modelled subglacial drainage configurations for five flux parameterizations (Table 2). With both synthetic and KAN surface melt forcing, the laminar and transition models tend to show more desirable behaviour than the turbulent models. The laminar and transition models result in higher winter water pressures, while still exhibiting late-summer water pressure minima below the winter baseline. These desirable features are more clear in the KAN scenario (Fig. 8), since the smooth melt forcing in the synthetic scenario results in muted seasonal pressure variations (Fig. 4). Given the consistently lower performance of the turbulent model across parameter settings intended to represent extremes of low and high channelization (Figs. S12, S14), these findings do not appear to be a consequence of the particular parameter values used throughout.

The laminar and transition models have generally similar performance in terms of summer and winter water pressure. However, the laminar model produces lower peak water pressures and fewer instances of water pressure above overburden than the transition model. These desirable features produced by the laminar model occur when water flow is demonstrably turbulent, and the laminar model therefore physically inconsistent (Fig. 7). This contradiction points to missing model physics that would act to increase drainage capacity, as accomplished by the laminar model, in the spring to reduce water pressures to closer to overburden.

Unrealistic modelled winter water pressure has previously been addressed using the turbulent 5/4 model by prescribing the sheet conductivity $k_{\rm {s}}$ as a linear function of surface melt rates to allow for reduced conductivity during winter and increased conductivity during summer (Downs and others, 2018). The result of this conductivity parameterization is a similar seasonal pattern of turbulent conductivity as reproduced by the laminar and transition models (Fig. 5j). The major difference between the laminar and transition models and the Downs and others (2018) parameterization is the magnitude of variation. Downs and others (2018) prescribe the conductivity to vary on the order of ${\cal O}( 10^4)$ in time but remain constant in space, whereas we have a variation of order ${\cal O}( 10^1)$ in time, and order ${\cal O}( 10^2)$ in space. These variations in our model results have not been prescribed, but emerge naturally as a result of the flux parameterizations, as has also been demonstrated for the SHAKTI transition parameterization (Sommers and others, 2023).

Seasonal pressure variations have been shown to depend on the evolving connectivity of distributed drainage elements, where low winter water pressure in connected bed regions may by compensated for by high pressure within disconnected bed regions (e.g., Andrews and others, 2014; Hoffman and others, 2016; Rada Giacaman and Schoof, 2023). By comparing a coupled hydrology–dynamics model to sliding speed, subglacial discharge, and ice thickness data from Argentière Glacier, Gilbert and others (2022) found a turbulent flow exponent $\alpha _{\rm {s}} \geq 5$ provided the best fit to observed velocities. The high value for the turbulent flow exponent was interpreted as possibly representing a switch in bed connectivity as a function of the water thickness h (e.g., Flowers, 2000; Helanow and others, 2021). In other words, Gilbert and others (2022) suggest that some of the net effects of changing bed connectivity can be included by increasing the sheet-flow exponent $\alpha _{\rm {s}}$. In this context, some of the poor performance of the turbulent model can be attributed to its failure to represent decreased hydraulic connectivity (i.e., taking $f_{\rm {D}}^{-1}$ as a proxy for connectivity) in winter. Based on these considerations, the possibility that $\alpha _{\rm {s}} > 3$ for sub-turbulent flows, in particular for the transition model, should be investigated if further data suggest that $\alpha _{\rm {s}} > 3$ can reproduce key features related to changes in bed connectivity.

These advantages in the laminar and transition models over the turbulent model come with minimal costs in terms of the difficulty of running the model and in the computational burden. Running the laminar model only requires a trivial change in parameters ($\alpha _{\rm {s}}$, $\beta _{\rm {s}}$, and appropriately scaling the conductivity $k_{\rm {s}}$ in Eqn. (1)). Running the transition model requires a simple modification of the model source code to replace Eqn. (1) with Eqn. (2). The laminar and transition models differ from the turbulent models in terms of computation time by up to ~20%, with the models with $\alpha _{\rm {s}} = 3/2$ taking the longest to run (Table S2).

4.2 Turbulent flow exponent

The turbulent flow exponent ($\alpha _{\rm {s}}$ in Eqns. (1) and (2), Table 1) has a secondary impact on modelled water pressure and drainage configuration relative to the primary control of the form of the flux parameterization. However, winter water pressure for the turbulent model is sensitive to the value of $\alpha _{\rm {s}}$, with the turbulent 3/2 model predicting higher (slightly more realistic) winter water pressure (e.g., Fig. 8h). Sensitivity is very low for the transition model, since the turbulent exponent $\alpha _{\rm {s}}$ only applies in fully turbulent ωRe ≫ 1 limit, which is rarely reached in our model configuration (Fig. 7).

Given that the fully turbulent limit is not reached in our model outputs (Fig. 7), the choice of $\alpha _{\rm {s}}$ for the turbulent and transition models can not be assigned strictly from Darcy–Weisbach pipe flow theory. However, the upwards slope of the envelope of modelled friction factors for the turbulent 5/4 model in Fig. 7 is inconsistent with the other flux models and with empirical friction factor curves, suggesting that $\alpha _{\rm {s}} = 3/2$ is a more reasonable choice than $\alpha _{\rm {s}} = 5/4$.

Our model outputs and theoretical considerations suggest that $\alpha _{\rm {s}} = 3/2$ yields marginally more realistic outputs than $\alpha _{\rm {s}} = 5/4$ (i.e., ~10% higher winter water pressure for comparable parameter values). For modelling studies that make the turbulent flow assumption, we recommend $\alpha _{\rm {s}}$ be treated as an uncertain parameter and tuned where possible (e.g., Gilbert and others, 2022) rather than prescribed as $\alpha _{\rm {s}} = 5/4$ based on precedent. Given the minimal sensitivity for the transition model, and since the turbulent exponent $\alpha _{\rm {s}}$ is only applied in the transition model in the true turbulent limit (ωRe ≫ 1), it should be appropriate to use the transition 3/2 model, instead of transition 5/4, by default.

4.3 Choosing an appropriate flux parameterization

Considering the discussion of both the form (Section 4.1) and turbulent exponent (Section 4.2) of the distributed flux parameterization, we recommend the following:

1. 1. Use the transition 3/2 model by default based on its theoretical (i.e., unlimited Re range of applicability; Fig. 7) and practical (i.e., desirable features in modelled water pressure compared to the turbulent model; Fig. 8) attributes.

2. 2. If only aggregate model outputs (e.g., spatio-temporally averaged basal effective pressure) are important, the laminar model may be appropriate as an approximation of the transition model. In this case, it should be verified that the modelled Reynolds number does not reach the turbulent regime, since the model is physically inconsistent and overestimates sheet capacity with ωRe > 1. However, if modelled water pressures far exceed overburden, it may be practically beneficial to use the laminar model precisely because of this overestimation of the sheet capacity compared to the transition model, while accepting the physical inconsistencies as a symptom of remaining model discrepancy.

3. 3. Avoid the turbulent model for seasonally varying subglacial drainage simulations, unless theoretical (i.e., modelled Reynolds number) and/or practical (i.e., demonstrated sensitivity of quantities of interest to the flux model) reasons are discovered that make its performance superior to the transition model. In this case, the turbulent 3/2 model is recommended over the turbulent 5/4 model, but sensitivity of any quantities of interest to the value of $\alpha _{\rm {s}}$ should be assessed.

4.4 Study limitations

4.4.1 Model geometry and domain

There are a number of limitations related to the idealized model setup utilized here. We have presented results for a flat bed, which is not broadly representative of topography beneath Greenlandic outlet glaciers (e.g., Morlighem and others, 2017). To address this, we additionally tested the sensitivity of model outputs to different realizations of bed topography, including a bed with a ~6 km-wide and 350 m-deep trough both along the centre of the domain and following a sinusoidal path, and U-shaped bed topography (Fig. S17). These tests show no difference in the relative performance of each model since topography has a similar influence on water pressure for all flux parameterizations (Figs. S18, S19, S20). These tests suggest that the performance of the parameterizations are not sensitive to the choice of synthetic bed topography. We have also assessed the sensitivity of our results to boundary conditions by imposing a full floatation boundary condition. Based on this test, the boundary condition imposed at the terminus has a limited impact on water pressure and channel discharge beyond the lowest ~5 km (Fig. S21). Future work should investigate differences in modelled water pressure and the performance of the flux parameterizations when applied to real topography, including for marine-terminating glaciers. Based on these synthetic tests, bed and surface topography or the outlet boundary condition are not expected to change the relative performance of each model.

4.4.2 Meltwater forcing

The synthetic and KAN surface melt forcings and uniform basal melt rate used here are simplified relative to realistic melt rates. The KAN surface meltwater forcing scenario, derived from daily mean temperatures recorded at the KAN_L AWS (How and others, 2022) and a simple temperature-index model, elicits a more realistic water pressure response than the unrealistic synthetic scenario. For the KAN forcing, we have used a single sample of temperature forcing measured at the KAN_L PROMICE station. This timeseries was chosen to be representative in terms of the total melt volume and melt season duration, however different temperature timeseries will result in different modelled water pressures. Given the consistency of the differences between the flux parameterizations between the synthetic (Fig. 4) and KAN scenarios (Fig. 8), it is unlikely the performance differences of the flux parameterizations are a function of the choice of temperature timeseries.

Since we are focussed on subglacial water pressure on seasonal timescales, we have chosen to omit diurnal variations in forcing the subglacial drainage model. We have also ignored supraglacial (e.g., Hill and Dow, 2021; Poinar and Andrews, 2021) and englacial (e.g., Andrews and others, 2022) hydrologic processes that impact the diurnal evolution of water pressure (e.g., Andrews and others, 2018). Neglecting diurnal oscillations has previously been shown to have only a limited impact on the seasonal development of the subglacial drainage system (e.g., Werder and others, 2013), and experiments with prescribed diurnal forcing show a minimal impact (Fig. S7).

We have prescribed a uniform basal melt rate of 0.05 m w.e. a⁻¹ in our simulations. Since GlaDS does not allow for melt by potential energy dissipation within the distributed drainage system (Werder and others, 2013), this rate represents melt by external heat sources including geothermal flux and basal drag (e.g., Harper and others, 2021; Karlsson and others, 2021). Based on tests with a lower basal melt rate, representing a low friction regime, the winter water pressure decreases for all models (Fig. S11), such that our specific conclusions do not depend on the chosen basal melt rate. The increased winter water pressure we observe with higher basal melt rates corroborates previous findings that accurate accounting of the basal heat balance (e.g., Harper and others, 2021; Karlsson and others, 2021) is important for reproducing high winter water pressure (Sommers and others, 2023).

4.4.3 Reynolds number and transition parameter

The partitioning between laminar and turbulent flow (Fig. 5) has been based on the Reynolds number computed using the distributed flux q, which represents the average flux through many subglacial cavities within each model element. It is therefore not exactly clear how representative this bulk-averaged Re metric is of flow through physical subglacial drainage elements comprising the ‘distributed water sheet’ as represented in models. The problem of determining a representative Reynolds number is shared by models of non-Darcy porous flow (e.g., Ward, 1964; Bear, 1972; Venkataraman and Rao, 1998). In this context, the problem can be partially addressed by direct numerical simulation of flow through a particular medium (e.g., Wood and others, 2020). Given the uncertainty in the exact form of subglacial drainage elements, this is not a question that can be answered within the framework of current subglacial hydrology models, but it is important to consider when assigning the transition parameter ω, since the Reynolds number cannot be interpreted as precisely as for simple flows. We have assumed that the transition from laminar to turbulent flow occurs at Re ≈ 2000, but it remains to be shown what transition threshold yields the best agreement with velocity or subglacial water pressure data in more realistic model settings.

5. Conclusions and recommendations

Subglacial drainage models are key to understanding the relationship between surface and basal melt, basal motion, and ultimately grounded-ice contributions to sea level (e.g., King and others, 2020). However, these models have important shortcomings when applied to ice-sheet-scale domains with seasonally varying melt forcing. Subglacial models (1) underpredict winter water pressures, (2) fail to capture the late-summer pressure minimum, and (3) require a priori assumptions about distributed flow being fully laminar or turbulent, among other limitations. We have demonstrated that these three problems can be measurably addressed by modifying the parameterization controlling water flux in the distributed (linked-cavity) drainage system while maintaining purely turbulent flow within subglacial channels.

We have tested five flux parameterizations (Table 2), including the standard turbulent model (e.g., Schoof and others, 2012; Werder and others, 2013), the fully laminar model (e.g., Hewitt, 2013; Gagliardini and Werder, 2018; Cook and others, 2022), and a parameterization that transitions between laminar and turbulent flow based on the local Reynolds number, for two values of the turbulent flow exponent ($\alpha _{\rm {s}} = 5/4,\; \, 3/2$) where appropriate. The flux parameterizations are tested within the GlaDS model (Werder and others, 2013) using synthetic and realistic seasonally varying air temperature forcing on a synthetic ice-sheet margin domain.

Laminar and transition models outperform turbulent models on all identified criteria. Winter water pressure is increased by 18–33% of overburden across the domain by using the laminar and transition models for comparable parameter values. When forced with realistic seasonally varying air temperatures, the laminar and transition models produce late-summer water pressures slightly below the winter baseline. Fundamentally, the turbulent and laminar models are inconsistent with their underlying assumptions over the full range of Reynolds numbers that must be represented by the sheet-flow parameterization (e.g., Fig 7).

We suggest using the transition ($\alpha _{\rm {s}} = 3/2$) model where possible based on its desirable features and physical consistency in representing flows over a realistic range of Reynolds numbers. The laminar model produces similar results for seasonal-scale simulations, but suffers from conceptual inconsistencies. However, it may be beneficial to use the laminar model despite its conceptual limitations when modelled spring water pressure is unexpectedly high. The turbulent model should be used with caution and an appreciation of its structural limitations.

The practical and conceptual improvements made by the transition model are encouraging for modelling transient subglacial water-pressure variations. However, a gap remains between models and observations, especially for winter water pressure and in explaining why the physically inconsistent laminar model sometimes performs best in the scenarios we have presented. Promising areas to direct efforts to improve modelled winter water pressure include the representation of spatiotemporal heterogeneity in basal hydraulic connectivity (e.g., Andrews and others, 2014; Hoffman and others, 2016; Helanow and others, 2021) and transmissivity (e.g., Gilbert and others, 2022) in subglacial models, and establishing routine two-way hydrology–dynamics coupling (e.g., Cook and others, 2022). Further aspects of the physics captured by subglacial drainage models that are open for development, and should lead to improved realism in model outputs, include coupling basal hydrology with englacial (e.g., Koenig and others, 2014; Andrews and others, 2022) and supraglacial processes (e.g., Das and others, 2008; Hill and Dow, 2021); the transition to open-channel flow when water pressure becomes negative (e.g., Röthlisberger, 1972; Hewitt and others, 2012; Sommers and others, 2018); representing spatially variable basal materials (e.g., Muto and others, 2019; Maier and others, 2021) and the corresponding impact on the appropriate cavity opening parameterization and the form of distributed water flow; and the physics of over-pressurization, including the mechanical response of the ice overhead (e.g., Tsai and Rice, 2010) and the englacial and supraglacial hydrologic implications (e.g., St Germain and Moorman, 2019; Andrews and others, 2022).