## Abstract

Changes in active layer thickness (ALT) over Arctic and permafrost regions have an important impact on rainfall-runoff transformation. General warming is observed across Svalbard Archipelago and corresponds to increases in ground temperatures. Permafrost thaw and changes in ALT due to climate warming alter how water is routed and stored in catchments, and thus impact both surface and subsurface processes. The overall aim of the present study is to examine the relationships between temporal changes of active layer depth and hydrological model parameters, together with variation in the catchment response. The analysis was carried out for the small unglaciated catchment Fuglebekken, located in the vicinity of the Polish Polar Station Hornsund on Spitsbergen. For hydrological modelling, the conceptual rainfall-runoff HBV (Hydrologiska Byråns Vattenbalansavdelning) model was used. The model was calibrated and validated on runoff within subperiods. A moving window approach (3 weeks long) was applied to derive temporal variation of parameters. Model calibration, together with an estimation of parametric uncertainty, was carried out using the Shuffled Complex Evolution Metropolis algorithm. This allowed the dependence of HBV model parameters on ALT to be analysed. Also, we tested the influence of model simplification, correction of precipitation, and initial conditions on the modelling results.

## INTRODUCTION

Hydrological modelling in areas with the occurrence of perennially frozen ground together with seasonally changing active layer thickness (ALT) requires a better understanding of the influence of soil processes on water balance. Groundwater movement in permafrost terrain is limited because the frozen soil is practically impermeable and so infiltration and groundwater recharge are possible only within the active layer that thaws in summer or in taliks (Woo 2012). The impact of climate change on the ground thermal regime may affect active groundwater and hydrogeological processes directly and indirectly (Streletskiy *et al.* 2015; Walvoord & Kurylyk 2016). Dramatic environmental changes including higher air temperature, changes of precipitation patterns, permafrost degradation, and longer melting seasons have all been occurring recently in cold environment catchments (Callaghan *et al.* 2011; Bindoff *et al.* 2013; Bring & Destouni 2013; Walvoord & Kurylyk 2016; Osuch & Wawrzyniak 2017, 2018). The hydrological cycle in the Arctic is currently undergoing alterations due to different impacts (Vihma *et al.* 2016) and affects the global climate system (Bindoff *et al.* 2013). Despite widely documented evidence, remote catchments are still characterised by the poor availability of hydro-climatological data from *in situ* measurements, so a modelling approach is especially required for checking consistency in datasets, validation of assumptions, description, and parametrisation of processes which, by all means, improve our understanding of natural processes.

Hydrological models are tools for studying local, regional, or global catchment behaviour and runoff variability over changing climate conditions. Application of a hydrological model allows the diagnosis of a hydrological system (Carey & Woo 2001; Pomeroy *et al.* 2007; Gupta *et al.* 2008; De Vos *et al.* 2010; Clark *et al.* 2015; Tetzlaff *et al.* 2015; Guse *et al.* 2016; Pianosi & Wagener 2016; Krogh *et al.* 2017). Understanding the dominant processes and their temporal variations could be performed by analysis of temporal variation of model parameters and their influence on runoff (Osuch *et al.* 2015; Guse *et al.* 2016; Pianosi & Wagener 2016). Such temporal variation of dominant processes is particularly evident in the case of Arctic catchments with well-pronounced seasonality. Mathematical modelling of hydrological processes in these catchments demands more information than only basic meteorological inputs and catchment attributes. It has to take into consideration the different levels of interaction between groundwater and surface runoff, which is a function of the seasonal change in the ALT of the catchment.

In this paper, the conceptual catchment runoff HBV (Hydrologiska Byråns Vattenbalansavdelning) model (Bergström 1995) was applied for modelling of hydrological processes within the small, non-glaciated Arctic catchment Fuglebekken. The HBV model was used as a diagnostic tool to provide information about dominant processes and water storage in the catchment. This variability is associated with differences in water storage in snow pack, soil moisture, fast runoff, and slow runoff, and these differences, in turn, have important implications for the ground thermal regime.

To analyse the temporal variability of such processes the model was calibrated using 3-week long moving windows in four study periods (2014–2017) with varying hydro-meteorological conditions. The HBV model calibration was carried out using a Shuffled Complex Evolution Metropolis algorithm (SCEM-UA; Vrugt *et al.* 2003) that allows for estimations of parametric uncertainty. A special goal of this paper was to investigate the dependence of the HBV model on hydro-meteorological indices, and in particular ALT, which is subject to large changes in the study area (Wawrzyniak *et al.* 2016). These analyses are especially important for simulations of future water balance while using climate conditions projections (Osuch & Wawrzyniak 2016).

As simulation studies strongly depend on the applied tools (hydrological models), we also investigated the influence of model simplification, correction of precipitation, and initial conditions on the modelling results.

## STUDY AREA

The unglaciated Fuglebekken catchment (1.27 km^{2}) is located near the Polish Polar Station Hornsund (PPS; 77°00′N, 15°30′E) on Spitsbergen and drains to the northern shore of the Hornsund fjord (Figure 1). The studied catchment elevations range from 522 m a.s.l. just below the summit of Fugleberget, to 4 m a.s.l. at the runoff measurement point. The catchment is heterogeneous regarding topography and includes the steep slopes of the Ariekammen-Fugleberget massif and relatively flat Fuglebergsletta plain with raised marine terraces.

