Interactive comment on “ A confined – unconfined aquifer model for subglacial hydrology and its application to the North East Greenland Ice Stream ” by Sebastian

In ‘A confined-unconfined aquifer model for subglacial hydrology and its application to the North East Greenland Ice Stream’, the authors develop a new model of subglacial hydrology based upon the idea that subglacial water flow can be approximated as a layer of variably permeable material (a linked cavity system or an actual till aquifer) topped with a mostly inpermable cap (the glacier). The paper makes the claim that this approach has lower computational requirements than more explicit models and can

arborescent networks; they usually develop over the summer season when a lot of melt water is available. It is assumed that these channelized or efficient drainage systems able to drain large amounts of water in short time spans are predominant in alpine glaciers and on the margins of Greenland, where substantial amounts of surface melt water are capable of reaching the bed (van den Broeke et al., 2017). In the interior of Greenland and also in most parts of Antarctica, the water supply is limited to melt due to the geothermal and frictional heating within the ice (Aschwanden et al., 2016) -a circumstance favoring 5 distributed systems.
Seasonal variations of ice velocity have been observed and attributed to the evolution of the drainage system switching between an efficient and inefficient state in summer and winter (Bartholomew et al., 2010). For this reason, a new generation of subglacial drainage models has been developed recently that is capable of coupling the two regimes of drainage and reproducing the transition between them (Schoof, 2010;Hewitt et al., 2012;Hewitt, 2013;Werder et al., 2013;De Fleurian et al., 2014; 10 Hoffman and Price, 2014). While these models demonstrate immense progress for modeling spontaneously evolving channel networks, it is still a challenge to apply them on a continental scale. A comprehensive overview of the various operational and newly emerging glaciological hydrology models is given in Flowers (2015).
Distributed or sheet structures can naturally be well represented using a continuum approach, while channels usually require a secondary framework, where each feature is described explicitly. Water transport in channels is a complex mechanism that 15 depends on the balance of melt and ice creep (Nye, 1976;Rothlisberger, 1969), channel geometry, and network topology.
Additionally, the network evolves over time which further complicates modeling of this process. When simulating channel networks particular care must be also taken to prevent the emergence of instabilities due to runaway merging of channels (see the discussion in Schoof et al. (2012)). This leads to increased modeling complexity and high computational costs. An exception to this is the work of De Fleurian et al. (2014), where both systems are represented by Darcy flow through separate porous 20 media layers. The layer representing the channels has its parameters (namely hydraulic conductivity and storage) adjusted to exhibit the behavior of an effective system.
We take this idea even further and only use a single layer of Darcy flow with locally adjusted transmissivity of the layer at locations where channels form, while, at the same time, accounting for a decrease in storage. This means that we approximate the channel flow as a fast diffusion process similarly to work in De Fleurian et al. (2014); however, a single Darcy flow layer with 25 spatially varying parameters (effective hydraulic conductivity) accounts for both drainage mechanisms. Similar approaches are known to have been applied to modeling of fracture networks in rock Van Siclen (2002). This reduced complexity model does not capture channels individually but represents their effect by changing specific local properties. Since our model aims to simultaneously represent the main properties of both drainage mechanisms (efficient and inefficient), special care must be exercised when choosing the model parameters and relating them to the physical properties of a specific scenario. While this 30 strategy may not help to advance the precise understanding of channel formation processes, it captures the overall behavior and allows to examine the complex interactions on larger spatial and temporal scales. Currently we focus on channelised hydrology, whereas in future we will consider to expand the model towards cavity formation and closure. This is primarily due to the fact that cavity formation is not considered a dominant effect underneath the large ice sheets (Fowler, 1987).
In addition, we introduce a new Confined-Unconfined Aquifer Scheme (CUAS) that differentiates between confined and unconfined flow in the aquifer (Ehlig and Halepaska, 1976). While the assumption of always saturated -and therefore confined -aquifers may be true for glaciers with large water supply, it does not hold in areas with lower water input. Especially in locations far from the coast, the water supplies are often insufficient to completely fill the aquifer. Ignoring this leads to significant errors in the computed hydraulic potential and unphysical, i.a. negative, water pressure. This problem has been 5 analyzed in detail by Schoof et al. (2012), but here we study the effect in the context of equivalent aquifer models using unconfined flow as a possible solution.
Large scale ice flow models often compute the basal velocity using a Weertman-type sliding law, where the inverse of the effective pressure (difference between ice overburden pressure and water pressure) determines the velocity at the base.
Low effective pressure leads to high basal velocity. Without subglacial hydrology models, the ice models simply take the 10 ice overburden pressure as effective pressure completely neglecting water pressure. This is a major reason why these models struggle to represent fast flowing areas such as ice streams. The effective pressure computed by our model can be easily coupled to an ice sheet model and improve results for fast flowing areas.
The following work is structured as follows. In the next section, we present the one-layer model of subglacial aquifer and briefly describe its discretization. In Sect. 3 the model is applied to artificial scenarios, and the sensitivity to model parameters 15 and stability are investigated. In addition, results for seasonal forcing are presented there, and we show how the model evolves over time. Section 4 demonstrates the first application of the proposed methodology to the North East Greenland Ice Stream (NEGIS), which is the only interior ice stream in Greenland. It penetrates far into the Greenland mainland with its onset close to the ice divide, so sliding apparently plays a major role in its dynamics. A short conclusions and outlook section wraps up the present study.

