2.1 Model framework

The analytical modeling framework of Botter et al. (2007c) for probabilistic characterization of rainfall-driven daily discharges is based on a soil moisture model originally proposed by Rodriguez-Iturbe et al. (1999). This point-scale model represents the dynamics of soil moisture as the result of a deterministic, state-dependent loss function, combined with stochastic increments triggered by rainfall events. Rodriguez-Iturbe et al. (1999) showed that the corresponding spatially averaged soil moisture s(t) can be obtained from the water balance equation as follows:

where −ρ[s(t)] is the loss function, due to evapotranspiration, surface runoff and deep percolation, and where ξ_t represents the stochastic instantaneous increments due to infiltration from rainfall.

Botter et al. (2007c) proposed describing the dynamics of daily streamflow with a similar stochastic differential equation, supposing that rainfall acts as a stochastic forcing for discharge production and that, at the catchment scale, the water is released following a linear decay:

where Q is the daily streamflow, k is the inverse of the time constant associated with the loss function and ${\mathit{\xi }}_{\text{t}}^{\prime \prime }$ is the stochastic process associated with discharge-producing precipitation events (i.e., the sequence of events that trigger a flow response in the river).

Figure 2Annual cycle of discharge and air temperature for three selected catchments representing three different hydrologic regimes (pluvial, snow-dominated and glacier). The mean monthly values computed over the entire observation period for each catchment are shown (see Table 1).

Download

It is assumed here that discharge Q is the result of a series of rainfall inputs (${\mathit{\xi }}_{\text{t}}^{\prime \prime }$) that deliver enough water to fill the water deficit in the soil, i.e., that deliver enough water to raise the soil moisture level above its retention capacity. The excess water is removed from the soil as subsurface runoff and becomes river discharge. This description of the streamflow response neglects any direct surface flow.

The overall rainfall forcing ξ is modeled as a marked Poisson process with frequency λ_p and exponentially distributed rainfall depths with average α (average rainfall on rain days). Nevertheless, not all the rainfall events trigger streamflow responses, and the discharge-producing precipitation events ξ” are also modeled as a marked Poisson process with a decreased frequency λ<λ_p and the same α. λ is influenced by the soil storage capacity and soil drying time and can be written as (Botter et al., 2007a)

where Γ(a,b) is a lower incomplete Gamma function with parameters a and b, $\mathit{\eta }=E/\left(n{Z}_{\text{r}}\right(s_{\mathrm{1}}-s_{\text{w}}\left)\right)$, $\mathit{\gamma }={\mathit{\gamma }}_{\text{p}}n{Z}_{\text{r}}(s_{\mathrm{1}}-s_{\text{w}})$ and ${\mathit{\gamma }}_{\text{p}}=\mathrm{1}/\mathit{\alpha }$. E is the maximum evapotranspiration rate and nZ_r(s₁−s_w) synthesizes the soil volume liable to be filled by water before drainage starts; n is the porosity of the soil, Z_r is the effective soil depth, s₁ is the retention capacity and s_w the permanent wilting point.

Figure 3Location of the case study catchments in Switzerland. The six biogeographical regions of Switzerland (Federal Office for the Environment, 2004) are summarized here into three main regions. Data source as follows: digital elevation model, SwissTopo (2005); catchments, Helbling (2016).

Download

As discussed in detail by Botter et al. (2007c), this framework results in the following probability distribution of daily discharges at the catchment scale:

where A is the catchment area. This corresponds to a Gamma distribution with shape parameter λ∕k and a scale parameter αkA. The corresponding expected mean discharge equals $\stackrel{\mathrm{‾}}{Q}=\mathit{\lambda }\mathit{\alpha }$. The model is suitable for steady-state conditions, at the annual or seasonal scale, depending on the temporal variability of the model parameters (Botter et al., 2007a).

Nonlinear storage–discharge relationships at the catchment scale are commonly observed (Botter et al., 2009; Brutsaert and Nieber, 1977; Mutzner et al., 2013). Accordingly, Botter et al. (2009) proposed an extension of the above modeling framework assuming that

