Evaluating the effects of alternative model structures on dynamic storage simulation in heterogeneous boreal catchments

Storage and release of water are essential catchment functions, with large effects on modulating hydrological extremes such as drought and floods (Creutzfeldt et al. 2012). Water entering catchments is stored in different forms, including snow, ice, lakes, soil moisture, and groundwater (Riegger & Tourian 2014). Depending on soil physical properties, water can be retained in the unsaturated zone, move through the soil column and percolate to saturated layers as groundwater recharge. Accordingly, there are various pathways that water can take when it flows between unsaturated and saturated zones before leaving the catchment, depending on soil physical properties, catchment characteristics, and climate (Jutebring Sterte et al. 2021).

The temporal dynamics of stream discharge depend largely on catchment storage (Kirchner 2009; Creutzfeldt et al. 2012; Brauer et al. 2013). In periods when evapotranspiration and precipitation are negligible, it can be assumed that the only driver of discharge is the amount of water stored in the catchment (Moore 1997; Kirchner 2009). Water content freely available for drainage and flow is known as ‘mobile storage’ (Farrick & Branfireun 2014; Staudinger et al. 2017). However, not all catchment storage contributes to the runoff generation process, which is sometimes referred to as ‘hydraulically decoupled storage’ (Dralle et al. 2018). Changes in the amount of these hydrologically decoupled storages can have a negligible effect on runoff generation, at least in the short term, because of deep groundwater storage, or have no effect, such as canopy interception that leaves the basin through evaporation without entering the soil.

In many previous studies, catchment storage at different spatial scales has been quantified using a range of different methods: water balance approaches (Sayama et al. 2011; Wang & Alimohammadi 2012), hydrometric analysis (Kirchner 2009; McNamara et al. 2011; Brauer et al. 2013; Amvrosiadi et al. 2017a, 2017b), analysis of stable isotope tracers (Soulsby et al. 2009; Sprenger et al. 2018), gravimeter measurements (Creutzfeldt et al. 2012; Kennedy et al. 2014), and application of hydrological models (Staudinger et al. 2017; Sterte et al. 2018; Karlsen et al. 2019). One of the most common methods for direct estimation of catchment water storage involves using a network of piezometer wells and soil moisture probes to provide detailed information on whether a sufficient number of sensors has been installed (Kalbus et al. 2006). However, due to the high spatiotemporal variability, collecting massive hydrological data across a broad range of scales is time-consuming and costly (Kirchner 2009). Hence, efforts to reduce field investigations have resulted in various modeling attempts to understand better the complexities of the hydrological behavior in heterogeneous catchments.

There are several definitions of catchment water storage (Staudinger et al. 2017). Here, we focus on the part of the storage that is relevant for catchment discharge. This leads to an operational definition where the storage is zero when streamflow ceases completely. Various terms have been proposed for this discharge-generating storage; these include ‘dynamic storage’ (Kirchner 2009; Sayama et al. 2011; Staudinger et al. 2017), ‘active storage’ (McNamara et al. 2011), and ‘direct storage’ (Dralle et al. 2018). Here we use the term dynamic storage as a critical metric to assess the sensitivity of a catchment to extreme weather events such as flooding or its resistance to drought by maintaining low flows (McNamara et al. 2011). Revealing the dynamic storage-discharge relationship can also help make more accurate predictions of streamflow changes in response to global warming.