Confined-Unconfined Aquifer Scheme
The vertically integrated continuity equation in combination with Darcy's law leads to the general groundwater flow equation (see e.g. Kolditz et al. (2015)): 25 with S the storage coefficient (change in the volume of stored water per unit change of the hydraulic head over a unit area), h the hydraulic head (water pressure in terms of water surface elevation above an arbitrary datum; piezometric head), T transmissivity of the porous medium, and Q the source term. For a confined aquifer, T = Kb, where K is the hydraulic conductivity, and b is the aquifer thickness. S = S s b with specific storage S s given by with material parameters for the porous medium (porosity ω, compressibility α) and water (density ρ w , compressibility β w ).
In order to consider the general form covering both cases (confined and unconfined), we follow Ehlig and Halepaska (1976) and write the general form for the confined-unconfined problem: Now the transmissivity and the storage coefficient depend on the head and are defined as is the the local height of the head over bedrock z b and effective storage coefficient S e is given by with This means that as soon as the head sinks below the aquifer height, the system becomes unconfined, and therefore only the 10 saturated section contributes to the transmissivity calculation. This also prevents the head from falling below the bedrock as detailed in Section 3.2. Additionally, the mechanism for water storage changes from elastic relaxation of the aquifer (confined) to dewatering under the forces of gravity (unconfined). The amount of water released from dewatering is described by the specific yield S y . Since this amount is usually orders of magnitudes larger than the release from confined aquifer (S y S s b), it is useful to introduce a gradual transition as in Eq. (6) controlled by a user defined transition parameter d. 15 Note that the transmissivity is not homogeneous making Eq. (3) nonlinear. This fits with our approach to describe the effective system (channels) by locally increasing the hydraulic conductivity. The benefit of this approach is discussed in Sect. 3.2.
Water pressure P w and effective pressure N are related to hydraulic head as and with g acceleration due to gravity, P i = ρ i gH the cryostatic ice overburden pressure exerted by ice with thickness H and density ρ i .

Opening and closure
Opening and closure of channels is governed by melt at the walls due to the dissipation of heat and the pressure difference 25 between the inside and outside of the channel leading to creep deformation. We follow de Fleurian et al. (2016) in using the classical channel equations from Nye (1976) and Röthlisberger (1972) to scale our transmissivity in order to reproduce this behavior. In contrast to de Fleurian et al. (2016), we do not evolve the aquifer thickness but the conductivity K, which leads to 5 and v creep = 2An −n |N | n−1 N K with L the latent heat, r roughness factor, A the creep rate factor depending on temperature, and n the creep exponent, which we choose as n = 3. Depending on the sign of N , creep closure as well as creep opening can occur. Negative effective pressure over prolonged time is usually considered unphysical, and the correct solution to this would be to allow the ice to separate from is still applicable because this is how a channel would behave for N < 0. In Sect. 3.1, we test the sensitivity of K and N to the magnitudes of r and A.