where k_n and a are the constants of the nonlinear recession. As for the linear model, it is possible to obtain an equation for the pdf of the daily discharges:

where C is a normalizing constant (Botter et al., 2009).

Figure 4Examples of the temporal variation of the model parameters over the course of a year. The parameters are calculated for 90-day intervals beginning at the calendar day for which the value is plotted; for a given time window, the data points corresponding to this window in all available civil years are pooled together. (a–c) Residence time τ_k and mean daily precipitation depth α. (d–f) Precipitation frequency λ_p and discharge-producing frequency λ.

Download

2.2 Parameter estimation 1: forward estimation

We use the term “forward parameter estimation” to emphasize that the parameters are estimated directly from observed data, without calibration. This method is generally used in the context of this modeling framework for the estimation of the parameters related to the stochastic inputs (λ_p, α, λ), and this method is always used for these parameters in the present paper. However, the recession parameters (k, k_n and a) are either estimated in a forward mode (Botter et al., 2007c, 2009; Ceola et al., 2010; Schaefli et al., 2013) or in an inverse mode (Ceola et al., 2010) (see Sect. 2.3).

The computation of the precipitation parameters first involves the computation of a reference catchment-scale precipitation time series (here obtained from gridded data; see Sect. 3). Then interception losses (I) are subtracted from the observed daily precipitation depths. These losses are in fact evaporated (or sublimated in case of snow) before participating to soil moisture dynamics. Following Rodriguez-Iturbe et al. (1999), previous model applications generally assumed that these losses are accounted for when the frequency of precipitation events is corrected to the frequency of discharge-producing events. In view of understanding how the model parameters vary in space, it was decided here to treat interception losses explicitly with minimal assumptions about this process: different maximum interception depths are attributed to four different land covers: 4 mm for forests, 2 mm for low vegetation, 1 mm for impervious areas and 0 mm for water bodies (Gerrits, 2010). The catchment-scale maximum interception depth is obtained as the land-use-weighted average of these values, but a minimum interception depth of 1 mm is imposed. This catchment-scale interception depth is subtracted from daily precipitation depths, assuming that at a daily time step, all intercepted water re-evaporates during the same time step.

Figure 5Temporal variation of the residence time (${\mathit{\tau }}_k=k^{-\mathrm{1}}$) for the 25 catchments. The temporal variation is obtained as in Fig. 4.

Download

Instead of correcting the frequency of precipitation events λ_p according to Eq. (3), the frequency of discharge-producing events λ is estimated directly from the theoretical relationship between the mean discharge and the precipitation parameters, $\stackrel{\mathrm{‾}}{Q}=\mathit{\lambda }\mathit{\alpha }$ (see Eq. 4). Replacing the mean modeled discharge $\stackrel{\mathrm{‾}}{Q}$ with the mean observed discharge $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{Q}}$, it follows that

Estimating λ from the above equation rather than directly from the soil properties as in Eq. (3) has been shown by Ceola et al. (2010) to provide much better results, and this method has been used by the majority of studies since then (e.g., Ceola et al., 2010; Botter et al., 2013; Basso et al., 2015).

The recession parameter for the linear model is calculated directly from observed daily discharge based on a classical Brutsaert–Nieber recession analysis (Brutsaert and Nieber, 1977; Biswal and Marani, 2010, 2014; Mutzner et al., 2013), considering, however, only discharges below a certain threshold, fixed to 95 %. The nonlinear recession parameters k_n and a are also obtained from a recession analysis, using the same discharge threshold via linear regression of the logarithm of ($-\text{d}Q/\text{d}t$) versus the logarithm of Q, where a is the slope and k_n the intercept.

Figure 6Difference between λ and λ_p as a function of mean catchment elevation.

Download

2.3 Parameter estimation 2: inverse estimation