Bucket-type models, also called conceptual models, are suitable tools for simulating runoff and water availability based on storage–discharge relationships and represent in general terms how precipitation results in groundwater recharge, evapotranspiration, and discharge (Hrachowitz & Clark 2017; Sitterson et al. 2018). These models consist of several reservoirs, which can be linear or nonlinear and be connected in serial or parallel (Stoelzle et al. 2015; Parra et al. 2019a). The most important advantage of these bucket-type models, compared to more detailed and complex fully physically-based models, is a lower number of free parameters, which means that the parameter values of the models can, in principle, be estimated by calibration only. However, even for such models, model parameter uncertainty has to be considered when deriving storage estimates by modeling. Furthermore, there is model structure uncertainty as there are usually various possible alternatives to arrange the different buckets. For example, the proper choice of linear and nonlinear storage-discharge relationships has been debated for a long time (Stoelzle et al. 2015). Some researchers have concluded that groundwater release is a linear process (Chapman 1999; Fenicia et al. 2006), while others have demonstrated that nonlinear relationships are more appropriate (Wittenberg 1999; Mishra et al. 2003; Botter et al. 2009; Maneta et al. 2018; Rezaei-Sadr 2019). These findings show that choosing appropriate model structures is challenging and uncertain. In other words, the mathematical formulation of hydrological functions that control the transformation of rainfall to runoff might vary by type of catchment and climate, and that a ‘one size fits all’ is not easily found (Hogue et al. 2006; Ajami et al. 2007).

Many studies have used bucket-type models to improve the understanding of hydrological functioning in catchments, such as estimation of water storage (Krasnostein & Oldham 2004; Mendoza-Sanchez et al. 2013; Staudinger et al. 2017) as well as streamflow signatures (Hailegeorgis & Alfredsen 2015; Ledesma & Futter 2017; Teutschbein et al. 2018). These studies have also demonstrated that bucket models can provide suitable representations of catchment hydrology and result in good runoff simulations for both calibration and validation periods. However, as pointed out in previous studies, achieving an acceptable value of model performance based on a good agreement between modeled and observed streamflow data does not necessarily correspond to a good simulation for other hydrological functioning (Gupta et al. 2012; Teutschbein et al. 2015; Lane et al. 2019; Seibert et al. 2019). Therefore, choosing a model structure with an appropriate degree of complexity is a crucial step in catchment modeling.

Hence, despite a good model performance for a certain calibration or validation period, model outputs can still be unreliable when extrapolated beyond calibration conditions (Seibert 1997). Therefore, a proper analysis of how the model parameters and the model structures are correlated and a detailed assessment of their effects on the model output can improve the model reliability. To our knowledge, few previous studies have comprehensively analyzed the most important parameters controlling the amount of storage in the catchment reservoirs (Tetzlaff et al. 2015; Teutschbein et al. 2015; Ledesma & Futter 2017). This is especially true in northern latitude catchments defined by long-lasting winters and snowmelt-dominated hydrological conditions.

The main purpose of the approach lies in the evaluation of the performance of different model structures, which only differ in the number of storage reservoirs, in estimating dynamic storage using a bucket-type model for a large number of catchments with contrasting landscape characteristics but similar climatic conditions. In addition, we aimed to better understand the relationship between dominant catchment characteristics and dynamic storage variability across a heterogeneous boreal landscape. We addressed the following research questions: (1) How does the uncertainty of parameters in the response routine of a bucket-type model affect the dynamic storage estimation? (2) Is there a major difference between dynamic storage estimates obtained using different model structures? If yes, which model structure is most reliable at different catchment scales? We then used this approach to evaluate how simulated dynamic storage varies in both time and space for different subcatchments. Furthermore, we studied potential relationships between the different types of storage zones and dominant catchment characteristics.

The study was conducted in the Krycklan catchment located in northern Sweden, about 50 km west of the Baltic Sea coast (64°25′ N, 19°46′ E), with an area of approximately 68 km². The catchment area comprises 14 nested subcatchments with elevation ranging from 130 to 370 meters above sea level, as shown in Figure 1. The higher elevations of the catchment are dominated by till (58% of total area) and peat soil, and lower elevations are covered by sediment soils (Laudon et al. 2021). Forests cover 87% of the entire catchment and are dominated by Scots pine (Pinus sylvestris) and Norway spruce (Picea abies).

The hydrological characteristics of the catchment are well studied and are monitored by 14 streamflow gauging stations (named C1-C20, with some of the original stations abandoned for various reasons) (Table 1). The climate in this area is characterized by relatively cold winters with consistent snow coverage for four to six months of the year, between November to April. The annual average precipitation was 614 mm for 1981–2010, and the average daily mean temperature was 1.8 °C. The mean temperature was −9.5 °C in January and +14.7 °C in July (Laudon et al. 2013).