The Fuglebekken runoff is characterised by seasonal variability, with freezing from late autumn to early spring, and with typical large variation in the runoff during the ablation season. The hydrologically active flow season can be divided into three periods: high runoff during the snowmelt (until mid-July), decrease in the middle of summer, and large runoff in the late summer–autumn (September) period due to intense precipitation, as recently observed (Polkowska *et al.* 2011; Wawrzyniak *et al.* 2017). The results of ground thermal modelling at the study area were presented by Wawrzyniak *et al.* (2016), including the derived subsurface temperature variability together with the evolution of ground thermal regime, which show an increase of the ground temperatures in two recent decades. Permafrost and active layer conditions at a catchment are affected by the local boundary conditions between the atmosphere and the ground, including air temperature, soil properties, geothermal heat, and snow depth. Their high spatial and temporal variability influences ALT and makes it hard to fully understand the small-scale complexity of hydrogeological processes.

Meteorological data used for modelling purposes in this study come from the meteorological site at PPS Hornsund (10 m a.s.l.), where data have been collected continuously since July 1978. Datasets include all basic meteorological parameters, such as air temperature, humidity, precipitation, wind speed, wind direction, atmospheric pressure, solar radiation, snow depth, snow water equivalent, and ground temperatures. To model in this study we used 1-minute interval data on air temperature (sensor HMP155) and precipitation (OTT Parsivel), together with potential evapotranspiration calculated using the Hamon method (Hamon 1961) and ground temperature up to 50 cm (3 hours interval).

Runoff measurements at Fuglebekken in 10-minute intervals started in June 2014, and since then have been continuously conducted during each season over periods without snow cover. Snow and ice in the river channel create disturbances for runoff measurements, especially during a brief but intense snowmelt period. Runoff was measured in the same profile (77.0053°N, 15.5561°E) using a portable flowmeter Nivus PCM-F, with Active Doppler sensor (KDA-KP 10). The final value of velocity and water level uncertainty was ±1 and ≤0.5%, respectively.

The original data used for further analysis are presented in Figure 2. A comparison of air temperature input indicates small differences between the four seasons. In July 2015, the maximum long-term air temperature was recorded at Hornsund station. In the case of precipitation, a higher amount was observed in September than in the two previous months. In each study season from the second half of July to the second half of August there was almost no precipitation. The highest peak of precipitation was recorded in September 2017.

Observed runoff is presented in the lowest panel of Figure 2 in logarithmic scale in order to highlight low flows and the recession curve. The peak in runoff corresponds well with peaks in precipitation, and the highest values are observed in September. The driest of the four analysed seasons was 2014, with lack of runoff in the second half of August. As in the case of precipitation, the highest runoff in the analysed period was in September, except 2016.

## METHODS

### Catchment runoff models

In this study, the conceptual HBV model was applied, which is a well know rainfall-runoff model introduced by the Swedish Meteorological and Hydrological Institute (Bergström & Forsman 1973; Bergström & Lindström 2015) and currently available in different forms and versions used for an increasing number of applications and purposes (Bergström 1995; Lindström *et al.* 1997; Booij 2005; Booij & Krol 2010; Romanowicz *et al.* 2013; Osuch *et al.* 2015). A detailed description of the applied version of the model is presented by Piotrowski *et al.* (2017). A scheme of the applied HBV model is presented in Figure 3. It consists of five conceptual storages (snow pack, melt water, soil moisture, fast runoff, and slow runoff) and the dominant process of rainfall-runoff transformation. The required input data are precipitation, air temperature, and potential evapotranspiration. As an output time-series of runoff, evapotranspiration, snowpack storage, fast and slow runoff reservoir, and soil moisture are calculated; however, usually just runoff simulations are taken into account. In this study, the simulation of water storage was also analysed.

In the full version, the model is parametrised with 14 parameters, for which values are estimated by calibration. Six parameters are related to snow processes: *TT*, *TTI*, *CFMAX*, *CFR, DTTM,* and *WHC*. Three parameters, *FC*, *LP,* and *β,* are related to the third conceptual storage (soil moisture). *CFLUX* represents the capillary transport between the fast runoff reservoir and the soil moisture reservoir. *KF* and *α* parameters describe the fast runoff. Percolation is the transport between the fast and slow runoff reservoirs and is described by the parameter *PERC*, while the slow runoff reservoir is represented by the *KS* parameter. In addition, an initial condition of the five storages (snow pack, melt water, soil moisture, fast runoff, and slow runoff) should be quantified as their values influence modelling results.

The applicability of the HBV model to simulate runoff in the Fuglebekken catchment was previously tested and presented by Wawrzyniak *et al.* (2017). The results of that study indicated that the calculated fit of simulated to observed runoff depends on the year, time step, and data averaging. The best results were obtained for the model from the year 2015 for 3 and 6 hours using averaged input data. Therefore, in this study, we use a 3-hour averaged time step.

### Modelling procedure

Model calibration was carried out using the SCEM-UA algorithm (Vrugt *et al.* 2003). This method combines controlled random search, competitive evolution, and complex shuffling together with the Monte Carlo Markov Chain (MCMC) algorithm for estimation of posterior parameter distributions. The settings used for the calibration run are based on the recommendation of Vrugt *et al.* (2006) for complex problems (number of complexes/sequences is five, population size is 250, the maximum number of iterations is 100,000). MCMC was carried out using an informal likelihood measure, i.e. the sum of squared errors (denoted as SSE). The procedures of SCEM-UA (the evolution and shuffling) were repeated until the Gelman–Rubin convergence diagnostic for each of the parameters confirmed convergence to a stationary posterior distribution (Gelman & Rubin 1992). Convergence was achieved in all analysed cases.