Experiments with artificial geometries
Testing out equivalent layer model and finding parameters for it is not straightforward because there are no directly comparable 15 physical properties. Moreover, observations and measurements of subglacial processes are in general difficult and sparse.
We address this by testing the model with some of the benchmark experiments of the Subglacial Hydrology Model Intercomparison Project (shmip.bitbucket.io).
The proposed artificial geometry mimics a land-terminating ice sheet margin measured 100 km in the x-direction and 20 km in the y-direction. The bedrock is flat (z b (x, y) = 0 m) with the terminus located at x = 0 while the surface z s is defined by a square root function z s (x, y) = 6 (x + 5e3) 1/2 − (5e3) 1/2 + 1. Here, we use the SHMIP/B2 setup, which includes 10 5 moulins with constant in time supply. Boundary conditions are set to zero influx at the interior boundaries (y = 0, y = 20, x = 100) and zero effective pressure at the terminus. All experiments start with initial conditions that imply zero effective pressure and are run for 50 years to ensure that they reach a steady state.

Parameter estimation and sensitivity
SHMIP is primarily intended as a qualitative comparison between different subglacial hydrology models, where results from 10 the GlaDS model (Werder et al., 2013) serve as a "common ground". Here, we use it as a basis for an initial tuning and a study of the sensitivity of our model with regard to parameters. The upcoming results from the SHMIP are also the reason why we do not show a comparison to other models in this study but refer to the manuscript in preparation instead. In Table 1, we show the physical constants used in all setups and runs. The values in the lower half are properties of the porous medium and are only estimated. Since they are utilized in the context of the equivalent layer model this is not an issue. 15 Table 2 contains the model parameters in the upper part and the variables computed by the model in the lower part.
We divide the sensitivity analysis into a general block investigating the sensitivity to the amount of water input into moulins, the layer thickness b, the confined / unconfined transition parameter d, grid resolution dx (Fig. 3) and a block that examines the parameters directly affecting channel evolution such as channel roughness factor r, creep rate factor A, and bounds for the allowed conductivity K min and K max (Fig. 5). In Table 3, we list values that lead to the best agreement with the SHMIP 20 benchmark experiments and thus are used in the following as the baseline for our sensitivity tests.
In Figs. 3a and b, the model's reaction to different amounts of water input through the moulins is shown. With deactivated conductivity evolution (K = const., dashed lines), larger water inputs lead to higher water pressure, hence lower effective  (2014) pressure N . In this case, a moulin input of 18 m 3 s −1 leads to negative values of N . With activated evolution of K, the conductivity adapts to the water input: as more water enters the system through moulins, the conductivity rises. Vertical gray bars show the location of moulins along the x-axis, and the most significant increase in K occurs directly downstream of a moulin. This happens because the water is transported in this direction leading to increased melt. At the glacier snout (x = 0), the ice thickness is at its lowest so almost no creep closure takes place; hence, the conductivity reaches its allowed maximum 5 of 0.5 ms −1 for all tested parameter combinations. Significant development of effective drainage is visible for inputs above 2.25 m 3 s −1 (yellow line). The resulting effective pressure decreases with rising water input as the system becomes more efficient at removing water. Up to ca. 35 km distance from the snout this results in almost identical values of N for all forcings above 2.25 m 3 s −1 . The system adapts so that it can remove all of the additional water efficiently. In Figs. 3i and j, the twodimensional distributions of N and K are shown for the baseline parameters; in addition, we detail the temporal development 10 of K in Fig. 4.
"Channels" (indicated by regions of high conductivity) form downstream from moulins and continue straight towards the ocean. The effective pressure drops around water inputs and along the "channels".
The layer thickness b directly influences the transmissivity according to Eq. (4). A thicker layer can transport more water out of the system and therefore leads to lower water pressure and higher effective pressure (Figs. 3c and d). At the same time, 15 a thinner layer leads to increased water pressure and higher pressure gradients which results in higher v melt and can induce negative effective pressure and creep opening.  Table 3. Selected baseline parameters for all experiments unless otherwise noted. These parameters best match the SHMIP targets.

