Hydrology and runoff routing of glacierized drainage basins in the Kongsfjord area, northwest Svalbard

Freshwater discharge from tidewater glaciers modulates fjord circulation and impacts fjord ecosystems. There can be significant delays between meltwater production at the glacier surface, and discharge into the fjord. Here, we present a hydrological analysis of the tidewater glaciers around Kongsfjorden, northwest of Svalbard, examining the pathways of glacier surface melt to the glacier fronts. To simulate discharge hydrographs at the outlets of the major drainage basins in the Kongsfjord area we use 1) a simple, heuristic routing model and 2) the physically-based model HydroFlow to route runoff 5 derived from a coupled surface energy balance – snow model. Plume observations at one of the tidewater glacier outlets and measurements of proglacial discharge of a land-terminating glacier are used for model calibration. Our analysis suggests that the local subglacial topography diverts a substantial amount of water from the drainage area of the glacier Kongsbreen to the neighboring glacier Kronebreen, across the border of their surface catchments. This is supported by the relative sizes of the plumes observed at the respective glacier fronts. Runoff from the glaciers on the south side of the fjord is one order magnitude 10 lower than runoff from the glaciers on the east and north sides of the fjord, reflecting differences in the size of the glaciers. We derive discharge hydrographs at all the major outlets of Kongsfjord basin, presenting here a detailed analysis of two of the glacierized basins. The average annual discharge period from the tidewater glaciers due to surface runoff was 105±10 days. The largest discharge comes from Kronebreen, which is equivalent to around 40% of the total freshwater flux to the fjord.


Surface and bed digital elevation models
We use a 5-m gridded digital elevation model (DEM) of the glacier surface, constructed from aerial photographs taken in 2009 and 2010 (Norwegian Polar Institute, 2014), which is resampled onto a 250-m grid (Fig. A2a). The basal topography of the tidewater glaciers in the Kongsfjord basin has been mapped previously with airborne and ground-based ice-penetrating radar (Lindbäck et al., 2018) to derive a 150-m bed DEM, which is resampled to the same 250-m grid as the surface DEM ( Fig.   100 A2b).

Plume data
Plume observations were made with a time-lapse camera over the period 2014-2016. The camera was installed on Collethøgda, a mountain between the termini of Kronebreen and Kongsbreen, and took pictures at every hour. The extent of the plume at the fjord surface was digitized manually from cloud-free images, as described in How et al. (2017). Subsequently, plume extents 105 were georectified with the PyTrx toolset as outlined in How et al. (2020), using a set of ground control points, camera model, and the surface DEM described previously (Norwegian Polar Institute, 2014).

Bayelva discharge data
The total area of the Bayelva basin is 33 km 2 , of which 17 km 2 is glacierized. In early summer, discharge is from snowmelt runoff of the non-glacierized area, and runoff of snow and ice melt from glaciers. In late summer, rainfall and runoff from 110 the glaciers contribute to water flow in the river. Discharge measurements at the Bayelva station ( Fig. 1) are conducted automatically; the water level in a concrete-floored weir is measured at hourly intervals by a pressure transducer, and a float and wire system (Killingtveit et al., 2003). The system is calibrated periodically to derive a rating curve that converts water 4 https://doi. org/10.5194/tc-2020-197 Preprint. Discussion started: 21 September 2020 c Author(s) 2020. CC BY 4.0 License. level to discharge. Over the period 1990-2000, only daily data are available, whereas over the period 2000-2015 data were stored hourly. Because the monitoring is unattended, the discharge data have periods with erroneous readings, caused by ice or 115 sediment build-up at the sensor; however, the timing of discharge events (i.e. the peaks) is in general not affected.