The model parameters were drawn from uniform distributions within predefined upper and lower parameter boundaries (Table 1) based on the previous findings described by Wawrzyniak *et al.* (2017). The HBV model was calibrated using available precipitation, air temperature, PET, and runoff observations at 3-hour time steps from four study seasons (2014–2017). To analyse temporal variability in the model parameters, a 3-week long moving window was applied changing on a daily basis, i.e. the window moved day-by-day between two consecutive periods, and, therefore, the subperiods overlapped. Due to the different lengths of the datasets of runoff measurements in subsequent years, there are a different number of periods for each study season: 57 (2014), 58 (2015), 60 (2016), and 67 (2017). In all cases, the first 7 days were not taken into account for the calculation of objective function to avoid the influence of initial conditions on simulated runoff.

No. . | Parameter (unit) . | Description . | Minimum . | Maximum . | Simulation type . |
---|---|---|---|---|---|

1 | FC (mm) | Field capacity, maximum soil moisture storage | 10 | 1000 | 7, 8, 9, 15 |

2 | β (–) | Non-linearity parameter | 0.01 | 7 | 7, 8, 9, 15 |

3 | LP (–) | Factor limiting potential evapotranspiration | 0.1 | 1 | 15 |

4 | α (–) | Non-linearity parameter | 0.1 | 3 | 7, 8, 9, 15 |

5 | KF (1/period^{a}) | Fast runoff parameter | 0.0005 | 0.3 | 7, 8, 9, 15 |

6 | KS (1/period) | Slow runoff parameter | 0.0005 | 0.3 | 7, 8, 9, 15 |

7 | PERC (mm/period) | Rate of percolation | 0.001 | 6 | 7, 8, 9, 15 |

8 | CFLUX (mm/period) | Rate of capillary rise | 0 | 4 | 7, 8, 9, 15 |

9 | TT (°C) | Threshold temperature for snowfall | –3 | 3 | 15 |

10 | TTI (°C) | Threshold temperature interval length | 0 | 7 | 15 |

11 | CFMAX (mm/°C period) | Degree day factor, rate of snowmelt | 0 | 20 | 15 |

12 | DTTM (°C) | Value to be added to TT to give the threshold temperature for snow melt | 0 | 1 | 15 |

13 | CFR (–) | Refreezing factor | 0 | 1 | 15 |

14 | WHC (mm/mm) | Water holding capacity of snow | 0 | 0.8 | 15 |

15 | wprecip (–) | Weight to correct precipitation | 0 | 5 | 8, 9, 15 |

16 | gw1 (mm) | Initial water storage in slow runoff reservoir | 0 | 1000 | 8, 9 |

No. . | Parameter (unit) . | Description . | Minimum . | Maximum . | Simulation type . |
---|---|---|---|---|---|

1 | FC (mm) | Field capacity, maximum soil moisture storage | 10 | 1000 | 7, 8, 9, 15 |

2 | β (–) | Non-linearity parameter | 0.01 | 7 | 7, 8, 9, 15 |

3 | LP (–) | Factor limiting potential evapotranspiration | 0.1 | 1 | 15 |

4 | α (–) | Non-linearity parameter | 0.1 | 3 | 7, 8, 9, 15 |

5 | KF (1/period^{a}) | Fast runoff parameter | 0.0005 | 0.3 | 7, 8, 9, 15 |

6 | KS (1/period) | Slow runoff parameter | 0.0005 | 0.3 | 7, 8, 9, 15 |

7 | PERC (mm/period) | Rate of percolation | 0.001 | 6 | 7, 8, 9, 15 |

8 | CFLUX (mm/period) | Rate of capillary rise | 0 | 4 | 7, 8, 9, 15 |

9 | TT (°C) | Threshold temperature for snowfall | –3 | 3 | 15 |

10 | TTI (°C) | Threshold temperature interval length | 0 | 7 | 15 |

11 | CFMAX (mm/°C period) | Degree day factor, rate of snowmelt | 0 | 20 | 15 |

12 | DTTM (°C) | Value to be added to TT to give the threshold temperature for snow melt | 0 | 1 | 15 |

13 | CFR (–) | Refreezing factor | 0 | 1 | 15 |

14 | WHC (mm/mm) | Water holding capacity of snow | 0 | 0.8 | 15 |

15 | wprecip (–) | Weight to correct precipitation | 0 | 5 | 8, 9, 15 |

16 | gw1 (mm) | Initial water storage in slow runoff reservoir | 0 | 1000 | 8, 9 |

The last column denotes simulation type (number of model parameters).

^{a}Three-hour interval in this study.

The selected version of the HBV model includes 14 parameters and could suffer from equifinality problems (Beven 2006). To simplify the model, the results of parametric uncertainty using the SCEM-UA method were applied. Due to the different ranges and units of the parameters, a comparison of relative standard deviation (RSD, the ratio of the standard deviation to the mean) was performed. The smallest values of the RSD were obtained for the *FC*, *β*, *α*, *KF*, *KS*, *PERC*, and *CFLUX* parameters. Therefore, these seven parameters were selected, and the HBV model with seven parameters was analysed. Snowmelt and snow accumulation are very important for the modelling of hydrological conditions at Svalbard, although only for a limited period as presented in Osuch & Wawrzyniak (2017) and Kepski *et al.* (2017). The period between July and September is without snow cover in the Fuglebekken catchment, so the influence of parameters connected with snowmelt and snow accumulation on the modelling results is irrelevant in the summer season. In the simplified versions of the HBV model the parameters related to snowmelt and snow accumulation were assumed constant and were not optimised.