Discharge observations from 14 subcatchments within the Krycklan catchment were used (Karlsen et al. 2019). Water levels, measured using automatic stage loggers, were possible year-round for five gauging stations in heated houses (C2 and C4 have been heated since 2011, C5 and C13 since 2012, and C7 since 1981). The monitoring program began in 1981 inside a heated hut (with an hourly resolution), and in 2003, gauging started in all 14 subcatchments. The lengths of all hydro-climatic data series in our study were adjusted to match the span of subcatchments with the least amount of data, covering a period of seven hydrologic years (1 October–30 September) from 2011 to 2017. Frequent manual water-level measurements (monthly during winter, minimum bi-weekly during the rest of the year) were made to calibrate automatic water level data, and stage-discharge relationships were defined using manual flow gauging (Karlsen et al. 2019). Specific discharge, defined as streamflow per unit drainage area, was calculated for each catchment.

Air temperature, humidity, net radiation, and wind speed for estimating potential evapotranspiration (PET) via the Penman equation were measured in the central part of the Krycklan catchment at the Svartberget research station (Laudon et al. 2021). All climate data were assumed to be uniform for the whole catchment area. Precipitation from one centrally placed weather station was used (64°14′ N, 19°46′ E, 225 m a.s.l) for all catchments as there is no significant elevation gradient observed for the region by the network of the Swedish Meteorological and Hydrological Institute (SMHI) (Karlsen et al. 2016).