Name Value Units
Kmin The confined-unconfined transition parameter d does not show noticeable effects on the results (Figs. 3e and f) because the experiment has sufficient water input so that all cells are confined in the steady state.  and h). However, coarse resolutions are unable to resolve the steps that appear at the moulins.
In Figs. 5a and b, we show the results for different values of K min . K min ensures that the system is never completely watertight and determines how much water can pass through independently of the system's current state. Higher values lead to more water transport, lower water pressure, and higher effective pressure. In this experiment, values higher than 0.01 m s −1 5 show no reactions of K to moulins upstream of 10 km.
K max (Fig. 5e and f) has no visible impact on the resulting pressure distribution. The differences in K are only large at the snout, where the maximum is always reached because the ice thickness is too low to counteract the melt term. No significant influence on the upstream area has been detected.
The creep rate factor A determines the "softness" of the ice and therefore effects the creep term in Eq. (9). Larger values of 10 A imply warmer ice; hence, more creep closure (see Figs. 5e and f). Note, that this also effects creep opening if N < 0. The roughness factor r is meant as a measure for small bumps and imperfections, where a rougher "channel" (larger r) would experience more melt because of the larger contact area and more turbulent mixing. Larger values of r lead to higher conductivity and more water transport resulting in lower P w and higher N .
3.2 The benefit from treating unconfined aquifer As described above, the confined-unconfined aquifer approach is advantageous for obtaining physically meaningful pressure distributions. In the example illustrated in Fig. 6, we use a slightly modified geometry, where the bedrock rises towards the upstream boundary forming a slab z b (x, y) = max (3((x + 5e3) 1/2 − (5e3) 1/2 ) − 300, 0). The supply is constant in time and space, and we choose a low value of 7.93e-11 m/s (≈2.5 mm/a) to compare our improved scheme to the simple confined only 5 case. Fig. 6 shows a comparison of the steady state solutions: For the confined-only case, the hydraulic head drops below the bedrock at the upstream region. This results in negative water pressure for these regions. Addressing this by simply limiting the water pressure to zero would result in inconsistencies between the pressure field and the water supply. Our new scheme limits the transmissivity when the head approaches the bedrock and by this means ensures p w ≥ 0 in a physically consistent way. Additionally, the confined-only solution completely depends on boundary conditions and supply terms, basal topography 10 has no influence in this case (apart from governing dK/dt). The possibility of the aquifer to become unconfined captures the expected behaviour much better: At high water levels, water pressure distribution dominates water transport, while at low levels the bed topography becomes relevant.

Seasonal channel evolution and properties
In order to understand our model's ability to simulate the seasonal evolution of subglacial systems, we selected the setup 15 SHMIP/D and ran it with different values of key model parameters. This experiment does not include any moulins but prescribes a non-uniform spatial distribution of supply instead that also varies seasonally. A simple degree day model with varying temperature parameter dΘ provides water input rising from the downstream end (lowest elevated) of the glacier towards the higher elevated areas over summer: Q dist (z s , t) = max (0, (z s LR + Θ(t))DDF) + Q basal .
Here, yr = 31536000 s denotes the number of seconds per year, LR = −0.0075 Km −1 the lapse rate, DDF = 0.01/86400 mK −1 s −1 is the degree day factor, and Q basal = 7.93 × 10 −11 ms −1 is additional basal melt. The resulting seasonal evolution of the sup-5 ply is shown in Fig.7a. The model is run for 10 years so that a periodic evolution of the hydraulic forcing is generated. Here, we present the result for one parameter set only since the model is not very sensitive in this setup. For the prescribed water supplies, the conductivity only reaches its lower limit K min for very short periods, and the upper limit K max is never reached at all (apart from the terminus, where it does not have any significant influence on the upstream behavior of the system as we have shown in the previous experiment). 10 We  The observed behavior is expected and indicates that our model is able to represent the seasonal evolution of the subglacial water system. Increasing water supply over the year leads to rising water pressure and dropping effective pressure. When the conductivity rises in response, the effective pressure goes up again despite the supply not yet falling again because the more efficient system is able to transport the water away. For the cases, where no visible change in K occurs such as dΘ = −6 (blue line in Fig. 7b), the effective pressure directly follows the supply at the terminus, while at the center position (dΘ = −2, cyan 30 line, Fig. 7c), the minimum is offset by the time needed for the supply to reach that location. The maximum in conductivity K is reached later because once the system becomes efficient, increased water transport stimulates melting that opens the system even more. This self-reinforcing process is only stopped when enough water is removed and the reduced water flux reduces the melt again. We assume that this leads to similar locations of the conductivity maxima for different dΘ and the resulting similar reemerging of winter conditions in N .
In this experiment, N becomes negative during the seasonal evolution, which is not physically meaningful. We attribute such behavior to a lack of adjustment of water supply to the state of the system. In reality, the supply from runoff or supraglacial drainage would cease as soon as the pressure in the subglacial water system becomes too high; here we simply continue to 5 pump water into the subglacial system without any feedback. This then leads to negative values of N . It is also consistent with the finding that N becomes negative earlier in the season in cases of higher supply. This deficiency will be addressed in future work.

