Estimating meltwater retention and associated nitrate redistribution during snowmelt in an Arctic tundra landscape

Nitrogen availability in Arctic ecosystems is a key driver for biological activity, including plant, growth and thereby directly linked to the greening of the Arctic. Here, we model the redistribution of meltwater following spring snowmelt as well as the accumulation of meltwater and dissolved nitrate at landscape scale. By combining snow mapping with unmanned aerial systems, snow chemistry, and hydrological modelling, we argue that the majority of nitrate in the snowpack is flushed out of the landscape due to the limited storage capacity of meltwater in the early growing season frozen soil. We illustrate how landscape micro-topography is a crucial parameter to quantify storage capacity of meltwater at landscape scale and thereby the associated pool of soluble compounds such as nitrate. This pool will be available for plants and may be important for plant diversity and growth rates in the wettest part of the landscape. This study illustrates that the evenly distributed nitrate input during the Arctic winter may be redistributed during the initial snowmelt and lead to marked differences in biologically available nitrate at the onset of the growing season, but also that the majority of deposited nitrate in snow is lost from the terrestrial to the aquatic environment during snowmelt.


Introduction
There is a current focus on understanding the varied responses of Arctic terrestrial ecosystems to climatic changes, specifically the net effect of an increasing (Forkel et al 2016) or a decreasing carbon sink strength (Nauta et al 2015, Dahl et al 2017. Nitrogen (N) in a plant available form is a key component for these responses since N is considered a limiting factor for plant growth in the Arctic (Shaver et al 2006). Future warming and changes in precipitation patterns are likely to catalyze mineralization of stored N, but also potentially increase other sources such as atmospheric deposition, which could lead to alterations in Arctic vegetation dynamics (Rousk et al 2017a, Bokhorst et al 2018, D'Imperio et al 2018. Main sources of plantavailable N at catchment scale in most pristine Arctic ecosystems include redistributed internal N from mineralization of soil organic matter and inputs from fauna, and external input from atmospheric deposition and N 2 fixation (Skrzypek et al 2015). Mineralization is limited in poorly aerated soils and atmospheric N 2 fixation may therefore play an important role, particularly later in the season when fixation rates increase with higher summer temperature (Rousk et al 2017b). It also means that early season N availability in landscape depressions may rely on external N sources in contrast to aerated soils relying more on internal N (Semenchuk et al 2015).
A marked part of the atmospheric N is deposited into a winter snowpack as plant-available nitrate (NO 3 − ), particularly at higher latitudes where a substantial part of the precipitation (often >50%) falls as snow (Przybylak et al 2003). In addition, a higher NO 3 − concentration has been reported in precipitation during winter than during summer in the whole Arctic (Hole et al 2009), and locally in e.g. Svalbard (Kühnel et al 2011). This increase is partly related to anthropogenic sources, which is evident in the Greenland Inland Ice Sheet snowpack (Hastings et al 2009) and in cores from Arctic lakes (Holmgren et al 2010, Holtgrieve et al 2011. When focusing on the spatial distribution of NO 3 − in the snowpack in Greenland, concentrations in the Inland Ice Sheet seems to increase towards the southern and western regions due to an enhanced effect of anthropogenic sources, with the highest degree of variation found in low accumulation zones (Burkhart et al 2009). When observed at a smaller scale as coastal-inland gradients, Curtis et al (2018) reports decreasing concentrations of deposited NO 3 − in snow towards the coast, yet a higher net deposition than inland due to higher precipitation. The winter snowpack melting and the subsequent runoff is a decisive annual hydrologic event in Arctic terrestrial ecosystems, during which soluble ions in the meltwater are redistributed over time and space. Studies of snowmelt and subsequent runoff at catchment scale traditionally combine in situ observations of snow parameters or snow distribution model outputs with energybased snowmelt models (Kane et al 1997, Bruland et al 2001, Liston and Elder 2006 and raster-or tile-based hydrologic routing schemes (Zhang et al 2000, Dornes et al 2008. Two advantages of such approaches are estimates of snowmelt and runoff over time, and the possibility to add levels of soil infiltration. Modelling the winter snow, including redistribution, requires sophisticated coupling of meteorology, topography and energy-balance data. Moreover, it is demanding in processing power, particularly when using high spatial resolution, and would have the inherent uncertainty of pure model studies. And while the simplest of snowmelt models can be forced with temperature data only, most of them require complete energy-balance datasets (Pohl et al 2005). Similarly, runoff-routing with infiltration capability require processing power and relies on detailed mapping of soil characteristics over time, particularly soil moisture and structure; datasets which are often scarce in remote areas. This calls for an efficient way to actually measure snow depths and an optimized method to quantify runoff routing.
When linking meltwater runoff with nutrient reallocation, an important component is the ionic pulse occurring during snowmelt (Rasch et al, 2000, Elberling et al 2007, Boucher and Carey 2010, which reflects that most water-soluble ions in a snowpack are flushed out in the early phase of snowmelt leaving the remaining snow drifts with almost distilled water and solid particles behind. The specific ratio of percent snowmelt to flushed ions will be site-and species-specific. Due to this ionic pulse, the movement of early meltwater and water storage capacity in the landscape during early spring with completely frozen tundra may be critical for the overall redistribution of winterdeposited N species. Previous studies from terrestrial Arctic sites have successfully linked the loss of N at catchment scale to snow meltwater runoff over time (e.g. Jones et al 2005, Boucher andCarey 2010). While these studies focus on the temporal component of loss of soluble ions such as NO 3 − , a combined spatiotemporal approach describing NO 3 − redistribution is for instance given by Nowak and Hodson (2015). However, Nowak and Hodson (2015) do not link the spatial distributions to specific measurements of deposited NO 3 − in snow or the landscape water retention capacity. To increase our understanding of this mechanism, access to high spatial resolution data of snow distribution is one key aspect, alongside detailed topographic representations to evaluate meltwater runoff patterns, and ground-based measurements of snow chemistry. Recent development in lightweight drones and increased computer processing power has revolutionized the use of unmanned aerial systems (UAS) in biogeosciences within the last decade. When focusing on small drones (<5 kg), aircrafts with high-quality consumer cameras can acquire overlapping aerial photos allowing the use of structure-from-motion techniques to photogrammetrically derive 3D point clouds of the photographed surfaces (James andRobson 2012, Lucieer et al 2014a). Subsequent digital surface models (DSM) and orthophotos are then created by interpolating planes between points, thereby reproducing topography typically at 1-20 cm spatial resolution. In remote areas such as the Polar regions, this data source is highly versatile for describing the heterogeneity and topography of natural landscapes and has for instance been used to map micro-topography of specific vegetation types (Lucieer et al 2014b) and evaluate surface heterogeneity at landscape scale (D'imperio et al 2017). Although such high-resolution DSMs are clearly advantageous to quantify runoff routing in permafrost landscapes where vegetation is absent or low, and would be a novel use of the technique, one limitation with traditional raster-based runoff approaches is that they are inefficient on high-resolution datasets.
When combined with cm-precision ground control points (GCP), the surface models have moreover been found to offer sufficient spatial accuracy to map snow depths using subtractions of co-located snowcovered and snow-free surfaces. Vander Jagt et al (2015) showed that this approach can provide snow depths with <10 cm error in a naturally vegetated local (<1 ha) alpine area. Harder et al (2016) reports a similar <10 cm accuracy in an alpine site, as well as for a prairie site, however, the snow depth accuracy was lower at sites with taller vegetation. This observation is supported by Bühler et al (2017) who achieve accuracies ∼7 cm in sparsely vegetated areas, and up to ∼30 cm accuracy in areas with vegetation of that particular average height. Consequently, auxiliary data on surface classes and vegetation heights are important for the precision of snow depth estimates.
In summary, the spatial fate and retention of snowpack NO 3 − in Arctic landscapes during snowmelt is uncertain and dependent on micro-topography.
High spatial resolution surface models and snow measurements are needed to capture this topography, but can be achieved with UAS. The subsequent runoff modelling and landscape retention capacity requires new efficient network-based models due to the high spatial resolution. Here we present a new integrated approach to assess the landscape retention capacity of meltwater and dissolved NO 3 − within a catchment.
This includes (1) assessing the spatially distributed concentrations of NO 3 − in the snowpack, (2) quantifying 3D snow depths using UAS, and (3) quantifying the run-off pattern and retention-capacity of topographic dips, using a high-resolution digital surface model coupled with a high-efficiency surface runoff model. Our aims are to show the ability of this approach as a precise and efficient screening tool, and show for an Arctic tundra site how much NO 3 − is likely lost from the melting snowpack.

Study site
The Blaesedalen valley at 69.26°N 53.47°W located on the southern part of Disko Island in West-Greenland has been ice-free since the last glaciation approx. 10 000 years ago and represents today a heterogeneous mesic tundra landscape. The southern part of the valley is elevated approx. 120 m above sea level, and the studied area is located 500 m north of the coast. Annual mean temperature is −3.0°C, the warmest month being July at 7.9°C and coldest being February at −14.0°C (data from 1991 to 2014; Hollesen et al 2015), causing the area to be underlain with discontinuous permafrost (Westergaard-Nielsen et al 2018), with a mean annual soil temperature in the upper 5 cm of −0.9°C (data from 1991 to 2004; Hansen et al 2006). Precipitation data from the area are very scarce, but the annual average has been estimated to ∼400 mm (Hansen et al 2006). The ratios of land surface classes in Blaesedalen are considered representative for greater parts of Greenland  and on broader scale substantial parts of the circumpolar Arctic with dominance of dwarf-shrub tundra.

Unmanned aerial system data and ground control points
The UAS-based datasets were collected in the snowfree and snow-covered period, respectively. The snowfree data was collected in 2014 and repeated in 2017 (all in July/August), and the snow-covered data during peak snow depths in 2018 (April 2nd). All flight routes were planned in Mission Planner (ArduPilot Dev. Team) with parallel flight lines and 75% image overlaps. All flights were conducted under clear-sky conditions. Table 1 lists the details of the equipment used and derived products.
GCPs were measured to correct the geometry of the DSMs using a differential global positioning system (DGPS, Two Trimble R8 receivers in a post-processing kinematic setup), and distinctive locations estimated to be visible on both the snow-free and snow-covered orthomosaics were selected (abrasion plateaus with distinct corners, centres of boulders, masts etc). A total of 14 GCPs were recorded in the snow-free landscape and initially 10 in the snow-covered landscape. On both occasions, the GCPs were distributed evenly across the study area. DGPS data was postprocessed in Trimble Business Center 3.82, with an average horizontal and vertical precision of 3.4 cm and 6.1 cm, respectively. We added 14 GCPs (giving a total of 24 GCPs to correct the snow-covered DSM) derived from the georeferenced snow-free model from distinct surface features that were also visible on the snow-covered model (barren rocks, masts etc), in order to increase the relative precision between the two surface models.
The imagery from the UAS was postprocessed with Agisoft Photoscan Pro 1.4.5 using a structurefrom-motion workflow with GCPs. The workflow can be summarized as: (1) image alignment and sparse point cloud generation (key points), (2) optimization of camera alignment using 14 ground control points (snow-free scene) and 24 ground control points (snow-covered scene), (3) dense point cloud generation, (4) DSM generation based on the dense point cloud and (5) orthomosaic generation based on the DSM.
The images were resampled to half the original pixel resolution during the processing using Photoscan Pro. Trials using full image resolution did not improve the results in this study but increase the processing time substantially. Both the snow-free and snow-covered DSM represented a full catchment, delineated topographically towards north, west, and east. We then subtracted the snow-free DSM from the colocated snow-covered DSM, resulting in a snow depth data layer covering 1.347 km 2 .

Vegetation classification
Using structure-from-motion on image data allows for developing a DSM which will not represent the true terrain, unless the barren surface is sufficiently observable in the images. We therefore classified the surface cover types in the area to allow a subtraction of average vegetation heights from the snow-free surface model. In principle, this will result in an actual terrain model. The classification was based on a pan-chromatically sharpened WorldView-2 scene covering the entire study region, acquired on the 12th of June 2012. The scene consists of 5 spectral bands (WorldView-2 band) 2-6 at 0.5 m spatial resolution and 1 pan-chromatic band (for details see www.digitalglobe.com). We used a Maximum Likelihood classifier based on statistics from 11732 training pixels distributed across 7 surface  3 for accuracy assessment).

Validation data of snow properties
Before the UAS survey, we sampled 22 snow cores of the entire snow column along a snow depth gradient using a Snow-Hydro plexiglass probe (www. snowhydro.com) and DGPS measurements for accurate positioning of the sample locations. We wore gloves to avoid contamination of the snow samples, and cleaned the probe with snow between each sample. After the UAS survey, we took 98 snow depth probe measurements using a MagnaProbe (Sturm and Holmgren 2018). The latter are limited by a maximum probe depth of 120 cm, and the inherent spatial uncertainty of an ordinary GPS (specific precision not reported by the MagnaProbe). Collected snow samples from the probe with a known volume were stored in plastic bags until melted and then weighed to derive snow water equivalents (SWE). Gloves and plastic bags were tested to be N-free. We took subsamples from each of the 22 melted and mixed snow cores, as well as 8 randomly distributed samples, for subsequent analysis of NO 3 − to get the total deposition regardless of depth (which can vary significantly in concentrations levels (Kühnel et al 2011)). All samples were immediately melted and filtered using N-free filters, grade 521, made by VWR International. The weight of meltwater and filters were used to calculate the water equivalent per area. NO 3 − concentrations were analysed using a Metrohm 820 Advanced IC dual-channel system med 838 Advanced IC sample processor (detection limit=0.005 mg l −1 , precision=0.002 mg l −1 ). The true samples were measured in between blanks and standard samples, with a replicability <5%. NO 3 − concentrations and water equivalents were finally used to estimate fluxes across the landscape, assuming an even spatial distribution at <5 km scale (Björkman et al 2013). In addition to the measured nitrate concentration, we did a sensitivity test considering a 10% increase or decrease in concentration, in agreement with longer-term stratospheric concentration changes (Duderstadt et al 2014).

Estimation of runoff streams and landscape water retention
The snow meltwater runoff is modelled following the method developed by Balstrøm and Crawford (2018): the terrain model was examined for the presence of bluespots (landscape depressions) acting as local retention depots for meltwater and soluble ions. Each detected bluespot has a contributing catchment defined by the bluespot's topographic watershed (figure 1). When a bluespot's retention capacity is exceeded, the bluespot spills over at its pour point and carries the excess water volume downhill in accordance with the steepest local path. While the lake in the center of the study area is in principle a bluespot, we set its volume to zero assuming that the water table was stable throughout the study period. When the hydrologic screening is completed the output consists of a number of bluespots, their local contributing watersheds, pour points and streams stored as feature classes in ArcGIS Desktop 10.7 (https://www.esri.com). If two streams merge a pseudo pour point with zero retention capacity is established for the sake of the upcoming tracing procedure. Next, the streams are turned into a topologically This novel approach means that beyond the raster based calculations of the overall flow directions and identification of sinks from the terrain model having ∼187·10 6 cells, all further computations are handled in a simplified vector space. Thus, the total computation time for this high-resolution dataset in a combined ArcGIS/Python programming environment is only 4-5 CPU minutes on a desktop computer, which is critical to allow iterative sensitivity analyses. Here we run the model five times (measured SWE, plus/minus the SWE uncertainty, and a sensitivity test with plus/ minus 50% change in SWE relative to the measured). The hydrologic screening model assumes a Hortonian flow, i.e. surface runoff with no infiltration, which agrees with studies of snow meltwater runoff in areas with permafrost and/or frozen ground during initial snowmelt (Kuchment et al 2000, Johansson et al 2015. In addition, the period with snowmelt without refreezing of meltwater in the snowpack is only expected to last <9 d (Westergaard-Nielsen et al 2017). Consequently, we consider the meltwater from the snowpack to be released as one instant pulse and the remaining snow to be depleted with respect to watersoluble ions. In 2018, the majority of snow had melted away on the 10th of June (http://g-e-m.dk), during which the soil temperature in 1 cm depth was 0.05°C on average (1st-10th of June 2018) (http://g-e-m.dk). Consequently, the active layer was still frozen at that time, supporting our assumption of no infiltration.

Evaluating bluespot locations
Time integrated normalized difference vegetation index (TINDVI) is a commonly used spectrally derived proxy for vegetation productivity (Jia et al 2009, however it can also serve as an indicator of standing water for which the signal would be low or negative. As an indirect validation of the spatial distribution of bluespots, we computed TINDVI based on the atmospherically corrected L2A products from the Sentinel-2A and 2B satellites (being part of the ESA Copernicus program). The Sentinels are sun-synchronous Polar-orbiting resulting in overpass times ∼3-4 d at 69°N per satellite. We computed TINDVI using Google Earth Engine (https://earthengine.google.com) as the sum of NDVI values from days with <20% cloud cover from 10th to 30th of June 2018, corresponding to the melt season and early spring. In total our dataset was based on 15 overpasses. We then extracted TINDVI data from the areas outside and inside the bluespots' coverage for subsequent analyses.

UAS data
The average horizontal error of the snow-free and snow-covered surface models was 32 cm after GCP corrections. From the 22 validation points of snow core measurements (ranging from 0.17 to 1.55 m of snow) with DGPS locations, we had to omit two since they were in the periphery of the UAS-based surface models, where uncertainty is high due to a lack of overlapping images. Using the remaining 20 cores, we estimated the vertical error of the snow depth estimates to a mean absolute error of 9.8 cm (figure 2). The histogram for the photogrammetrically derived snow depths (figures 3(a) and (b)) shows that most of the depths are in the range of a few centimeters to 2 m, with a right-skewed distribution and a mode around 0.5 m. The average UAS-based snow depth is 1.08 m across the mapped area.

Snow properties
Due to the spatial uncertainty of the MagnaProbe depth measurements, it is not possible to compare those with the UAS-based snow depths. However, the statistical distribution is comparable to the manual probe measurements, although the probe is limited to 120 cm ( figure 3(c)).
Using the 22 snow cores from the Hydro-Snow probe, we measured an average snow density of 318 kg m −3 ±SD 55. We found a statistically significant linear relationship between snow depths and the snow samples' water equivalents (figure 4(a)) allowing for a subsequent conversion of UAS-based snow depths into SWE per grid cell ( figure 4(b)). Alongside the snow-free terrain model, this SWE data set formed the input to the Hortonian flow model.
The concentration of dissolved NO 3 − in the 30 independent samples was on average 0.095 mg l −1 ±SD 0.021, and the reported standard deviation shows that potential spatial differences in the concentrations are neglectable. The total estimate of deposited NO 3 − in the snowpack was 12.9 mg m −2 ±SD 1.8 or 17.4 kg±SD 3.9 for the entire catchment. The sensitivity test considering 10% increase or decrease in NO 3 − concentration but the same total amount of snow resulted in ±1.3 mg m −2 (0.29 mg N m −2 ), respectively, of NO 3 − change in deposition levels.

Water runoff modelling
Bluespot capacities 100 l were eliminated in the hydrologic screening considering the digital elevation model's overall vertical accuracy. From the derived watersheds it was computed whether the local bluespots' capacities were big enough to retain the meltwater within the bluespots. If not, the spill over was accumulated downstream. The modelled drainage pattern is presented in figure 5. 12 688 individual bluespots were identified with an average volume of 0.58 m 3 meltwater. The total water volume contained in bluespots was 7295 m 3 , +4 m 3 and −10 m 3 resulting from the 9.8 cm mean absolute error in snow depths. The disproportionate uncertainty is due to the vast majority of bluespots being filled up to maximum capacity, so increased meltwater will result in increased surface runoff rather than a higher retention. The total surface area of bluespots was 0.116 km 2 . Assuming equally   From the snow-free terrain model we see that the study site is separated into one prominent drainage catchment concealed by a lateral moraine towards North, a river system to the West (Røde Elv), and bedrock to the South and East (figure 5). Following those delineations, we divided the study area into four regions which drain to different catchments once the active layer thaws and allows for lateral water transport in the soil.

Relative spatial differences in TINDVI
The extracted TINDVI from the bluespots was significantly different from the average TINDVI in the studied area without bluespots (t-test, α=0.05). The TINDVI outside of bluespots followed a gaussian distribution with a mode at 2.05. However, the distribution for TINDVI from the bluespots was bimodal, with a gaussian mode at −1.34 and 1.73, respectively. Negative NDVI refer to objects with a higher absorption in the near-infrared than in the red part of the electromagnetic spectrum; it usually occurs over water bodies or highly moist surfaces. We found no correlation between TINDVI and area or volume of individual bluespots.  Using UAS data to estimate snow depths introduce uncertainty at several levels, including the precision of the structure-from-motion technique. Being a numerical optimization technique, the derived products are rarely 100% reproducible, and the precision depends on camera and image quality, image overlap, image contrast, type of surface features, and sensor-to-object distance (James andRobson 2012, Lucieer et al 2014a). The volume of bluespots derived from the snow-free terrain model (i.e. the surface model subtracted by vegetation heights) is consequently affected by this. In addition, the uncertainty in vegetation classification and heights should be considered. In this study area, the vegetation is <20 cm, with the exception of Salix glauca reaching up to 70 cm in certain spots. The terrain features, which the low vegetation generally follows, often have a greater difference in elevation across distances of a few meters, and the total elevation difference measured by the UAS is 183 m, illustrating that the bluespots are not a product of the photogrammetric surface model's vegetation canopy but rather the general terrain. Moreover, it is the relative elevation difference between the snow-free and snow-covered period that is important, which is ultimately reported by the validation against ground truth measurements of snow depths.

Discussion
We found our uncertainty in snow depth estimation on par with previous studies (Vander Jagt et al 2015, Harder et al 2016, Bühler et al 2017. More importantly, there is no significant difference in the estimates of retained snow meltwater in the landscape when taking the uncertainty in snow depths into account in the SWE-modelling. This results from the bluespots in the area being filled to maximum capacity already from 2% of the estimated meltwater volume in 2018. The measured snow density was only acquired at one point in time, although it is known to increase throughout the winter season due to compaction and refreezing of meltwater (Pedersen et al 2016). Here, we measured late-season densities which may overestimate the annual mean. To meet this, we measured densities in a transect of snow depths with an expected variation in compaction level. The linear correlation between depths and water content (figure 4(a)) suggests that such compaction was likely linear when the samples were taken. While the sample range cover <580 mm SWE and the modelled range is 0-2203 mm ( figure 4(b)), our sample range is covering 83% of the modelled SWEs. Hence, we argue to have a representative SWE dataset. The precision of the hortonian flow model is fully dependent on the accuracy of the terrain model, which in our case is a combination of uncertainty related to photogrammetry, vegetation heights, and a surface classification. Here, we validate the result based on two factors. Firstly, the modelled streams generally follow the visible streams in the area (figure 5) as well as the micro topography (e.g. lower lying routes between elevated patches with exposed bedrock). Secondly, the bluespots resulting from the flow model are co-located with negative TINDVI in early spring, which we interpret as presence of standing water. The lake is drained towards West, into the major river system denoted as 'Main westward streams' on figure 5. From visual inspection the lake had spill-over through these streams during and after snowmelt in the study period, and precipitation data show that the previous autumn had normal rainfall (Zhang et al 2019). We therefore assume that the lake was filled to maximum capacity throughout the study period.
Climate model projections of precipitation patterns in the Arctic suggest higher winter atmospheric moisture levels, with a possibility of increased snowcover (Cohen et al 2012). Yet, a trend analysis from regional climate model data of snowfall (MAR v3.6.4; Fettweis et al 2017) did not show any statistically significant changes in the period 1985-2016 for the studied site. The SWE sensitivity test showed that a higher ratio of the total deposited NO 3 − would run off with a 50% increase in SWE. However, the retained meltwater volume would only increase by 8 m 3 and the associated increase of NO 3 − less than 4 g for the total catchment. That means more snow with an unchanged mg l −1 NO 3 − concentration would not lead to a marked increase in N input to the terrestrial ecosystem in a landscape such as Blaesedalen because >95% of the nitrate-pool in the snow pack is lost to streams, lakes and other aquatic environments. If we assume that the total pool of nitrate in the snowpack is relocated to bluespots, the estimated total mass of retained NO 3 − will be~0.69 kg, corresponding to 11.7 mg m −2 if we consider the surface area of the bluespots only. The 10% sensitivity test of NO 3 − con- It is important to note that the previously mentioned ionic meltwater pulse is likely to have much higher ion concentrations in the first quarter of the meltwater (Rasch et al 2000, Elberling et al 2007, Boucher and Carey 2010. This suggests that the initial meltwater filling the bluespots can have a higher concentration of e.g. NO 3 − than the total meltwater volume on average. Therefore we may underestimate the levels of retained NO 3 − using our integrated method. Future climate changes could also result in increased autumn rainfall as shown by Bintanja and Andry (2017). Such rain events can be important for input and redistribution of NO 3 − and other nutrients partly due to soil infiltration but cannot be treated as simple as the initial meltwater, since infiltration may influence the vertical and horizontal fluxes. The effect of rain events in terms of NO 3 − wet-deposition will greatly depend on the timing and intensity. Highintensity rain would result in substantial runoff on the surface, while low intensity rain during the growing season could result in increased plant uptake. Another scenario with increased occurrences of winter warming events (Pedersen et al 2015) could lead to brief melt and runoff during winter. Under such circumstances our approach is expected to perform well, since the runoff would be hortonian on top of a frozen active layer.
The presented new combination of field measurements with UAS-based surveys and subsequent efficient runoff modelling show, that the spatial redistribution of NO 3 − can be estimated based on short field campaigns, which can aid the closure of the N budgets at catchment scale. The approach was successful used at the study site with a moderate topography and several gradients with more than 10% slope. The approach may not perform similarly well in large-scale catchments with flat topography due to the challenges in collecting precise ground control points for the UAS-survey, as well as the chance for meltwater infiltration. We therefore encourage testing of our approach on several contrasting sites with differences in topography and snow-cover to assess cross-site variability.

Conclusions
In this study, we successfully estimated a 3D snowpack at an Arctic tundra site, using UAS. We achieved a sufficient accuracy to model surface runoff pathways and landscape retention capacity of snow meltwater during spring thaw using a novel highly efficient network-based flow model. We measured similar nitrate concentrations in the pre-melt snow to other reported levels, and upscaled to a total nitrate stock in the snowpack. We found that the majority of the meltwater is drained from the terrestrial zone to streams and lakes, and with that, more than 95% of the atmospherically deposited nitrate in snow and likely other soluble nutrients. From this we conclude that future potential increases in snow fall does not necessarily lead to increased input of plant-available nitrate, as opposed to changes in atmospheric deposition rates or autumn rainfall which could potentially increase inputs. The method is generally applicable in environments with snowmelt occurring before the soil thaws to locate surface stream pathways and aid a closure of annual N budgets. We encourage use of this method in other landscape types to assess the variability of nitrate loss during spring thaw in permafrost areas.