A LiDAR DEM with 2 m resolution was created from a point cloud with a point density of 15–25 points/square meter. This DEM was hydrologically corrected by burning streams and culverts across roads as described in Lidberg et al. 2017. Catchment areas were delineated from the hydrologically correct DEM using Deterministic-8 (D8) (O'Callaghan & Mark 1984). The Swedish property map (1:12,500, Lantmäteriet Gävle, Sweden) was used to calculate forest, lake and, wetland coverage for each catchment. The proportion of soil type cover was calculated for sediment soils, till, and thin soils using the quaternary deposits map (1:100,000, Geological Survey of Sweden, Uppsala, Sweden). Landscape characteristics of all catchments, including the main outlet (C16, in bold), are summarized in Table 2.

A modified version of a bucket-type, semi-distributed hydrological model, namely the HBV model (Lindström et al. 1997) in the version HBV light (Seibert & Vis 2012) was used to simulate dynamic storage at a daily time step in the catchment. The HBV model is a widely used bucket-type model for simulating runoff. One advantage of HBV is the limited demand for input data. The input data for the HBV model were rainfall, air temperature data, and derived PET. The observed streamflow data were used for the calibration of the model. In general, the model consists of four commonly used routines which represent: (1) snow by a degree-day method; (2) soil moisture, and (3) groundwater by three linear reservoir equations; and (4) channel routing by a triangular weighting function. In the snow routine, a threshold temperature is used to distinguish between rainfall and snowfall, and snowmelt is considered by a degree-day approach.

A more detailed description of HBV's routines can be found in Seibert & Vis (2012). HBV-light includes three alternative model structures based on a varying number of conceptual reservoirs (one, two, or three) and different shapes of the storage-discharge relationship (linear versus nonlinear).

We implemented the three model structures of the response routine. This included a one-bucket structure with a single groundwater reservoir and three linear outflows (Q0, Q1, and Q2) at three different thresholds; a two-bucket structure with two reservoirs in parallel with a nonlinear outflow (Q1) in an upper bucket and one linear outflow (Q2) in a lower bucket (basic version); and a three-bucket structure with three parallel reservoirs (STZ, SUZ, and SLZ) and one linear outflow in each reservoir (Figure 2). It should be noted that these model structures differ only in the response routine, i.e., the same snow and soil routines were used in all model variants.

The built-in Genetic Algorithm followed by a Powell optimization (GAP, Seibert 2000) was used to calibrate the models for each catchment. For calibration, warm-up periods of at least one year were used for each catchment based on data availability. The initially chosen possible parameter ranges were adapted from Uhlenbrook et al. (1999) where preliminary simulations indicated that suitable parameter values were close to the limits.

We considered parameterizations resulting in good runoff simulations in terms of Nash–Sutcliffe efficiencies (R_eff) as plausible. From the 100 calibration trials, we thus derived an ensemble of plausible simulations that result in a range of storage estimates for the upper (SUZ) and lower storage (SLZ) reservoirs. In the next step, the estimated storage from all reservoirs was combined and the difference between the minimum and maximum values was considered as dynamic storage.

A qualitative preparatory evaluation of parameter sensitivity indicated that automatic calibration could sometimes result in very small values for K2 (baseflow recession coefficient) with tiny model performance improvements for these very small K2 values (Figure 3 and Supplementary Material, Figure A). These very small K2 values, in turn, resulted in extremely large estimates of the derived catchment storage for some of the catchments. This optimization issue, which greatly affected the storage simulation results, occurred mainly with the three-bucket model structure. All but one of the catchments (C20) showed a strongly nonlinear relationship between K2 and storage in the lower zone. To avoid artefacts like this, we limited the lower range of K2 to 0.02. The summary statistics (min, max, and mean) of calibrated parameter values with limited K2 parameter values are summarized in Table 3.

To investigate the relationship between landscape characteristics and simulated dynamic storage, we used non-parametric Spearman rank correlations (Spearman 1904). The storage metrics used in this analysis were dynamic storage (the difference between the minimum and maximum storage), mean storage (mean value during the entire analysis period), storage in the upper and lower zone separately (SUZ and SLZ), and mean storage for each season. Instead of using a single model output, we used the average of all three model structures to calculate these storage metrics and then examined their relationship with landscape characteristics.

Based on the results, differences in model performance were found between catchments with different landscape characteristics. There are also performance differences between the different model structures. The range of the best 100 model performances (R_eff) of the three model structures for each subcatchment is shown in Figure 4. The plot suggests that model performance decreases as the mean elevation above the streams increases (r = −0.85, p < 0.05). Moreover, the three-bucket model gave good agreement between the simulated and observed streamflow in all the catchments (R_eff between 0.76 and 0.87). The one-bucket model performed well in the lake-influenced catchments, C5 and C6 (0.80 < R_eff < 0.87), while it performed worse for large catchments with deep soils, C14 and C16 (0.68 < R_eff < 0.75). This indicates that the one-bucket model could not represent the hydrological functioning as well as the three-bucket model in the large catchments with deep sediment soils. The one-bucket model resulted in the poorest simulations in C14 and the best simulation in C6 with maximum R_eff values of 0.74 and 0.87, respectively. Additionally, the two-bucket model with two storage reservoirs and only two outflows (intermediate and baseflow) performed similarly to the one-bucket model for C2, C4, C9, and C20. The highest model performances were obtained for C5, C6 and, C9, which have a high lake percentage and lower elevation above their stream network. The simulated discharge time series derived from the three model structures calibrated against observed discharge, for four representing catchments including C4, C5, C7, and C16, are shown in supporting information Figure B1–4. These catchments are also representing the highest (C4 and C5), average (C7), and lowest NS values (C16). We also provided a detailed analysis of snow water equivalent and soil moisture estimated using the three model structures for the above-mentioned catchments. Results showed that all model structures simulated snow water equivalent relatively similar for the four representing catchments, though C4 and C5 had a slightly higher amount (supporting information Figure C1-4).

The contribution (as a fraction of the total discharge) of outflow from the surface (Q0), upper (Q1), and the lower (Q2) flow paths for all the model structures are shown in Figure 5 (note that the standard version of HBV has only two outflows, Q1 and Q2). The boxplots are based on 100 streamflow characteristic values derived from 100 different parameter sets. When comparing the runoff components calculated by all three model structures, remarkable differences were seen between different subcatchments and the three HBV model structures. For example, the difference in the contribution of different Q components between the sediment and till-dominated catchments was apparent. It can be seen in Figure 5 that the calibrated contribution of discharge from the lower groundwater bucket (baseflow) was higher in larger catchments with more sediment soils (e.g., C14, C15, and C16). In contrast, in catchments C5 and C4, the fast flow (Q1) showed a more significant contribution to total runoff generation.

Simulation results from the three-bucket structure, which showed reasonably good performance across all catchments, demonstrate differences in storage components between catchments. The simulated maximum storage in the upper and lower groundwater bucket are plotted in Figure 6. It can first be noted that in wetland-dominated catchments, the maximum storage in the upper groundwater zone was higher than the other subcatchments. In contrast, catchments dominated by sediment soils retained a higher amount of water in their lower groundwater bucket. It should also be noted that the upper reservoir (SUZ) is nearly empty during long periods and the mean value is, therefore, small. The simulated daily storage in upper and lower reservoirs estimated by each model structure was then aggregated to compare changes in dynamic storage estimated by different model structures.

The dynamic storage estimates calculated as the difference between the maximum and minimum groundwater storages from the different model structures are presented in Figure 7. While the estimated storage volumes varied, a consistent overall pattern among the catchments could be observed. In general, the forested catchments had lower storage estimates than the catchments with a lake or large wetlands, and storage estimates were even larger for the sediment catchments. For some of the catchments (C5, C6, and C9), the estimated dynamic storages were similar for the different model structures. In contrast, the differences were higher for others (C14, C15, and C16).

In contrast to the overall range of dynamic storage, the mean storage value calculated by each model structure showed large variation among catchments (Figure 8, Table 4).

For most of the catchments there was only a slight increase in estimated dynamic storage from the one-bucket to the two-bucket structure, but a significant increase from the two-bucket to the three-bucket in most of the catchments (Figure 8). This pattern was opposite for the lake dominated-catchment (C5).

Simulated dynamic storage, however, exhibited more variation both within and among the catchments on shorter timescales (Figure 9). The variation in dynamic storage between catchments was highest in spring with a minimum value of 9 mm for C1 and a maximum value of 30.9 mm for C20, and lowest in winter with a range of 5 mm. According to the three-bucket model, C5 had the largest seasonal variation with the mean storage approximately 5 times greater in spring compared to winter. Additionally, we quantified the similarity degree between the mean seasonal storage estimates by the three model structures using the Spearman correlation coefficient (Table 5). The comparison of similarity metrics revealed that the correlation between the mean storages estimated by the three model structures was highest for spring and lowest for winter. Moreover, the one-bucket and three-bucket structures had the largest differences in estimating dynamic storage during all seasons. These differences were highest for C14 and lowest for C5.

The results of Spearman's correlation analysis between catchment storage and dominant landscape characteristics are summarized in Table 6. The test showed a positive relationship between catchment tree volume, till soils, and total water storage, ΔV. This indicates that as the amount of tree volume and till soils increase, the total water storage, ΔV also increases (r = 0.60, p < 0.05). In contrast, with increasing wetland percentage, the amount of total water storage, ΔV decreased (r = −0.59, p < 0.05). Dynamic storage (HBV), on the other hand, showed a negative correlation with tree volume (r = −0.77, p < 0.05) and till soil (r = −0.60, p < 0.05). It also showed a positive correlation with lake percentage. Positive correlations were also found between catchment area, mean storage during the entire period, and mean seasonal storage (winter, spring, and summer). Elevation above stream (EAS) and slope both showed a negative relationship with storage in the upper zone (SUZ) and a positive relationship with storage in the lower zone (SLZ). In contrast, there was no correlation between mean catchment elevation and any of the storage metrics.

Catchment area and sediment soil were positively correlated to mean SLZ, which means that catchments with a larger area, higher EAS, larger amounts of sediment soil, and steeper slope would have higher storage in the lower groundwater bucket. For the mean storage over the whole study period, a positive relationship was observed with the mean catchment area and the proportion of sediment soil. However, the mean storage was inversely affected by till soils and tree volume (Figure 10(c) and 10(d)). For the mean seasonal storage, an inverse correlation was found between till soils and the mean storage in autumn, spring, and winter, with a much weaker correlation in summer (r = −0.388, p < 0.05).

Tree volume also revealed a negative correlation with all seasonal storage, although this relationship was stronger in winter and summer. Furthermore, a positive correlation was observed between sediment soil and the mean storage in summer and winter. This highlighted the role of sediment soils in maintaining a relatively high flow, mainly during the dry periods. The mean slope showed a significant correlation only to the mean storage during summer. Contrary to what was expected, the mean catchment soil depth did not correlate with any storage metrics.

The analysis of model performances showed that the three-bucket model captures the hydrological behavior best in all catchments. This finding is similar to Stoelzle et al. (2015), who suggested the simple structures of one or two reservoirs are less efficient in simulating catchment baseflow. Fenicia et al. (2014) also found that single-reservoir models are too simplistic to identify the different runoff-generating processes and that models with multiple parallel flow paths performed better. In addition, for all model structures, the best model performances were achieved in catchments with lower slope and mean elevation above streams (EAS) (Figure 4). A possible explanation is that with increasing elevations above the stream, the spatial heterogeneity of flow paths and processes within the respective catchments increases. Therefore, models that treat catchments as a single unit might have difficulties in representing the hydrologic functioning of these catchments.

Our results also suggest that in large catchments with higher slope and EAS, the one-bucket and two-bucket reservoirs were not adequate, and better performance was achieved by multiple reservoirs and flow pathways. The underlying reason could be that for the larger catchments there also is a higher percentage of sediment cover and deeper soils, which both can increase spatial variability within a catchment and, as a result, affect model performance. Furthermore, the relatively poor performance of all three model structures in large sediment catchments suggests that capturing the storage–discharge functioning is more challenging than the other catchments, likely because of more complex hydrological responses (Jutebring-Sterte et al. 2021). However, our results contradict the findings of other studies. For example, Broderick et al. (2016) used the NSE criterion and found higher model performance for catchments with greater storage capacity and baseflow as they are less sensitive to storm events, and therefore produce a less variable flow series. Additionally, Van Esse et al. (2013) tested 12 different conceptual model structures on 237 catchments in France and found that conceptual models have higher efficiency in larger catchments.

The analysis of simulated runoff components suggests various dominant hydrological functioning in different catchments that can potentially be attributed to the spatial heterogeneity of soil properties in the studied catchment. The results of 100 calibration trials showed the importance of baseflow in sustaining runoff generation, especially in sediment soil catchments. This is in line with findings from Karlsen et al. (2019), Teutschbein et al. (2015), and Peralta-Tapia et al. (2015) that also demonstrated that the amount of baseflow draining into the streams has a positive correlation with increasing catchment area and sediment cover.

Conversely, in wetland-dominated catchments (C4 and C5), overland flow from peat soils dominates storm runoff generation. Laudon et al. (2007) and Seibert et al. (2003) have also pointed out that in peatlands, groundwater levels reach the soil surface and, combined with overland flow from snowmelt, results in higher specific discharge compared to other catchments. In contrast, larger catchments with deep sediment soils have deeper flow paths and higher contribution to baseflow.

We also illustrated and compared the differences in the simulated storage components and flow paths among the catchment categories. We observe differences in dynamic storage at various locations (upper and lower storage zones) within the catchment between different landscape classes. Since the studied catchments are partially nested and have a similar climatological regime, the variability in storage behavior should be mainly due to different physical attributes such as soils, landform, and vegetation cover.

The large sediment catchments in the lower part of Krycklan had a higher infiltration rate and deeper subsurface flow paths. This suggests that water is stored mainly in the lower storage zone and thus can help maintain baseflow during dry periods. Our results are consistent with those obtained by Peralta-Tapia et al. (2015), who used stable isotopes to show that the contribution of deeper groundwater flow paths to catchment discharge increases with the catchment size. Conversely, the wetland-dominated catchments had more storage volume in the upper groundwater zone during high flows. This supports the finding of Laudon et al. (2007) that in wetland-dominated catchments, the proportion of new or event water was much higher than forest on till catchments. Using isotopic data, they also concluded that there is a shallow pathway in wetland catchments close to the surface, which might be formed by a concrete frost layer inhibiting the infiltration of rain and meltwater.

The different model structures had a varying number of storage reservoirs, which could partly explain differences in model performance as they capture hydrological catchment functioning differently. However, our main question was how much the estimated storages differed between these structures. As expected, using multiple storage reservoirs, even with the same number of parameters (one-bucket and three-bucket structures), makes a large difference in some catchments. Typically, for the most lake-influenced catchments (C5, C6, and C9), as well as the sediment dominated catchment C20, the modeling results were similar for all model structures, which implies that regardless of the high model performance value, model structures with only one reservoir can adequately represent the dynamic storage behavior. The explanation for C20, which was an outlier of sediment categories, could be due to its smaller drainage area compared to other sediment-dominated catchments (C14 and C16).

Although applying different model structures resulted in different estimates of dynamic storage for each catchment, all dynamic storage estimates calculated for different model structures ranked similar among the catchments. For instance, all model structures resulted in less dynamic storage for till catchments, while a higher amount of dynamic storage was obtained for those most influenced by lakes. These results are in line with Karlsen et al. (2019), who calculated dynamic storage, as the difference between the observed daily minimum and maximum specific discharge for each catchment and found the lowest dynamic storage for C2 (20 mm), and the highest for C6 (75 mm).

Moreover, we found that the spatial variability of dynamic storage is higher in spring, possibly due to differences in soil frost extent during snowmelt, depending on the vegetation cover and soil type. The difference between the two dominant soil groups (till and sediment) in transferring snowmelt water inflow within the catchment may also account for this difference. Catchments with sediment soil have deeper, or probably more groundwater paths, while catchments with till soils are less permeable so that most of the snowmelt leaves the catchment as overland flow and through the upper storage box.

As mentioned above, the boreal ecosystems are characterized by long winter and large snow accumulation. Therefore, spring snowmelt has a dominant contribution to the annual storage-discharge magnitude in northern regions (Laudon & Ottosson Löfvenius 2016). The simulation results determined that although the precipitation and temperature were similar across all catchments, the differences in catchment topography, soil, and land cover resulted in differences in the snow accumulation among catchments. Compared to other forested similar size catchments, the higher amount of discharge and dynamic storage of lake and wetland-dominated catchments (C5 and C4) during spring snowmelt can be explained by the larger snowpacks in these open canopy catchments (Kozii et al. 2017). In seasons when rainfall has the dominant effect on runoff, the lake catchment, C5, had the smallest amount of dynamic storage compared to other catchments. This agrees with previous findings suggesting that deep glacial till soils, wetlands, and forests on sediment soils have more storage during low flow conditions, while shallower till soils and open-water wetlands have been shown to sustain less water as they already contain a lot of water and their storage capacity (active storage) is small (Meriö et al. 2019).

The three different model structures showed that although the pattern of dynamic storage (the difference between the minimum and maximum annual storage values) estimates among catchments were similar, the catchments ranking for mean storage values differed greatly. However, for every model structure, catchments with higher tree volume and more till soils had a lower amount of dynamic storage compared to sediment catchments with less tree volume. This is likely due to transpiration and interception losses that reduce groundwater recharge (Ilstedt et al. 2016; Bonnesoeur et al. 2019), or soil water uptake within the tree root zone (Allen & Chapman 2001).

To consider parameter uncertainty and to ensure that results are not affected by it, as a first step we calculated the conceptual dynamic storage using an ensemble average taken over 100 best simulations, to have more robust storage estimations. Second, the results revealed that although the three-bucket model performed relatively better for all catchments, for some catchments the one-bucket model with less complexity also yielded a high model performance. Additionally, our findings indicated that the three-bucket model, despite its high performance, had more uncertainty in simulating dynamic storage than the one-bucket model, which gained the lowest model performance. Most studies, on the other hand, have focused on the impacts of model structures on discharge simulation (Uhlenbrook et al. 1999; Van Esse et al. 2013; Fenicia et al. 2014; Parra et al. 2019a), and considered the performance criteria to determine the predictive power of a model. With our focus on estimating a fraction of total storage that is actively controlling discharge release, and using it as a comparison of hydrological functioning across contrasting boreal catchments, makes the validation more difficult.

Therefore, we argue that using a simple ensemble averaging method that combines the prediction of each model equally would decrease the weakness of every single model in representing the hydrological functioning. This approach will also reduce the biases resulting from modelers' personal preferences

The results from the Spearman rank correlation tests suggest that the largest mean storage estimates were found in larger catchments, and catchments characterized by a large proportion of sedimentary soil, especially during winter. This finding is in accordance with other studies in Krycklan. For example, Karlsen et al. (2016), who applied partial least square regression to quantify the linkage between recession characteristics and catchment properties found that with increasing catchment size, the relative contribution of groundwater, especially during winter baseflow, increased. This is also in agreement with Tiwari et al. (2014) and Peralta-Tapia et al. (2016), who have shown that deep groundwater contribution in the catchment increases with catchment size.

When comparing our findings with Jutebring Sterte et al. (2021), a similar pattern was seen between travel time and dynamic storage among the catchments. For example, the modeled mean travel time (MTT) increased with increasing slope, sediment soil proportion, and catchment size. In their study, C20 was seen as an outlier in most of the correlation analyses, and the longest MTTs were found for that catchment despite its relatively small size. A direct association of sedimentary soil on prolonged groundwater storage-release processes has also been found in a study conducted in south-central Chile (Parra et al. 2019b), while Staudinger et al. (2017) found no correlation between catchment area and any of the storage metrics in alpine catchments.

We also found a strong correlation between catchment tree volume and dynamic storage, which is in line with findings of Barrientos & Iroumé (2018), who explored the impacts of forest management on dynamic storage in 15 forested catchments and concluded that catchments with lower forest cover (biomass volume and plantation density) have higher storage volumes. This might be due to the increasing evapotranspiration rate by deep-rooted trees with water uptake in sub-surface storage and therefore results in reducing the available amounts of dynamic storage. Furthermore, we assumed that the positive association between the total water storage (estimated by the water balance method) and till soils could be because densely forested catchments store significant amounts of input water, and this large quantity of immobile storage was included in the calculations of the water balance method.

In this study, we used the HBV model with three model structures, each with a different number of storage reservoirs, to evaluate differences in dynamic storage simulation. First, we found that there is a high variability among different parameterizations. This means that to achieve robust results for further analyses, it is important to consider parameter uncertainty, for instance, by an approach such as the 100 calibration trials used in our study. Second, we found high variability in dynamic storage estimates not only between the catchments but also between the different model structures.

Moreover, this study highlighted the importance of a three-bucket conceptual model that generates runoff using multiple storage reservoirs, especially in large groundwater-dominated catchments with deeper sediment soils. In contrast, for the lake-influenced catchments, a single-reservoir structure performed equally well, indicating that these simple structures can adequately represent the hydrological functioning in such catchments due to the lower storage capacity and shallow groundwater levels.

We also concluded that an ensemble average of simulations can be used as a point of hydrological functioning comparison between different catchments. Considering that the model validation using real data in these large-scale heterogeneous catchments is not possible in practice, this can reduce the time, cost, and risk associated with uncertainty of each single model structure.

Finally, as demonstrated by our results, the dynamic catchment storage could be explained by landscape features such as drainage area, slope, soil characteristics, and tree volume.

The GIS and hydroclimate data are publicly available from (https://www.slu.se/Krycklan). The HBV-light software and its tutorial are freely available from (https://www.geo.uzh.ch/en/units/h2k/Services/HBV-Model.html). The source code of the different response routine model structures can be requested from the second author.

The authors declare that they have no conflict of interest.

The GIS and hydroclimate data are publicly available from https://www.slu.se/Krycklan. The HBV-light software and its tutorial are freely available from https://www.geo.uzh.ch/en/units/h2k/Services/HBV-Model.html. The source code of the different response routine model structures can be requested from the second author.