Subglacial hydrology of NEGIS, Greenland
The role of subglacial hydrology in the genesis of ice streams in general is not well understood yet. NEGIS is a very distinct 10 feature of the ice sheet dynamics in Greenland; thus, the question about the role of subglacial water in the genesis of NEGIS is critical. The characteristic increase in horizontal velocities becomes apparent about 100 km downstream from the ice divide (Vallelonga et al., 2014). Further downstream, the ice stream splits into three different branches: the 79 • North Glacier (79NG), Zacharias Isbrae (ZI), and Storstrømmen. Thus far, large scale ice models have only been able to capture the distinct flow pattern of NEGIS when using data assimilation techniques such as inverting for the basal friction coefficient (see e.g. horizontal 15 velocity fields in Goelzer et al., 2017). It is assumed that most of the surface velocity can be attributed to basal sliding amplified by basal water instead of ice deformation . This means that the addition of a subglacial hydrology might have the potential to improve the results considerably. While many glaciers in Greenland have regularly draining supraglacial lakes and run-off driving a seasonality of the flow velocities, little is known about the effect at NEGIS (Hill et al., 2017).
Because of this lack of data, to avoid an increased complexity, and to focus on the question if basal melt alone can account for 20 the development of an efficient system, we do not include any seasonal forcing into our experiment.
Our setup includes the major parts of this system. The pressure adjusted basal temperature Θ pmp obtained from PISM (Aschwanden et al., 2016) is utilized to define the modeling region. We assume that for freezing conditions at the base (T pmp < 0.1 K) basal water transport is inhibited and take this as the outline of our model domain. Fig. 8 shows the selected area and PISM basal melt rates used as forcing. 25 For the ice geometry, we use the bed model of Morlighem et al. (2014) interpolated on a 1.2 km grid. Boundary conditions at lateral margins are set to no flux, whereas the termini at grounding lines are defined as Dirichlet boundaries with a prescribed head that implies an effective pressure of zero. This means that the water pressure at the terminus is equal to the hydrostatic water pressure of the ocean assuming floating condition for the ice at the grounding line. Parameters used for this experiment are the same as those for the previous experiments (Table 3) but with K max reduced to 0.3 ms −1 to speed up the computation.  The resulting distributions of effective pressure and conductivity are shown in Figs. 9a and b, respectively. As expected, effective pressure is highest at the ice divide and decreases towards the glacier termini. Conductivity is low for the majority of the study area with the exception of areas in the vicinity of grounding lines and two distinct areas that touch in between 79NG and ZI. The northern area (marked I in Fig. 9b) is located at the northern branch of 79NG and has no direct connection to the snout. The second area (marked II in Fig. 9b) emerges in the transition zone between the southern branch of 79NG and higher effective pressure than 79NG and ZI, which is in accordance with lower observed horizontal velocities for that glacier (Joughin et al., 2010). At the onset of the NEGIS, the effective pressure is high, and no relationship to the flow velocity can be observed.
To further examine the possible influence of our hydrology model to basal sliding, we investigate the impact on the sliding law. We chose to compare our computed N CUAS to the reduced ice overburden pressure, defined in Huybrechts (1990) as 15 N HUY = P i + ρ sw g(z b − z sl ) for z b < z sl and N HUY = P i otherwise.
We show the quotient of H HUY and N CUAS in Fig. 9c. This demonstrates where the application of our hydrology model would increase basal velocities.
In order to demonstrate the effect of the modeled subglacial hydrology system on the NEGIS ice flow, we setup a simple, oneway coupling to an ice flow model. Here, we use the Ice Sheet System Model (ISSM, Larour et al., 2012), an open source finite 20 element flow model appropriate for continental scale and outlet glacier applications (Bondzio et al., 2017;Morlighem et al., 2016). The modeling domain covers the grounded part of the whole NEGIS drainage basin. The ice flow is approximated with the Shallow Ice Approximation (SSA, MacAyeal, 1989;Morland, 1987) within a 2D plan-view model, which is appropriate for fast flowing ice but not for the slow flowing parts in our model area. Since we aim do demonstrate that our addition of the hydrology model improves results in fast flowing areas where the SSA is valid, this is not a huge concern. As we use the SSA 5 we do not perform a thermo-mechanical coupling but prescribe a depth-averaged hardness factor in Glens flow law. Model calculations are performed on an unstructured finite element grid with a high resolution of 1 km in fast flow regions and coarse resolution of 20 km in the interior. The basal drag τ b is written in a Coulomb-like friction law: where v b is the basal velocity vector tangential to the glacier base, N the effective pressure, and k 2 a positive constant. We run  Fig. 10a and c, respectively. Additionally, we show the observed velocities (Fig. 10d, Rignot and Mouginot, 2012) and the PISM surface velocities (Fig. 10b, Aschwanden et al., 15 2016). Note that the latter is a PISM model output on a regular grid interpolated to the unstructured ISSM grid.
Velocities computed with the reduced ice overburden pressure are generally too slow and do not resemble the structure of the fast flowing branches at all. The result from PISM shows distinct branches for the different glaciers, which display a relatively sharp separation from the surrounding area. Note, that PISM also uses a basal hydrology model as described in Bueler and van Pelt (2015). Velocities are slightly lower than observed velocities, especially for Zacharias Isbrae and in the area, where 20 ZI and 79NG are closest. In the upper part towards the ice divide, the ice stream structure is not visible in the velocities. The ISSM model using effective pressure computed by CUAS produces high velocities towards the ocean that closely resemble N .
The transition between the ice streams and the surrounding ice is poorly reproduced. While the stream structure is way too diffused, the velocity magnitude for the glaciers appears reasonable. The inland part is similar to observed velocities but -as in the PISM simulation -the upper part where NEGIS is initiated is not present. The onset of NEGIS is thought to be controlled 25 by high local anomalies in the geothermal flux , which PISM currently does not account for. Higher geothermal flux would lead to more basal melt, hence, water supply in the hydrology model. However, the consequences for the modeled effective pressure would require further experiments which are not in the scope of this paper.
In Tab. 4 we show the root mean square error (l2-norm), Pearson correlation coefficient r 1 and ∆v (l1-norm) between the modeled and observed velocities. 30 We find it impressive that even without extensive tuning, we can considerably improve the velocity field in ISSM by our simple one-way coupling to the hydrology model. However, the results in this section are to be understood not as a thorough study of the NEGIS but as a first application of the model to a real geometry. A complete study requires extended observations in order to determine the optimal model parameters. However, we are confident that our results represent the general aspects of  (Rignot and Mouginot, 2012). Herein RMS denotes the root mean square error or l2-norm, r 2 is the Pearson correlation coefficient and ∆V is the l1-norm. The onset of NEGIS is not well reproduced in the PISM simulation as well as in our ISSM result. Since the ice is slow in the PISM results in that area, basal melt rates are low, and, since we use these as input in our hydrology model, it is expected that our model computes low water pressure here. In our opinion this another point of having a real two-way coupling between the ice model and the basal hydrology model in order to obtain good results. These results could then in turn be used to guide further optimization of the modeling parameters in our hydrology model in the future. 10