In addition, we tested the influence of the precipitation time series correction on the calibration results and temporal variability of the model parameters by analysis of the HBV model with and without a precipitation correction parameter. In the former case, observed precipitation was multiplied by an additional model parameter denoted as *wprecip*. It was assumed that the values of this parameter are constant in the analysed moving window, which means that they are not attributed independently to each precipitation event (Kavetski *et al.* 2002, 2006). Precipitation correction by multiplication causes no change in the case of observed zero precipitation. This approach provides an important new way to estimate and correct areal average precipitation over the catchment, information that is very important for hydrological modelling. Precipitation time series for the Fuglebekken catchment were developed on the basis of point measurements conducted at the meteorological site at PPS Hornsund, located some 500 m from the gauging station. As mentioned before, the catchment is heterogeneous in terms of topography and includes a mountain massif (maximum elevation 522 m a.s.l.), so the precipitation gradient has to be taken into consideration. As it is not available, the catchment average precipitation was a subject of correction by weight optimisation.

To avoid the influence of initial conditions, the proposed procedure assumes that simulations from the first week in each moving window were not taken into account during calculation of the objective function. To confirm this assumption (1 week excluded), we tested a version of the HBV model with initial conditions of slow runoff reservoir (*gw1*) taken as an additional parameter. This parameter was calibrated together with the others. Further comparisons of the calculated fit and model parameters allowed for the determination of the influence of *gw1* on the results. These outcomes indicated that the influence of initial storage in other reservoirs (e.g. snow pack and melt water) on the model performance is visible in a shorter period than 1 week. Therefore, these values were not optimised.

Altogether, four versions of the HBV model were applied:

- A.
Simplified model with seven parameters

*FC*,*β*,*α*,*KF*,*KS*,*PERC*, and*CFLUX*. - B.
Simplified model with eight parameters. Seven parameters as in case A and one additional parameter (

*wprecip*) for correction of precipitation. - C.
Simplified model with nine parameters. Eight parameters as in case B and one parameter that represents initial conditions in slow runoff reservoir (

*gw1*). - D.
The full version of the HBV model (14 parameters) with one additional parameter (

*wprecip*) for correction of precipitation.

These model versions were selected to analyse the influence of the model simplification, correction of precipitation, initial conditions on model performance, and also the temporal variability of model parameters.

### Hydro-meteorological indices

To explain the temporal variability of the HBV model, a set of hydro-meteorological indices was selected following the literature review (Merz *et al.* 2011; Woo 2012). The most important variables are precipitation and air temperature, but also ALT. Using air temperature and precipitation measurements at PPS Hornsund, four indices were calculated: mean air temperature, total precipitation, number of days with precipitation, and variability of precipitation in a given moving window. As observations of ALT are not available for Hornsund for all seasons (lack of observations below 1 m depth), the temporal variability of ALT was calculated by a simulation method. For this purpose, the CryoGrid 2 (Westermann *et al.* 2013) model was applied. The CryoGrid 2 model describes ground temperatures based on time series of air temperature and snow depth in accordance to conductive heat transfer in the soil and the snowpack. The model does not take into account the movement of water or water vapour in the ground, so the liquid water content can only change over time due to freeze-thaw processes. The CryoGrid 2 is a one-dimensional model; therefore, only heat flow in the vertical direction is considered, neglecting horizontal heat flow. In order to simulate ground temperatures at the PPS Hornsund, the CryoGrid 2 was calibrated and validated against borehole temperatures (Wawrzyniak *et al.* 2016). The results of a previous modelling study allowed for simulation of ground temperatures in the Fuglebekken catchment and estimation of ALT.

As calculated values of ALT are characterised by additional uncertainties, observed mean ground temperatures at 5, 10, 20, and 50 cm as potential variables that shape hydrological model parameters were also included. To assess the influence of water on the variability of model parameters, water storage in a slow runoff reservoir was also taken into account. The values of this index were quantified by the parameter *gw1* of the HBV model.

For each 3-week long time window, the following hydro-meteorological indices were determined: (1) mean air temperature, (2) precipitation total, (3) standard deviation of precipitation, (4) number of days with precipitation, (5) mean ground temperature at 5 cm, (6) mean ground temperature at 10 cm, (7) mean ground temperature at 20 cm, (8) mean ground temperature at 50 cm, (9) mean ALT (estimated using CryoGRID2 model), and (10) slow runoff reservoir storage (estimated as value of the *gw1* parameter of the HBV model).

### Dependence of hydrologic model parameters on hydro-meteorological indices

To assess the dependence of the HBV model parameters on hydro-meteorological indices, a weighted version of the Spearman correlation coefficient (WSCC; Bailey *et al.* 2018) was applied to take into account parametric uncertainty and nonlinearity of relationships. The weight structure was constructed by converting the standard errors of parameter estimates calculated by the SCEM-UA method into inverse-variance weights.

In this way, information on the uncertainty of the HBV model parameters was included. Weighted Spearman correlation coefficients (WSCC) are very similar to the Spearman correlation coefficient (SCC) and characterised by the same properties. Both range from −1 to 1. A value of 0 indicates no correlation while values of −1 or 1 indicate a perfect negative or positive correlation, respectively.