To objectively compare the potential of the linear and the nonlinear model formulations to capture observed flow-duration curves, the recession parameters for the linear model (k) and for the nonlinear model (k_n,a) are also estimated in a classical inverse estimation mode where the model parameters are obtained by maximizing the likelihood function formulated for the model. For the linear model, the likelihood function is obtained from the model as follows:

where the probability p(.) is obtained from Eq. (4) and $\mathit{\theta }=[\mathit{\alpha },A,\mathit{\lambda }]$ is the parameter vector containing all parameters that are estimated directly from observed data (i.e., not maximized). For the nonlinear case, the likelihood is obtained analogously by replacing k with k_n and a and using p(.) from Eq. (6).

2.4 Model evaluation criteria

To objectively compare different models, we propose using the Kolmogorov–Smirnov distance between the cdfs corresponding to different models (Ceola et al., 2010; Schaefli et al., 2013), i.e., the maximum difference between the values of the empirical and the modeled cumulative distributions:

where $F\left(\stackrel{\mathrm{̃}}{Q}\right)$ corresponds to the empirical cumulative distribution of the discharges and F(Q) to the modeled cumulative distribution of the discharges. A good model should have a low c^KS value.

Figure 7Modeled cdfs with forward and inverse parameter estimation for the three selected catchments. The shaded area is located between the cdf envelopes and represents the natural variability of the daily discharges.

Download

This comparison of the cdfs overcomes an important limitation inherent in the comparison of analytic pdfs and empirical pdfs. In fact, the choice of the number of classes for the calculation of the empirical pdf from observed data (via a so-called frequency polygon; Naghettini, 2016) can change the shape of the empirical pdf. The problem does not arise for cdfs given their cumulative nature.

Since the nonlinear model formulation has an additional parameter, the linear and the nonlinear models are also compared based on the Akaike information criterion (Burnham and Anderson, 2004; Laio et al., 2009; Ceola et al., 2010):

where n is the number of parameters of the model and ln$\left(\widehat{\mathcal{L}}\right)$ is the logarithm of the maximum likelihood function obtained by maximizing Eq. (8). As for c^KS, a good model should have a low c^AIC value.

Based on the above criterion, we measure the relative performance increase from the linear to the nonlinear model as follows:

where ${c_{\text{n}}}^{\text{AIC}}$ is the Akaike criterion for the nonlinear model and $c_{\text{l}}^{\text{AIC}}$ for the linear model. Taking the opposite of the relative difference between the Akaike criteria ensures that a higher r^AIC value indicates a stronger performance increase (recall that the Akaike criterion is to be minimized).

Figure 8As Fig. 7 but for the nonlinear model.

Download

In addition, to assess the performance difference between different models, the obtained models are compared to the natural variability of the observed discharge cdfs. Therefore, an empirical long-term cdf is constructed, obtained by ranking the observed data in ascending order and dividing the rank numbers by the total sample size. Furthermore, to assess the natural yearly variability, individual cdfs are constructed for each summer season of each civil year (Vogel and Fennessey, 1994). From this collection of annual cdfs, envelopes are based on the maximum and minimum values of discharge for each probability class of the annual cdfs. A reliable model should yield a cdf contained between these curves and should be as close as possible to the long-term cdf.

Figure 9Performance of the linear model and nonlinear model as a function of mean catchment elevation. The performance measure shown, c^KS, is zero for a perfect model.

Download

4.1 Discharge regimes and parameter seasonality

To gain further insights into the hydrological processes underlying the different regimes, Fig. 4 shows the within-year variability of the model parameters obtained by estimating the parameters in forward mode for moving and overlapping 90-day windows. The precipitation parameters α and λ_p do not show strong seasonal patterns, except for a few catchments such as the Goldach River (Fig. 4a). For snow and glacier regimes, the frequency of discharge-producing events, λ, increases strongly at the beginning of spring (Fig. 4b and c), which indicates the release of water from snowmelt or ice melt.