Flow paths and drainage basins
We use flow algorithms in the MATLAB package Topotoolbox (Schwanghart and Scherler, 2014) to determine flow direction and drainage basins. The land-terminating glaciers are constrained within relatively narrow valleys, and observations show 120 that runoff routing from these glaciers is controlled predominantly by surface topography alone for most of the glacier area, therefore, we use the surface DEM to evaluate flow directions and drainage basins. The tidewater glaciers are thicker (Lindbäck et al., 2018) and runoff routing is influenced by both surface and subglacial topography (Vallot et al., 2017).
Hydraulic potential is the sum of the basal water pressure and elevation potential; the former is approximated as a fraction, k, 125 of the ice overburden pressure via: where ρ i is the density of ice (916 kg-m −3 ), ρ w is the density of water at 0 • C (1000 kg-m −3 ), g is the gravitational acceleration in m-s −2 , H is the ice thickness in m, and Z i is the ice surface elevation in m (Shreve, 1972). The hydraulic head, h (in metres), can be expressed as follows: When k equals 1, the subglacial water pressure is equivalent to the ice overburden pressure, a situation we may expect in winter, in the absence of low-pressure channels (Zwally et al., 2002). If k equals 0, the subglacial water is at atmospheric pressure and water moves according to the bed topography alone. To calculate hydraulic potential, we assume a spatially and temporally constant value of k, and route water accordingly. 135 We use the 250-m grids of surface and bed elevations to derive the hydraulic potential, and calculate subglacial drainage networks with Topotoolbox (Schwanghart and Scherler, 2014). The model calculates the slope of each grid cell, with respect to the adjacent eight grid cells, and water is assumed to flow along the steepest direction. Depressions in the hydraulic potential surface are filled to remove sinks. Grid cells that are not receiving any flow from upstream lie on the boundary of the basin.
Each basin has a single outlet where all the grid cells of the basin drain, such that these grid cells demarcate the respective 140 drainage basin area. We investigate the effect of changes in water flow, and thereby the catchment areas, by varying k values between 0 and 1.
The hydraulic potential may depend on the resolution of the DEMs, where uncertainties in hydraulic potential increase with the coarser resolution. Therefore, we further calculated hydraulic potential for different DEM resolutions of 100-m and 150-m.
To check the robustness of the analyses, we conducted 10000 Monte-Carlo simulations in each cases of hydraulic potential calculation; i) with spatially distributed random k values and ii) with randomly perturbed ice thickness.

Runoff routing models
We use two approaches: the first uses a simple routing model where water moves at uniform speed to calculate the discharge hydrograph at the outlet points of each catchment. The second routing approach uses a topographically-controlled linear reservoir model, HydroFlow, to calculate discharge hydrographs at the outlet points of the studied glaciers.

Simple routing model
The main idea behind the simple routing model is that distance and water wave speed are the first-order factors affecting runoff routing (Kohler, 1995). While the wave speed is known to be a function of many factors, including surface slope, water storage, presence of crevasses and moulins, and seasonal changes in supraglacial, englacial, subglacial channel dimensions and roughness, we here apply a uniform and constant water movement speed to calculate the discharge hydrograph at the outlet 155 points of all the drainage basins Slater et al., 2017).
Time-varying runoff from the energy balance -snow model (Sect. 3.2) is used as the input to the routing model using 250 m gridded DEM of the surface and bed elevation. We use the hydraulic potential (calculated in Sect, 4.1) to derive flow paths and the distance of each grid cell to its basin outlet point. We assign a uniform water speed to calculate a time delay associated with each grid cell, and sum the delayed runoff from all the grid points in an individual glacier basin to derive the discharge 160 hydrograph at its outlet point. We test different wave speeds (0 -1 m-s −1 ) with intervals of 0.1 m-s −1 to consider the effect of the slower or faster flow of water, and use observational plume and discharge data to determine an optimal wave speed.
The primary factor controlling plume generation is runoff from glaciers. Plume area extent depends on several factors (such as wind speed and direction, terminus depth, discharge, and fjord stratification and circulation), and there does not appear to be a simple correlation between the discharge and plume area on a seasonal timescale (How et al., 2017). Our main assumption here 165 is that the timing of maximum plume extent coincides with maximum discharge at the front. We derive discharge hydrographs at Kronebreen outlet for different water speeds, and compare the peaks in the delayed runoff to the peaks in plume extent. In the absence of any direct discharge measurements at the tidewater glacier front, the water speed for the Kronebreen glacier runoff route is calibrated by computing the normalized cross-correlation between the high-frequency components of the modelled discharge and plume area data (Fig. B1, discussed in detail in the Appendix B).