## RESULTS

This section presents the results of the calibration and validation of model parameters performed for the selected study seasons, together with the relationships between the hydro-meteorological indices and the model parameters optimised for moving 3-week subperiods.

As we investigated the dependence of HBV model parameters on hydro-meteorological indices, time-variability of indices was analysed first, followed by the calibration, temporal variability of the HBV model parameters, and the discussion of the relationships between the HBV model parameters and hydro-meteorological indices.

### Temporal variability of hydro-meteorological indices

A comparison of temporal variability for four hydro-meteorological indices is presented in Figure 4. The upper left panel shows mean air temperature in four study seasons: 2014, 2015, 2016, and 2017. The similarity of these time series is visible with only slight differences between seasons. At the beginning of the study period, higher air temperatures were observed in the year 2016 (6.5 °C) than the other three seasons. There were differences in seasonal maximum mean air temperature between years: at the beginning of July in 2014 and 2016, and around the middle of July in 2015 and 2017. In three of the analysed seasons, a high decrease in mean air temperatures in September was visible, except in 2017 when September was warmer than in previous years. Differences in mean ground temperatures (lower left panel) and the active layer depth (lower right panel) were strongly related to the differences in mean air temperature between seasons. The highest ground temperatures and the ALT were noted in 2016.

The four study seasons differ from each other when taking into account the temporal variability of precipitation (upper right panel of Figure 4). The second part of June and the first part of July in 2014 and 2016 were characterised by higher precipitation amounts than in 2015 and 2017. No precipitation was recorded in the second part of July and August in 2014 and 2017. High precipitation in September was observed in 2014, 2015, and 2017.

Further, we studied the influence of these varying hydro-meteorological conditions on the HBV model parameters.

### HBV model calibration

*et al.*2009) were also calculated. This criterion was calculated according to the following formula: where subscripts

*s*and

*o*denote simulated and observed time series,

*KGE*is the linear correlation coefficient,

_{r}*KGE*is a measure of relative variability in the simulated and observed time series, and

_{α}*KGE*is the ratio between the mean simulated and observed values. This decomposition into

_{β}*KGE*,

_{r}*KGE*, and

_{α}*KGE*allows for detection of the main sources of problems during calibration. In the ideal case, the three components of decomposed KGE are

_{β}*KGE*= 1,

_{r}*KGE*= 1, and

_{α}*KGE*= 1. The results of the decomposition show that the

_{β}*KGE*values in the calibration period are always below one, i.e. the variability is systematically underestimated (Gupta

_{α}*et al.*2009). Good calibration results are required to derive potentially good parameter values and assess their temporal variability and correlation with climatic variables.

The comparison of the calibration results from 3-week moving windows for four study seasons is presented in Figure 5. It is shown that the performance depends on: moving window, season, and applied version of the HBV model with different numbers of parameters. The best calibration results were estimated for the model with nine parameters in study seasons 2014, 2015, and 2016 with mean KGE values of 0.5043, 0.4569, and 0.0368, respectively. In the case of season 2017, the largest mean KGE (0.2997) was estimated for the model with 15 parameters.

An influence of correction of observed precipitation on model performance was assessed by comparison of the results for models with seven and eight parameters. The largest improvement (an increase of the KGE values) was obtained for seasons 2014 (for the first 15 moving windows and from 40 to 57), 2015 (from 34 to 58), 2016 (from 31 to 47), and 2017 (from 47 to 56 and from 62 to 67). The periods with differences in model performance correspond to periods with high amounts of precipitation. A small decrease in model performance was found for just a few moving windows in the study season 2014.

The comparison of obtained KGE values between simulations of models with eight and 15 parameters indicates that in most cases model simplification does not influence the KGE values, except the beginning of the study season in years 2014 and 2017. This may be a result of snowmelt contribution from the upper part of the catchment. In some cases, better results were obtained for the full model (2016 from 31 to 42 and 2017 from number 51 to 56 of the moving windows). The opposite situation was calculated for 2014 from moving window number 40 to 57, 2015 from 34 to 47 and from 53 to 58, and 2017 from 62 to 67. In general, the results of the full model with 15 parameters were characterised by a large variation in model performance, while the simplified model gives a smaller variability of fit.

Influence of initial conditions in slow runoff reservoir on the model fit was analysed by comparison of models with eight and nine parameters. The model with nine parameters that included parameter *gw1* improved the fit for the following moving windows: in the year 2014 from 14 to 50, in 2015 from 19 to 28 and from 45 to 52, in 2016 from 1 to 40, and in 2017 from 51 to 60. The results for other time windows are almost the same, except a few moving windows in the study season 2014 when the application of the model with eight parameters gave better results than the model with nine parameters.

### Temporal variation of HBV model parameters

The posterior distributions of HBV model parameters were used to estimate the mean and standard deviation of the obtained parameters for each study season, each simulation type, and each 3-week period.

The temporal variability of the HBV model parameters and their uncertainty for the Fuglebekken catchment in 2015 are shown in Figure 6. The results for the remaining three study seasons and simulation types are presented in the Supplementary material (available with the online version of this paper). Panels represent results for each parameter. The dots denote the mean values of posteriors and whiskers represent ± standard deviation of posteriors. The X-axis shows the dates of the beginning of the 3-week long moving windows.