Conclusions
We present the first equivalent aquifer layer model for subglacial hydrology that includes the treatment of unconfined water flow. It uses only a single water layer with adaptive conductivity. Since extensive observations of the subglacial system are rare, its relative simplicity and empirical nature can be an advantage.
We find strong model sensitivity to the lower limit of conductivity K min , grid spacing dx, and the parametrization of melt 15 v melt and creep v creep , while the sensitivity to the upper limit of conductivity K max and the confined-unconfined transition parameter d is low. Our model robustly reproduces the seasonal cycle with the development and decline of the effective system over the year. Another positive aspect of our modeling approach was the complete absence of instabilities similar to those arising due to runaway channel enlagement -even at high flow rates. We attribute this fact that our combined confined/unconfined aquifer model quickly transports excess water away from confined aquifer parts, where high water pressure could lead to steep 20 increases in effective conductivity (via negativ v creep term).
In our NEGIS experiments, we find the presence of a partial efficient system for winter conditions. The distribution of effective pressure broadly agrees with observed velocities, while the upstream part is not represented correctly. When coupled to ISSM, our hydrology model notably improves computed velocities.
A number of aspects of the proposed model can be further developed; those include improved parametrizations of several physical mechanisms (e.g. adding feedback between pressure and water supplies), changing the hydraulic conductivity coefficient to a tensor-valued on to better represent the anisotropy of channel networks, and, last but not least, transition to a mixed formulation of the Darcy equation discretized on an unstructured mesh in order to preserve mass conservation and to improve resolution in the areas of interest.