170
We also compare the discharge hydrograph at Bayelva to that modelled in the Brøggerbreen glacier catchment. A single general value for the water speed of the simple routing model is tuned for each of the cases; subglacial water routing of tidewater glaciers and supraglacial water routing for land-terminating glaciers. For Bayelva, the Nash-Sutcliffe coefficient (NSC) is calculated over the year 2000-2010, quantifying the agreement between the modelled and observed hydrographs. We calculate the time difference between the peaks of no-delay and delayed discharge, thereby referred to as delay-time, which is 175 governed by the assumed value of the water speed and the distance.

HydroFlow model
The second routing model, HydroFlow, uses a linear-reservoir approach to calculate discharge hydrographs at the outlet points of the glaciers (Liston and Mernild, 2012). HydroFlow has been applied previously in different glacierized and river catchments (Mernild et al., 2015(Mernild et al., , 2017. It was developed by Liston and Mernild (2012) and first tested in the Mittivakat glacier catchment 180 in southeast Greenland. The model is described in detail by Liston and Mernild (2012), and a brief summary is given below.
HydroFlow is a gridded linear reservoir model, where each grid cell acts as a linear reservoir, transferring water to the steepest adjacent cell. The routing is conducted in four routines. First, the model calculates the topographically controlled flow networks from the surface DEM. Second, it calculates the individual watersheds inside the domain of interest. In the third part, it assumes that there are two different components associated with water transport: a slow-response and a fast-response system.

185
The slow time-scale accounts for distributed runoff through the grid cell matrix, for example flow within the snow, ice, and soil, whereas the fast time-scale represents water transport through the channelized system and includes supra-glacial, sub-glacial, or en-glacial flow. Here, we used hydraulic potential instead of surface DEM to make water flow according to the hydraulic potential gradient for tidewater glaciers. All parameters of the model in this study are adopted from Liston and Mernild (2012), except the timescale α associated with the fast-flow, which we calibrate, as described in detail below.

190
HydroFlow is a linear-reservoir model and does not consider explicitly subglacial hydrology in its flow routine. For landterminating glaciers of this region, supraglacial channels play a major role to transport water, but for tidewater glaciers, subglacial hydrology is important. We calibrated the fast-time scale coefficient α of the HydroFlow model to derive discharge hydrographs for Bayelva and Kronebreen, following a similar procedure for the simple routing model (sect. 4.2.1); we compute NSC and the normalized cross-correlation between modelled and measured discharge for Bayelva, and for high-frequency 195 components of the modelled discharge and observed plume for Kronebreen, respectively. Two separate coefficients were derived for the tidewater and land-terminating glaciers.

Subglacial hydrology and drainage basins
Based on the surface DEM, there are 114 basins draining into Kongsfjorden, of which 15 are either completely or partially 200 glacierized, and are named after the respective glaciers (Fig. 1). From the hydraulic potential at the base of the tidewater glaciers, there are five major subglacial catchments with an area of 50 km 2 or larger (Table A1). We vary k values between 0 and 1 to find that the subglacial drainage delineation deviates from the outlines of the surface catchments at the ablation zone of the glacier, implying subglacial transfer of water across the boundaries of surface catchments, a phenomenon referred to as water piracy (Anandakrishnan and Alley, 1997;Lindbäck et al., 2015). Water piracy is apparent between the two largest adjacent 205 ice fields of the basin, Holtedahlfonna and Isachsenfonna, which would drain to Kronebreen and Kongsbreen, respectively, according to the surface flow. We find that for a k value between 0.5 and 1, a substantial part of Isachsenfonna drains to Holtedahlfonna subglacially, whereas the reverse happens for k values between 0 and 0.1 (Fig. A3). We also found that for k values between 0 and 0.2, Kronebreen drains through an outlet at the southern side of the glacier which merges with the Kongsvegen outlet and drains as a single outlet point (Fig. A4, A5). Monte-Carlo simulations of hydraulic potential with 210 spatially distributed random k values show around 94% of water piracy cases are from Isachsenfonna to Holtedahlfonna, whereas only 5% of cases are from Holtedahlfonna to Isachsenfonna, and below 1% of cases occur with no water piracy between these two glaciers (discussed in detail in Appendix Sect. A2). Fig. 2 shows the changes in drainage catchment for surface DEM and for hydraulic potential corresponding to k = 0.1, 0.5, and 0.9, representing substantial changes in drainage basins. For k values between 0.5 and 1, a substantial amount of water from Isachsenfonna flows to Kronebreen. Although the 215 catchment area of Kongsbreen decreases by 75% for this water piracy, the total area is still 82 km 2 , the fourth largest catchment in this area. This water piracy affects the water influx to the fjord from Holtedahlfonna and Isachsenfonna in late summer only, since the switching between Holtedahlfonna and Isachsenfonna takes place in a small area situated in the upper part of these two glaciers, where meltwater is only produced in late summer. Because of this, water influx to the fjord through the outlets of these two catchments will differ substantially only in peak summer (July-August), when runoff occurs from the upper part of

Discharge hydrograph 225
Normalized cross-correlation is calculated for different water speeds of the simple routing model and for different α values of the HydroFlow model, comparing the discharge and plume area data (Fig. B2). From the visual inspection of Normalized cross-correlation, we found the best fits to the data for water speeds of 0.6 ms −1 , and for α values of 0.2, which are used for the final model run. The same water speed of 0.6 ms −1 is also optimal for Bayelva, yielding the highest Nash-Sutcliffe coefficient in all years (Fig. B3a). However, the best α values, as determined by the highest Nash-Sutcliffe coefficients are 0.15 and 0.2 230 (Fig. B3b). We observe, however, that α = 0.15 results in a smoothed hydrograph, whereas α = 0.2 gives sharper peaks which better resemble the measured discharge peak. We therefore choose α = 0.2 as the optimal value.    Kronebreen, comprising 39±1% of the total discharge to the fjord. Annual runoff from tidewater glaciers is about one order of magnitude greater than that from land-terminating glaciers, and is corresponding to their area proportion. discharge being around five times higher than Kronebreen discharge. The potential subglacial water capture from Isachsenfonna to Holtedahlfonna is supported qualitatively by satellite and terrestrial time-lapse imagery (Fig. A6), which consistently shows a significant plume in front of Kronebreen (How et al., 2017), but only a small plume at Kongsbreen (Schild et al., 2018). In addition, mammals and birds are observed to forage preferentially at the front of Kronebreen (Lydersen et al., 2014;Urbanski et al., 2017), an indication of vigorous upwelling induced by subglacial discharge. To further examine the question of subglacial water piracy, simulations using a process-resolving model of subglacial hydrology are required which more realistically elucidates spatio-temporal evolution of subglacial water pressure and associated hydraulic potential. Frequent and continuous monitoring of plume area extent at those two outlet glaciers is essential to infer more information about time-evolving subglacial drainage as long as no other reliable technique to monitor plume discharge is

Discharge hydrograph
The subglacial hydrology of glaciers and ice sheets is difficult to model due to the complexity of the processes, and the relative lack of data at the bed. A number of hydrological models with different degrees of complexity have been developed (Hewitt, 2011;Werder et al., 2013;Bueler and Van Pelt, 2015;Flowers, 2015). It is an ongoing debate among hydrological modelers whether a higher degree of complexity improves model accuracy, especially in predicting freshwater discharge at tidewater 300 glacier outlets (Li et al., 2015;De Fleurian et al., 2018). Here we are mainly interested in calculating discharge hydrographs at the tidewater glacier outlets on sub-daily timescales. We argue that the simple routing model allows adequate quantification of the discharge hydrographs at outlet points around the Kongsfjord basin, given that the main intent is to derive time delays at the outlets. The discharge hydrographs derived from the simple routing scheme compares well to those from the more physically realistic linear-reservoir model HydroFlow. We used plume data as a proxy to discharge to calibrate the parameter of the simple 305 routing model. We found a little difference in the time period of discharge and plume data for Kronebreen for one year. In 2014, some plume activity was visible in late September, even though the simulated runoff was close to zero. This late plume activity could be caused by runoff from basal melting or from a late rain event occurring at lower elevations, not captured by the model (Pramanik et al., 2018(Pramanik et al., , 2019. Kronebreen is polythermal and fast-flowing glaciers, thus there can be substantial basal melting in early summer as well as in late autumn. The shallow water depth of Kronebreen front makes the runoff, even with low magnitude, appear at the surface. Although we assign an optimum water speed for the routing model, we propose that the water wave speed will play a crucial role when discharge is calculated over a sub-hourly timescale. We did not try to find an accurate water speed value as our model presents the complex physical process through simple parameterization. Our discharge hydrographs are relatively less sensitive to the choice of water speed due to the relatively coarse temporal resolution of the modelled discharge data. We observe that 315 a simple routing model performs equally well with a physically-based linear-reservoir model when discharge is simulated over sub-daily timescale. However, we suspect that a major difference between a single parameter simple routing model and a complex model will be apparent when simulation is done over sub-hourly timescale, as complex models are expected to capture details of the hydrology and water routing. However, a complex model may even fail to produce an accurate discharge hydrograph due to a lack of temporal observational data from the subglacial environment. In all these models, uncertainty comes 320 in simulating discharge timing, but the magnitude of discharge remains unaffected. Nonetheless, with a lack of observational data, a simple model serves the basic purposes of biogeochemistry and oceanographic studies (Sundfjord et al., 2017;Everett et al., 2018;Schild et al., 2018;Halbach et al., 2019).
The five tidewater glaciers of the Kongsfjord basin contribute most of the freshwater flux to the fjord. This has significance as mixing of freshwater from glaciers with the fjord seawater influences circulation (Sundfjord et al., 2017), and enhances 325 primary and secondary production, thus playing an important role in the fjord ecosystem (Lydersen et al., 2014).
Here, we used a time-independent modelling approach to route meltwater to the outlet. We assume stationary and uniform values for hydraulic potential, water travel speeds and fast-response timescale for the entire season. However, in nature, these values would depend on several factors e.g., evolution of channels, surface slope, bedrock property, and vary temporally.
Further uncertainty comes from calibration using plume area observations as a proxy for subglacial discharge. Plume extent is 330 not a direct signal of outflow, and the routing model could further be improved by incorporating more time-dependent physical processes. Future studies could combine detailed hydrology and routing in a model framework to get a more robust estimate of discharge at the glacier outlet points.

Conclusions
Freshwater influx from glaciers substantially impacts fjord circulation and fjord ecosystem. Glacier hydrology and water rout-335 ing play the major role in controlling discharge hydrographs at fjord inlets. We analysed the subglacial hydrology and water routing of the entire glacierized area of the Kongsfjord basin in order to simulate discharge hydrograph of freshwater influx at different inlets points of the fjord. Subglacial hydrology of this region is poorly understood, and here, using steady-state hydraulic potential analyses, we hypothesize a structure of the subglacial environment of the basin's tidewater glaciers. We suggest that there is a higher possibility of subglacial water piracy between Isachsenfonna and Holtedahlfonna, where the latter 340 receives substantial water subglacially from the former. Our hypothesis is supported by the relative size of the plumes at the two tidewater glacier fronts. Furthermore, we implement a simple routing model to derive discharge at the different drainage catchments around Kongsfjord. The discharge hydrograph of simple routing model is compared with the ones derived from Hy-droFlow. We conclude that, with lack of observation data of subglacial conditions, the simple routing model performs equally well with HydroFlow in simulating discharge hydrographs over sub-daily timescales.

345
Appendix A: Hydrology

A1 Uncertainty
To determine the uncertainty in the hydraulic head, we estimated the uncertainty σ h using standard analytical error propagation method. The uncertainty in σ h is calculated as Where σ h is standard deviation in ice thickness H, σ Zi is the standard deviation in surface DEM Z i , σ ρi is the standard deviation in the density of ice ρ i , σ HZi is the covariance of H and Z i . The uncertainty in ice thickness DEM was calculated to be ±24 m (Lindbäck et al., 2018). The surface elevation dataset has a standard uncertainty of 2-5 m (Norwegian Polar Institute, 2014). The uncertainty of the density of water is very small, hence neglected. The combined estimated uncertainty in the hydraulic head is ±22 m.

A2 Sensitivity (Monte Carlo)
To check the robustness of hydraulic potential analyses, we conducted sensitivity tests with different DEM resolution and monte-carlo simulations with randomly distibuted k values and with random Ice thickness perturbations.

A2.1 Different DEM resolution
We conducted the hydraulic potential analysis for bed and surface DEM of 150-m and 100-m resolution to check the depen-360 dency of drainage delineation and water piracy results on the DEM resolution, and applicability of the analyses in routing model. We found that the drainage delineations are little sensitive to resolutions, however, the water piracy between Holtedahlfonna and Isachsenfonna is aptly represented in all the DEMs (Fig. A3-A5).

A2.2 Random spatially distributed k values
We did 10000 Monte-Carlo simulations for hydraulic potential with spatially distributed random k values ranging between 0 365 to 1. Thereby, we investigated the changes in subglacial drainage delineations. We found 94.36% cases where Isachsenfonaa is draining to Kronebreen, 0.38% cases with no subglacial water piracy and 5.25% cases where Holtedahlfonna is draining to Kongsbreen-North. Between Kronebreen and Kongsvegen, we found 29.72% cases where Kongsvegen and Kronebreen drain separately whereas rest 70.28% cases show a single outlet for Kronebreen and Kongsvegen.
We also did 10000 monte-carlo simulations for hydraulic potential with random DEM perturbations. We perturbed the Ice thickness by adding random noise within the error of Ice thickness measurement, which is ±24 m (Lindbäck et al., 2018). We found that in all the runs, Isachsenfonna drains to Holtedahlfonna, thus shows that the water piracy between Holtedahlfonna and Isachsenfonna is not sensitive to the DEM perturbation within the error limit of bed DEM.

B1 Calibration of parameters of Routing model
We use a median filter of 60-hours to filter out the low-frequency part of the signals from both modelled discharge and plume data (Fig. B1). Thereafter, we calculate normalized cross-correlations between the high-frequency part of the plume and modelled discharge for different water speeds. We found the least lag between modelled discharge and plume data for the water speed value of 0.6 ms −1 , and use this value for the final model run. We choose a cut-off frequency to filter out low-380 frequency components and compare the high frequency of plume with the high frequency of discharge data. We did not find good correlation of plume and discharge data for all cut-off frequency. We used different cut off frequencies and calculated normalized cross-correlation. A good correlation is observed for the cut-off frequency of 59 hours and above. Therefore, we use that as the final cut-off frequency to extract high-frequency signal to calibrate the water speed of the routing model.