The outcomes indicate the relatively small variability of the optimised parameters, especially for the first part of the study season 2015. For *FC*, *β,* and *CFLUX,* a small decrease of the mean for moving windows starting with August 3rd and later was found. For the same period, a small increase in *wprecip* was estimated (change from 2.5 to 3.0). In the case of *α* and *KF,* larger changes (decreases) were visible. The results for the *KS* parameter showed two-modal behaviour. The first maximum is at the beginning of the study season and the second is in the fifth to last moving window. An explanation of the temporal variability of the estimated parameters is presented in the following sections.

Application of the SCEM-UA method allowed for estimation of parametric uncertainty that is visualised by whiskers that represent ± standard deviation. In most cases, relatively large parametric uncertainty was found for the relatively short window (3 weeks long) used in this study. The smallest uncertainty was estimated for *KS* parameter. A comparison of results for four years and four simulation types indicated similarities between them. We also tested the significance in mean and variance of the estimated parameters, taking into account autocorrelation and parametric uncertainty. For this purpose, three tests were applied: two sample t-test with unequal variances, the Wilcoxon rank sum test, and the two-sample F-test. The tests were performed for the selected independent samples (moving window 1 and 40, 20 and 50). The results of these tests indicate rejection of the null hypothesis, the presence of statistically significant differences in mean and variance at the 5% significance level. Similar results were obtained in all tested cases: for different parameters, models, and seasons.

### Temporal variability of *KS* parameter

Temporal variability of *KS* parameter was studied in detail as it is characterised by the smallest parametric uncertainty and potentially highest correlation with the ALT, according to the literature review (Woo 2012). A comparison of the estimates of *KS* values between four tested versions of the HBV model for the study season 2015 is presented in Figure 7. In most cases, results were very similar except for outcomes for the model with nine parameters, where larger uncertainty and variability of *KS* were found. These differences were probably the result of choosing *gw1* as an additional model parameter.

A comparison of KS values estimated for four study seasons, 2014, 2015, 2016, and 2017, using the HBV model with eight parameters is presented in Figure 8. In most cases, *KS* values vary from 0.0005 to 0.0100 (1/period), except one moving window in the year 2014. The shape of these runs differs between study seasons probably due to changing ALT and meteorological conditions.

### Estimates of water storage in the slow runoff reservoir

The HBV model contains five conceptual storages. One of them is a slow runoff reservoir that has a single influx from percolation and a single outflow as slow runoff (*q _{s}*), which is parametrised with

*KS.*Runoff from this storage also depends on water storage in this reservoir (

*gw1*). The calculated water storage in this slow runoff reservoir was assessed in four study seasons and is presented in Figure 9. The results strongly depend on the season, because the timing and availability of water vary from year to year; however, a similarity in the annual run was also found. Usually, two maxima of

*gw1*are visible. The first maximum is observed at the beginning of the study season due to snowmelt and ice-rich permafrost thawing. The second maximum is noted in the second part of August and in September as a result of the high precipitation amount.

A comparison of the estimated values of *gw1* indicates large differences between years (Figure 9). An increase in time for these four seasons is visible. The smallest *gw1* was simulated in the year 2014 while the largest was in the year 2017. These differences correspond to precipitation amount in these seasons. Owing to drought in 2014 (for a long period there was no precipitation at all), *gw1* is the smallest in this season. In the case of 2017, the opposite situation occurred, and the highest amount of precipitation was recorded.

### Temporal variability of *wprecip*

Three versions of the HBV model use correction of precipitation (*wprecip*) as an additional parameter. The precipitation time series were corrected by multiplication of the observed time series by *wprecip*. As presented in Figure 5, this parameter improves model performance. The present modelling study allows for the quantification of such correction. Temporal variability of *wprecip* calculated with the help of the model with 8, 9, and 15 parameters for the year 2015 is presented in Figure S5 in the Supplementary material (available with the online version of this paper). It is shown that for the first part of the study season similar values of *wprecip* (around 2.5) were estimated. In the second part of the season, *wprecip* values varied from 2.5 to 3.4; however, consistency of the results between simulation types is preserved. Similar estimates of *wprecip* were obtained for other study seasons (Figures S4, S6 and S7 in the Supplementary material).

### Dependence of model parameters on hydro-meteorological indices

In the next step, the dependence of the HBV model parameters on hydro-meteorological indices using the weighted SCC (WSCC) was estimated. The analyses were performed for four study seasons and four versions of the HBV model for each of them (altogether 16 cases). A comparison of the WSCC indicated that there were larger differences between seasons than for different HBV versions for the same season. Therefore, we present outcomes for all seasons for the HBV model with nine parameters only. The results estimated for the study season 2015 and hydro-meteorological indices are presented in Table 2. The outcomes for the other three study seasons are available in the Supplementary material.

The boxes are colour coded based on the degree of the positive or negative correlation coefficient. Please refer to the online version of this paper to see this table in colour: http://dx.doi.org/10.2166/nh.2019.031.

^{a}The values of these hydro-meteorological indices were estimated using CryoGrid 2 and HBV models, respectively.

The results shown in Table 2 indicate that all parameters correlated with hydro-meteorological indices, even for parameters with small temporal variability like *FC*, *β*, *PERC*, and *CFLUX*.

There is no simple explanation for the temporal variability of the parameters due to the complex interactions between hydro-meteorological indices and water cycle elements in the catchment. In addition, the results from the four study seasons differ between each other probably due to the variability of hydro-meteorological conditions and their impact on ground water storage. Relatively similar results were obtained for the years 2014, 2015, and 2017, while 2016 was different. These differences in 2016 may be explained by higher mean air temperatures, mean ALT, and lack of high precipitation at the end of the study season.