5
Appendix A: Parametrization of evolution of conductivity We use the same parametrization as de Fleurian et al. (2016) detailed here using the notation in Cuffey and Paterson (2010).

Creep term
Nye (1976), found for the closure on channels due to creep that 10 with R c denoting the channel radius and A c the channel area (= πR 2 c ) (notation as in (Cuffey and Paterson, 2010, Eq. 6.15)). Multiplication by 2πρ i R 2 c = 2ρ i A c on both sides, leads to Rewriting the left side to area, using the chain rule ( ∂Ac ∂t = 2π ∂Rc ∂t ) yields (A3) 15

Melt term
Heat produced over ds in unit time is Q w G and pressure melting point effects are ρ w Q w c w B dPi ds , which leads tȯ (Cuffey and Paterson, 2010, Eq. 6.16), whereṀ represents the melt rate (mass per unit length of wall in unit time) and the magnitude of gradient of the hydraulic potential is given by Neglecting the PMP effects we geṫ using Q w = qb (confined case, unconfined would be Q w = q(h − z b )) and q = K∇(h) (ommiting the minus, because we need the magnitude here) this is which can be rewritten tȯ Opening and closure The conduit expands when there is more melt than ice inflow due to creep (Cuffey and Paterson, 2010, Eq. 6.42). InsertingṀ from Eq. A8 and dividing by ρ i results in (A10)

Change of area to conductivity
This is a purely geometrical argument, since we are interested in the transmissivity T = Kb (T = Kh for unconfined) and therefore, a change in K is equivalent to a change in the aquifer thickness b. 10 Hence, we add a roughness term r to be able to investigate the sensitivity of our model to the melt term and write and considering unconfined case as well, it is Our reasoning behind scaling K and not b are unintended side effects of changes in b on storage and the behaviour when 15 the aquifer becomes unconfined. A thicker aquifer in areas which represent channels would imply larger storage and also less transmissivity when the area eventually becomes unconfined (because the unconfined case triggers at a higher head already, limiting transmissivity). However, we do not expect the storage to increase, because channels usually do not store a lot of water and pressure changes travel faster than in the inefficient system. Also, the effective system should still have higher transmissivity when the water pressure drops (as long as it does not empty completely).

Appendix B: Discretization
We discretize the transient flow equation (Eq. (3)) on an equidistant rectangular grid using an explicit forward in time central in space (FTCS) finite-difference approximation. For sake of completeness, we give the equations for a non-equidistant grid here.
For the spatial discretization, we use a second-order central difference scheme (e.g., Ferziger and Perić, 2002) leading to the spatial discretization operator for the head L h : where half-grid values of T denote harmonic rather than arithmetic averages computed using Eq. (4), where denote central, forward, and backward differences, respectively. Re-writing this more compactly in compass notation We use the explicit Euler method for the time discretization, where the next time step m + 1 is computed from the previous time step m (∆t = t m+1 − t m ) (using a somewhat sloppy notation) Conductivity is updated in each time step by Eq. (9): where we use a combined forward-backward-difference scheme for the discretization of (∇h) 2 in Eq. (10): Compared to central differences, this stencil is more robust at nodes with large heads caused by moulins.
The time step is chosen sufficiently small that the discretization error is dominated by the spatial discretization (the Courant-Friedrichs-Le condition is always satisfied) Additionally, we check that the time step is small enough for the unconfined component of the scheme to become active by restarting the time step with a decreased ∆t if at any point h < z b .
All variables are co-located on the same grid, but the conductivity K is evaluated at the midpoints between two grid cells 5 using the harmonic mean due to its better representation of conductivity jumps (e.g. at no-flow boundaries).
A disadvantage of this discrete formulation is that it is not mass-conservative (see, e.g. Celia et al. (1990)). The solution to this is to use a mixed formulation for Darcy flow in which also the Darcy velocity is solved for. However, in our application, the resulting error is very small, and we plan to implement the mixed formulation approach in future work.
Competing interests. The authors declare that they have no conflict of interest.