The inverse of the linear recession coefficient $\mathit{\tau }=k^{-\mathrm{1}}$ shows a coherent annual cycle for all catchments, independent of the underlying discharge regime (Fig. 5). This seasonal pattern with consistently low τ values during summer for all catchments clearly justifies the choice of a common summer season (June, July, August) for all regimes. The amplitude of the annual cycle (the difference between high and low τ values) is stronger for snow or glacier regimes, which reflects the fact that in these regimes, parts of the catchment are effectively dormant during the winter (Schaefli et al., 2013).

Table 2Parameter values and performance indicators for all the catchments for summer with linear model and forward estimation, summer linear model and inverse estimation, summer nonlinear model and forward estimation, winter nonlinear model and inverse estimation and winter linear model and forward estimation. $\stackrel{\mathrm{‾}}{Q}$ stands for the mean observed discharge, P_s the mean total precipitation during summer, $\stackrel{\mathrm{‾}}{T_{\text{s}}}$ the mean temperature during summer, I for interception depth and c^KS for the Kolmogorov–Smirnov distance. The indices stand for f, forward estimation; i, inverse estimation; l, linear model; n, nonlinear model.

Download Print Version | Download XLSX

4.2 Linear model

All estimated parameters for both forward and inverse estimations are summarized in Table 2, together with the values of the performance indicator c^KS. It can be noted that for 11 catchments (i.e., Rein da Sumvitg, Dischmabach, Alpbach, Grosstalbach, Rhône à Gletsch, Massa, Verzasca, Riale di Calneggia, Krummbach, Poschiavino and Ova da Cluozza), λ exceeds λ_p, contradicting the original description of the model (Botter et al., 2007b), which states that the discharge-producing frequency λ is smaller than the precipitation frequency λ_p. Such an exceedance of λ over λ_p should only happen in catchments or seasons with an additional source of water (in addition to rainfall), which in the present case is snowmelt or ice melt.

The exceedance of λ over λ_p increases with mean catchment elevation (Fig. 6), the limit of λ=λ_p being at around 1500 m a.s.l. This important result is further discussed in Sect. 5.

The cdfs obtained from all estimated parameters are presented in Fig. 7 for the three example case studies. For the catchment with rainfall-driven discharges (GOL), it can be seen that the probabilities of occurrence of low flows are largely overestimated with forward estimation (Fig. 7a). This is a typical indication that the recession timescale is underestimated. The model values even exceed the envelopes that represent the natural variability of the discharges. In the presence of snow, the linear model in forward estimation mode tends to underestimate low flows, with satisfactory results for some cases, such as the Dischmabach (Fig. 7b).

Overall, there is a strong increasing trend of the linear model performance with mean catchment elevation (Fig. 9a). Despite this, the results of the linear model are not satisfactory for the forward estimation method for any of the regimes.

The inverse estimation of the model parameters improves the results significantly, but the c^KS performance indicator shows relatively high values and the curves are visually not accurate, especially for pluvial regimes. This suggests that the model with a linear discharge decay is overall not suitable for the studied catchments.

4.3 Nonlinear models

The results obtained from inverse parameter estimation for the nonlinear model are very good (Fig. 8, Table 2), and the nonlinear model outperforms the linear model for all catchments, both in terms of the KS performance and in terms of the Akaike criterion (Table 2). The relative model performance increase (as measured by r^AIC) shows furthermore a strong inverse trend with mean catchment elevation (Fig. 10), which results from the increasing performance of the linear model with increasing elevation (Fig. 9b).

It is noteworthy that the two catchments for which the performance increase in the nonlinear model over the linear models exceeds 20 % are the two catchments that have a strongly karst-influenced regime (Scheulte at Vicques and Venoge at Ecublens).

As for the linear model, the forward estimation mode gives worse results than the inverse estimation mode. For some catchments (i.e., Murg-Wängi, Gürbe, Sense, Ilfis, and Grosstalbach), the forward estimation mode gives nevertheless very good results with c^KS below 0.1. In general, for the catchments where the discrepancies between modeled and observed cdfs are due to an underestimation of τ, the nonlinear model yields a significant improvement. For catchments where the recession timescale is overestimated with the linear model, the nonlinear model in forward mode leads to a performance decrease.