In the years 2014, 2015, and 2017, two parameters characterising fast runoff routine, *KF* and *α,* were highly correlated with the sum of precipitation. There was also a high correlation between *KF* and mean air temperature, mean ground temperatures, and mean ALT.

Parameter *KS* is correlated with water storage in the slow runoff reservoir in the years 2015 (−0.8711), 2016 (−0.8736), and 2017 (−0.8736). In season 2014, the correlation coefficient was smaller, probably due to the small amount of water in this reservoir. A comparison of the correlation coefficient between parameter *gw1* and mean ALT between years gives some information about the quality of simulations of ALT by CryoGrid 2. The significantly higher correlation was obtained for the years 2014 and 2015 with relatively small ground water storage in comparison to the years 2016 and 2017. As CryoGrid 2 does not take into account changes of water content in the ground due to percolation, the estimates of ALT for the years 2016 and 2017 may be underestimated.

Estimated values of percolation (*PERC*) between fast and slow runoff reservoirs were negatively correlated with mean ground temperatures in seasons 2014, 2015, and 2017.

A comparison of the correlation coefficient between parameter *wprecip,* which is responsible for correction of precipitation, indicates large differences between seasons. In seasons 2015, 2016, and 2017, *wprecip* was negatively correlated with air temperature and *gw1*.

In order to explain the temporal variability of *KS* parameter in four study seasons, the interactions of *KS* and hydro-meteorological indices were studied. Figure 10 shows *KS* variability associated with differences in mean air temperature, the sum of precipitation, mean ground temperature at 50 cm, and simulated ALT. A clear relationship of *KS* with mean air and ground temperatures and associated mean ALT is visible. These indices are characterised by well-established seasonal increase and decrease, which are also reflected in the relationship with *KS* (part of a loop).

Large differences in these relationships between seasons were found, which probably result from complex interactions within the catchment system. There is a similarity in the results for *KS* and ALT for seasons 2014 and 2015, and for 2016 and 2017. We suspect that ground water storage was significantly higher in 2016 and 2017 than in 2014 and 2015, and this strongly influences these relationships.

## DISCUSSION AND CONCLUSIONS

In this study, we present the analysis of the temporal variability of hydro-meteorological conditions and their influence on hydrological processes in the small Arctic, unglaciated catchment in four study seasons, 2014–2017. For this purpose, the conceptual HBV model was applied. This model was calibrated by the SCEM-UA Bayesian optimisation method for each season using 3-week long moving windows. The applied method allows for estimation of uncertainty that arises from a variety of sources, such as data input error, calibration accuracy, parameter sensitivity, and parameter uncertainty. The proposed approach is crucial for analysis of the temporal variability of model parameters as well as catchment response to changing ALT.

In addition, we studied the influence of model simplification, correction of precipitation, and initial conditions on the quality of the modelling results. The obtained results indicate that both the full and the simplified models give almost the same values of the objective function. Calibration results of the full model were characterised by a large variation of model performance, while the simplified model gives a smaller variability of fit. It can be explained by an overparametrisation problem and the decrease of the power of the optimisation method due to the increase in model dimensionality (Vrugt *et al.* 2003). The obtained results (lack of the influence of parameters related to snow melting and accumulation) can be explained by the period under study. Snowmelt and snow accumulation are very important for the modelling of hydrological conditions at Svalbard but only for a limited time. The duration of the snow-free period in the Fuglebekken catchment is short and lasts typically from the beginning of July until the end of September (Kepski *et al.* 2017), so in this period the influence of parameters connected with snowmelt and snow accumulation on the modelling results is irrelevant, and the model was simplified.

Estimation of the average precipitation in the mountainous catchments is the subject of the largest uncertainty compared to the other hydro-meteorological variables. Catchment total precipitation is very important for the quality of hydrological models. The application of precipitation correction in the HBV resulted in improvement of the model fit. The calculated values of the precipitation correction parameter (*wprecip*) were consistent between moving windows and equal to around 2.5. The precipitation time series for the Fuglebekken catchment were developed on the basis of point measurements conducted at the PPS Hornsund meteorological site, located at 10 m a.s.l., 500 m from the gauging station. The catchment is heterogeneous regarding topography and includes a mountain massif (maximum elevation 522 m a.s.l.), so the precipitation gradient should be taken into account (Killingtveit *et al.* 2003). It has to be mentioned that not only does the precipitation gradient forced by topography influence the time series of precipitation but also other factors such as wind-induced undercatch, wetting loss, evaporation loss, and underestimation of trace precipitation amounts. The study seasons from June to September cover periods with a limited amount of precipitation. For most of the time it is trace precipitation only, therefore, derived precipitation correction results in relatively small absolute changes. The water fluxes and storages are simulated not observed, so precipitation correction also depends on the model and its formulation. Presented results confirm that hydrological models could help in the appropriate estimation of water balance elements (Vrugt *et al.* 2006).

Influence of the initial conditions on model performance was tested. The results indicated that in the case of the Fuglebekken catchment, the influence of initial conditions could be neglected taking into account a 1-week long warm-up period, except ground water storage (slow runoff reservoir). Incorporation of initial conditions of the slow runoff reservoir storage as an additional model parameter also improved obtained values of the objective function (model performance). The other option could be the prolongation of the warm-up period, but it was not possible for this study due to the limited time series of runoff. Analysis of temporal variability of HBV model parameters showed large consistency in the results, however, statistically significant differences in mean and variance between moving windows were found. The outcomes for particular moving windows were similar between study seasons regarding posterior means and standard deviations. The smallest parametric uncertainty was found for the slow runoff parameter (*KS*). Our results indicate that *KS* is the main factor influencing runoff in the Fuglebekken catchment. Parametric uncertainty in the case of other parameters was larger due to their smaller influence on the modelling results. These parametric uncertainties may result from the relatively short period (3 weeks) at the calibration stage.

Temporal variability of the *KS* parameter was characterised by the presence of two maxima: the first maximum at the beginning of the study season (due to snowmelt and ice-rich permafrost thawing) and the second maximum in September (a result of high amount of precipitation). Changes of hydro-meteorological conditions in the four study seasons were sources of differences between *KS* runs. The correlation analysis showed that the values of *KS* were associated with mean air temperature, the sum of precipitation, mean ground temperature at 50 cm, and ALT. Large differences in these relationships between seasons were found, probably as a result of complex interaction in elements of the catchment system. We suspect that significantly higher ground water storage (*gw1*) in 2016 and 2017 than in 2014 and 2015 strongly influenced these relationships. Temporal variability of other parameters was smaller and usually within their parametric uncertainty.

A dependence of the HBV model parameters on hydro-meteorological indices was analysed using a weighted version of the Spearman correlation coefficient that allows for an analysis of the influence of parametric uncertainty on the derived correlations. We compared our results to those described by Merz *et al.* (2011), Osuch *et al.* (2015), and Sleziak *et al.* (2016). It seems that the dependence of model parameters on hydro-meteorological conditions varied between catchments, however, there are some similarities. In the Fuglebekken catchment, the FC parameter is positively correlated with the temperature related indices. Similar results were obtained for 273 catchments in Austria analysed by Merz *et al.* (2011), 213 catchments in Austria described by Sleziak *et al.* 2016, Axe and Kamp catchments (Osuch *et al.* 2015), but not Wieprz and Wimmera in the same study. The interpretation of this is that higher temperature or PET will lead to a larger soil moisture storage capacity. Depending on snow cover occurrence in the catchment, the correlation between FC parameter with precipitation related indices is characterised by the different sign of correlation or lack of statistically significant correlation. The FC parameter is also negatively correlated with precipitation related indices. The more intense or longer the precipitation, the lower soil moisture storage capacity, in simple words: with an increase of soil moisture, less water can be stored. Osuch *et al.* (2015) found that the correlation between FC and precipitation related indices is visible only in catchments without snow cover due to the processes of snow accumulation that causes significant delays in the transition between precipitation and flow, resulting in the positive correlation between the FC parameter and precipitation related indices. Such a relationship is also visible in the Fuglebekken catchment despite the fact that it is an Arctic catchment with snow cover. An explanation of this derogation from the rule is that our study period lasts from mid-June to the end of September – the period with almost no snow – therefore the Fuglebekken catchment in the study period behaves like a basin without snow cover.

The application of the HBV model allowed for the estimation of water storage in five conceptual storages: snow pack, melt water, soil moisture, fast runoff, and slow runoff. The results of analyses for the four study seasons indicated large changes in the slow runoff reservoir. Water storage in this reservoir varies greatly with time. It is likely that these changes are related to changes in ALT and also precipitation amount. We suspect that higher correlation coefficients between *KS* and ALT would be obtained if exact measurements of ALT were available. The simulations of ALT were characterised by large uncertainties, which make the analysis of such dependencies much more complicated. To our knowledge, this is the first study that attempted to reproduce subsurface field conditions in an unglaciated catchment in Svalbard that includes groundwater storage, runoff, and influence of changes in mean ALT. However, it should be highlighted that such outcomes were obtained using the HBV model only, without *in situ* measurements of ground water levels. Further investigations are planned as well as the application of other hydrological models.

Increased research has focused recently on runoff generation mechanisms in cold-region environments, where warming of air temperature over the last three decades have resulted in permafrost degradation and led to an increased storage capacity of affected soils and a higher contribution of groundwater to river discharge. Similar studies were performed in other catchments in the Arctic region, for example in: Russia (Semenova *et al.* 2013; Streletskiy *et al.* 2015), Canada (Carey & Woo 2001; Janowicz *et al.* 2004; Quinton *et al.* 2005; Quinton & Baltzer 2013), and Alaska (McNamara *et al.* 1998; Petrone *et al.* 2006). Projections show further degradation of permafrost, which will provide greater reservoirs for subsurface storage of groundwater and soil moisture (Bring *et al.* 2016). Such changes would result in changes in hydrological processes such as the pathway, timing, and amount of runoff, therefore links between terrestrial hydrology and other components of the Arctic freshwater system become very important, especially in changing environmental conditions.

## ACKNOWLEDGEMENTS

Financial support for this work was provided by the Polish National Science Centre through grant No. 2016/21/B/ST10/02509. This work was also partially supported by the Institute of Geophysics, Polish Academy of Sciences within statutory activities No 3841/E-41/S/2018 of the Ministry of Science and Higher Education of Poland, and from the funds of the Leading National Research Centre (KNOW) received by the Centre for Polar Studies for the period 2014–18. Meteorological data were provided by the Institute of Geophysics, Polish Academy of Sciences from the Polish Polar Station Hornsund.