Integrated sensitivity analysis of a macroscale hydrologic model in the north of the Iberian Peninsula
Patricio Yeste, Matilde García-Valdecasas Ojeda, Sonia R. Gámiz-Fortis, Yolanda Castro-Díez, María Jesús Esteban-Parra
11 Integrated Sensitivity Analysis of a Macroscale Hydrologic Model in the North of the Iberian Peninsula
Patricio Yeste (1) , Matilde García-Valdecasas Ojeda (1) , Sonia R. Gámiz-Fortis (1) , Yolanda Castro-Díez (1) and María Jesús Esteban-Parra (1) (1)
Dept. Applied Physics, University of Granada, Spain
Corresponding author:
María Jesús Esteban-Parra ([email protected])
Abstract
Process-based hydrologic models allow to identify the behavior of a basin providing a mathematical description of the hydrologic processes underlying the runoff mechanisms that govern the streamflow generation. This study focuses on a macroscale application of the Variable Infiltration Capacity (VIC) model over 31 headwater subwatersheds belonging to the Duero River Basin, located in the Iberian Peninsula, through a three-part approach: (1) the calibration and validation of the VIC model for all the subwatersheds; (2) an integrated sensitivity analysis concerning the soil parameters chosen for the calibration, and (3) an assessment of equifinality and the efficiency of the calibration algorithm. The calibration and validation processes showed good results for most of the subwatersheds in a computationally efficient way using the Shuffled-Complex-Evolution algorithm. The sensitivity measures were obtained with the Standardized Regression Coefficients method through a post-process of the outputs of a Monte Carlo simulation carried out for 10 000 parameter samples for each subwatershed. This allowed to quantify the sensitivity of the water balance components
This is a preprint of a manuscript accepted for publication in Journal of Hydrology.The published article is available at https://doi.org/10.1016/j.jhydrol.2020.125230 to the selected parameters for the calibration and understanding the strong dependencies between them. The final assessment of the equifinality hypothesis manifested that there are many parameter samples with performances as good as the optimum, calculated using the calibration algorithm. For almost all the analyzed subwatersheds the calibration algorithm resulted efficient, reaching the optimal fit. Both the Monte Carlo simulation and the use of a calibration algorithm will be of interest for other feasible applications of the VIC model in other river basins.
Keywords
Duero River Basin; VIC model; calibration; validation; sensitivity analysis; equifinality Introduction
Water resources in the Mediterranean Basin have undergone dramatic changes during the 20th century as a consequence of the rising temperatures and the significant decrease of precipitation (García-Ruiz et al., 2011). The effects of climate change in this region are already noticeable and are expected to be much more pronounced by the end of the 21th century (IPCC, 2014). This fact, together with the increasing water demand for agriculture, industry and urban supply, makes the water scarcity problem of paramount importance, being its accurate identification essential for adopting adequate water management strategies and mitigation measures that ensure the sustainability of the water resources (Chavez-Jimenez et al., 2013; Garrote et al., 2016). As a part of the Mediterranean region, the Iberian Peninsula conforms a vulnerable area that has been identified as a hotspot (Diffenbaugh and Giorgi, 2012) where the streamflows have shown a marked reduction during the last half century (Lorenzo-Lacruz et al., 2012, 2013). Being able to identify the hydrologic behavior of a basin is necessary in order to assess the effects of climate and land changing conditions, and therefore a profound description of the main hydrologic processes governing the response of the basin is required. In this way, process-based hydrologic models are powerful tools that represent the underlying runoff mechanisms governing the streamflow generation for a given basin, and therefore constitute mathematical hypotheses on how the hydrologic system functions, characterizing the potential changes of the water resources using precipitation and temperature data as inputs variables. Calibration and validation of hydrologicmodels are required in order to develop reliable models (Savenije, 2009), and sensitivity analysis should be carried out for a better knowledge of complex models (Song et al., 2015). Moreover, the recognition of the equifinality concept, that is, the existence of many sets of parameters conducive to good adjustments to some target observations (Beven, 2006; 2012), is unavoidable and necessary (Beven and Freer, 2001). In the context of climate modeling, these models are usually called Land-Surface Models (LSMs) and are coupled to General Circulation Models (GCMs) and Regional Climate Models (RCMs) as the land scheme that allows to simulate the biophysical processes involved in the land-atmosphere interaction (Wood et al., 2011). Although there is a subtle difference between a hydrologic model and a LSM, this distinction has become blurred over time (Clark et al., 2015). In this respect, the Variable Infiltration Capacity (VIC) macroscale hydrologic model (Liang et al., 1994, 1996) has played the role of both LSM and hydrologic model in many previous studies. Melsen et al. (2016b) provided sufficient evidence in a meta-analysis of 192 peer-reviewed studies where the VIC model was calibrated and validated. Since its first development many efforts have been made in order to study the sensitivity of the VIC model, which has been explored in a broad sense in the following terms: -
Sensitivity to spatio-temporal variability : the impacts of the implemented spatial resolution in the simulated runoff and other water fluxes have been addressed in various studies, suggesting that there is a high influence of the sub-grid variability of the precipitation on the model performance (Haddeland et al., 2002; Liang et al., 2004). However, a critical spatial resolution under which a better model performance is not necessarily achieved (Liang et al., 2004) could exist. These impacts are also noticeable in calibration and validation exercises with an increase of the model accuracy at higher resolutions (Oubeidillah et al., 2014), although a high transferability of the calibrated parameters across the different resolutions may be an indicator of a poor representation of the spatial variability (Melsen et al., 2016a). Unfortunately, the time step of the calibration and validation has not kept up with the increasing spatial resolution, and this is a crucial aspect for the correct representation of the involved hydrologic processes (Melsen et al., 2016b). The fact that it is more difficult to transfer parameters across temporal resolutions than across the spatial dimension brings the need of a better representation of the spatial variability in macroscale hydrologic models (Melsen et al., 2016a). -
Sensitivity to boundary conditions : understanding the boundary conditions as the meteorological forcings that drive the simulations, the VIC model sensitivity to the boundary conditions has been studied through the application of different climate change scenarios and the analysis of the impacts of changing precipitation and temperature on the hydrology and water resources of several continental river basins (Nijssen et al., 2001), the Pacific Northwest (Vano et al., 2015), or more locally in the Colorado River Basin (Christensen et al., 2004; Bennet et al., 2018) and in the upper Ganga Basin (Chawla and Mujumdar, 2015). Also, these kind of studies sometimes are carried out in conjunction with other sensitivity analysis, i.e. the combined and segregated effects of climate change and land use changes on streamflow (Chawla and Mujumdar, 2015) or the parameter sensitivity under a changing climate (Bennet et al., 2018). -
Sensitivity to initial conditions : the question about if the hydrologic predictions are affected by the hydrologic initial conditions (i.e. the initial moisture state at the snowpack and the soil profile) or if the boundary conditions constitute the main contributor to the model simulations was studied in detail in Cosgrove et al. (2003), Wood and Lettenmaier (2009) and Li et al. (2009). It is known that the soil moisture content for the bottom soil layer of the VIC grid cells is the variable that commonly takes the longest time to reach the equilibrium, and although there is not a general agreement in how long the model spin-up period should be, since it highly depends on each particular application, it has been found that wetter states lead to faster spin-up times (Cosgrove et al., 2003; Melsen et al., 2016a). This issue is of relevance for a proper calibration and validation of the VIC model and is usually avoided by fixing a long-enough spin-up period previous to the simulation together with a wet initialization of the soil layers. -
Sensitivity to soil and vegetation parameters : since its generalization as a three-layer soil model (VIC-3L) in Liang et al. (1996) after its two-layer predecessor (VIC-2L, Liang et al., 1994), the sensitivity of the model to soil parameters has been studied at basin-scale and at global-scale through a cell-based approach. The basin-scale approach of Demaria et al. (2007) allowed to estimate the sensitivity of the simulated streamflows to the parameters that control the surface and subsurface runoff generation, and the global-scale study of Chaney et al. (2015) evaluated the efficiency of the VIC model for monitoring global floods and droughts under a parameter uncertainty framework. The sensitivity to land use changes and the vegetation parameters associated to the different vegetation classes (i.e. Leaf Area Index and albedo) have also been explored and expressed for the different components of the water and energy balances from the VIC model (VanShaar et al., 2002; Chawla and Mujumdar, 2015; Bennett et al., 2018). This work aims to contribute to the knowledge of the VIC model in a macroscale application over the headwater subwatersheds of an important basin located in the north of the Iberian Peninsula, the Duero River Basin. For this end, the hydrologic modeling exercise here developed has been divided into three interrelated parts: - The calibration and validation of the VIC model for the selected subwatersheds of the study area. - An integrated sensitivity analysis for all the subwatersheds focused on the soil parameters chosen for the calibration. - A final assessment of equifinality and the efficiency of the calibration algorithm that links the calibration and sensitivity analysis results. The Duero River Basin has been investigated in various previous studies and the main issues addressed are: the temporal trend of water supply and its relation to precipitation, temperature and plant cover changes (Ceballos-Barbancho et al., 2008); the hydrologic response to land-cover changes (Morán-Tejeda et al., 2010, 2012a), the impacts of different climate oscillations (Morán-Tejeda et al., 2010) and its response to the North Atlantic Oscillation (Morán-Tejeda et al., 2011a); the characteristics of the different existing river regimes (Morán-Tejeda et al., 2011b) and the effects of reservoirs on them (Morán-Tejeda et al., 2012b). However, all these studies were based on statistical analyses of different hydroclimatic and land-surface variables and contributed to a better understanding of the hydrologic behavior of the Duero RiverBasin. Therefore, the hydrologic modeling analysis carried out in this work can provide an added value to this set of issues since the potentialities of a macroscale hydrologic model such as the VIC model have been examined in detail for this river. In Sect. 2 and 3 the study area and the methods are described. Sect. 4 gathers the results of the three-part approach and Sect. 5 corresponds to the discussion of the key results. Finally, the main conclusions of this study are provided in Sect. 6.
2. Study area
The Duero River Basin constitutes the largest basin of the Iberian Peninsula with a surface of 98 073 km . It is a shared territory between Spain and Portugal, characterized by a high water contribution (~ 15 000 hm /year). The study is focused on the Spanish part of the basin (Fig. 1), which represents the 80% of the area (78 859 km ). Most of this territory constitutes a plain surrounded by mountainous chains, thus configuring two topographic areas well differentiated. The large depression is filled with sediments of the Tertiary and the Quaternary, constituting a complex hydrogeologic environment. The lithology of the northern mountains consists of siliceous, calcareous and carbonated rocks with local small aquifers to the west part and aquifers of greater capacity to the east. The south system harbors rocks of low permeability and is dominated by a granite batholith. Finally, the eastern mountainous areas hold a silicic core enclosed by carbonated rocks with a high presence of karstic aquifers. The basin presents a predominant Mediterranean climate with a mean annual precipitation volume of 50 000 hm which is mostly lost into the atmosphere through evaporative fluxes (~ 35 000 hm /year). Most precipitation is concentrated in the mountainous areas reaching values above 1500 mm/year to the north of the basin and values slightly below 1000 mm/year to the south and east. As for the most part of the Iberian Peninsula, precipitation exhibits a very irregular intra-annual distribution, being concentrated in spring and fall and almost nonexistent during summer. Winter months are cold with a mean temperature of 2ºC in January, while summer is soft with maximum temperatures occurring in July (~ 20.5ºC). The Duero River Basin is regulated by a total of 31 reservoirs where the streamflow records are estimated through a water balance of the daily storages and water releases. The streamflow monitoring is also carried out in a large network of ca. 200 gauging stations where the streamflow records are calculated through the rating curves.
3. Methods
The streamflow records were gathered in a monthly basis from the Spanish Centre for Public Work Experimentation and Study (CEDEX,
Centro de Estudios y Experimentación de Obras Públicas ) database for all the reservoirs and gauging stations of the Duero Basin. An analysis of the percentage of gaps in the time series revealed that the period from October 2000 to September 2011 presents less than 5% of missing values for all the time series. Therefore, it was chosen as the study period for this work. A reference hydrologic network was then defined applying the criterion of the absence of upstream hydraulic structures (Whitfield et al., 2012). This selection reduced the number of reservoirs and gauging stations for the analysis to 16 reservoirs and 15 gauging stations covering the headwaters of the Duero River Basin (Fig. 1). Table 1 collects the main characteristics of these subwatersheds: area (km ), mean elevation (m), averaged annual precipitation ( P an , mm/year), potential evapotranspiration ( PET an , mm/year), and streamflow ( Q an , hm /year), for the study period. Precipitation, maximum temperature and minimum temperature data were extracted from two high-resolution (~ 5 km x 5 km) daily gridded datasets: SPREAD (Serrano-Notivoli et al., 2017) for precipitation data and STEAD (available at http://dx.doi.org/10.20350/digitalCSIC/8622) for temperature data. Both datasets cover Peninsular Spain and the Balearic and Canary Islands and were built with information from observed precipitation and maximum and minimum temperature at a varying number of meteorological stations provided by several administrations including the Spanish Meteorological Agency (AEMET) and some hydrologic confederations. Atmospheric pressure, incoming shortwave radiation, incoming longwave radiation, vapor pressure and wind speed data were taken from daily outputs of high resolution (0.088º,
10 km) simulations carried out with the Weather Research and Forecasting 0 (WRF) model driven by the ERA-Interim Reanalysis data for the spatial domain of the Iberian Peninsula (García-Valdecasas Ojeda et al., 2017).
For each time step Q d [L] is the surface (direct) runoff, P [L] is the precipitation, z [L] is the depth of the upper two soil layers, θ is their volumetric soil moisture content, θ S is their porosity, i m [L] is the maximum infiltration capacity, i [L] is the infiltration capacity that corresponds to the soil moisture at that time step and b i is the infiltration shape parameter. Baseflow is generated in the third soil layer following the Arno formulation (Franchini and Pacciani, 1991), and is expressed as: Here, Q b [L] is the baseflow for each time step, D m [L] is the maximum baseflow, D S is a fraction of D m , θ is the volumetric soil moisture content of the soil layer 3, θ S is the porosity in this layer and W S is a fraction of θ S . The baseflow recession curve is divided into two parts: a linear part for lower values of θ and a non-linear (quadratic) part for higher values of θ . The characteristic large grid size for macroscale models makes that the VIC snow model conceptualizes the snow processes partitioning each grid cell of the spatial domain into snow bands, thereby accounting for the sub-grid variability in topography, land uses and precipitation. The snow model is applied separately to each snow band and land class and the outputs consist of the snow depth and the snow water equivalent for each grid cell. The snowpack is represented as a two-layer model that solves the energy and the mass balance and determines whether the snowpack is subject to 2 accumulation or ablation, making the model suitable for applications in any part in the world (Liang et al., 1994, 1996). 3.2.2 Soil and vegetation parameters The required soil parameters for the application of the VIC model were obtained from SoilGrids1km (Hengl et al., 2014) and EU-SoilHydroGrids ver1.0 (Tóth et al., 2017), all of them with a spatial resolution of 1 km x 1 km. In both datasets the different soil properties are provided for seven soil depths up to 2 m (0, 5, 15, 30, 60, 100 and 200 cm). These soil parameters are: (1) bulk density and soil textural classes of the United States Department of Agriculture (USDA) from SoilGrids1km; and (2), field capacity, saturated hydraulic conductivity, porosity and wilting point from EU-SoilHydroGrids ver1.0. The VIC model handles land uses information as a set of vegetation parameters for the different vegetation classes specified in a vegetation library. Here the UMD Global land cover classification (Hansen et al., 2000) was chosen with a spatial resolution of 1 km x 1 km, and the vegetation parameters (i.e. Leaf-Area Index, roots depth and roots coverage) were fixed for each vegetation class following the recommendations of the GLDAS Project (https://ldas.gsfc.nasa.gov/gldas/vegetation-parameters). 3.2.3 Routing procedure The water balance mode of the VIC model allows to simulate the surface runoff and the baseflow for each grid cell of the spatial domain where the model is implemented. Since the conceptualization of the VIC model does not include the horizontal transport processes between contiguous grid cells, the runoff generated needs to be routed to a 3 certain outlet in order to determine simulated values of streamflow that can be compared with existing observations. An approach consisting of a monthly aggregation of the total runoff (the sum of surface runoff and baseflow) for all the grid cells within a given subwatershed, weighted by the fractional area of each cell inside the subwatershed, was selected for this work. This approach has the advantage of working with the real boundaries of the subwatersheds and is not dependent on the resolution of the grid cells. 3.2.4 Model implementation The VIC model was implemented in the water balance mode at a daily time step and with a spatial resolution of 0.05º (~ 5 km x 5 km) for the 31 studied subwatersheds (Fig. 2). The meteorological forcings were interpolated to the grid cells of the spatial domain following a nearest neighbor assignment. The soil parameters and the elevation wereaveraged for each grid cell, and the vegetation parameters were kept at the original resolution of 1 km x 1 km because the VIC model allows to consider the sub-grid variability of the land uses for each grid cell. As in the soil database, the depth of each grid cell was set at 2 m, fixing the thickness of the first soil layer ( d ) at 0.1 m and varying the thicknesses of the second ( d ) and third ( d ) layers during the calibration and the sensitivity analysis. The outputs of the VIC model were finally aggregated into monthly values for each subwatershed in order to accomplish the calibration, validation and sensitivity analysis.
4 The model was calibrated for the period from October 2000 to September 2009 choosing the Nash-Sutcliffe Efficiency (
NSE , Nash and Sutcliffe, 1970) as the objective function. The
NSE was calculated by comparing the monthly observations of streamflow with the monthly aggregated total runoff simulated by the VIC model. Table 2 indicates the selected parameters for the calibration process together with their upper and lower bounds. The calibration was carried out using the Shuffled-Complex-Evolution Algorithm (SCE-UA) of Duan et al. (1994). A spin-up period of ten years, previous to the calibration period, was simulated in order to ensure that the soil moisture content of the three soil layers reached an equilibrium, and therefore the initial conditions did not affect the calibration process. The model was then validated for the period October 2009 – September 2011. The goodness-of-fit for both calibration and validation exercises was evaluated through the
NSE and the coefficient of determination r . The Standardized Regression Coefficients (SRC) method (Saltelli et al., 2008) aims to study the propagation of uncertainty from model inputs to outputs. The SRC method is focused on the behavior of the model outputs in relation to a certain set of parameters once the boundary conditions (i.e. the meteorological forcings) and the initial conditions (i.e. the soil moisture content of the three soil layers) have been fixed. This sensitivity analysis method requires two elements: first, a Monte Carlo simulation where the model is run with a specified number of parameter samples; and second, a multiple linear regression of each model output of interest as a linear function of the parameters. 5 The sensitivity analysis was carried out for each subwatershed considered in the study area. As in the calibration process, the period October 2000 – September 2009 was chosen and ten years of spin-up prior to the study period were run. 3.4.1 Monte Carlo simulation A parametric space is defined through the selection of several parameters and their upper and lower bounds. Here, a 5-dimensional parametric space was established choosing the five calibration parameters and considering the upper a lower bounds specified in Table 2. A sampling method is then applied to the parametric space, extracting a large-enough sample of parameters for the Monte Carlo simulation. The Latin Hypercube Sampling (LHS) method (Iman and Conover, 1982) was applied for this step extracting a total of 10 000 random samples. This process allowed to define a sampling matrix, Θ , of order m x n , where m represents the number of samples ( m = 10 000) and n the number of parameters for the analysis ( n = 5). The model was finally run for each parameter combination (i.e. row) of Θ . The Monte Carlo simulation was also used for assessing equifinality and the efficiency of the calibration algorithm by studying the response given by each parameter sample in terms of the
NSE.
The results were compared with the
NSE determined during the calibration period. 3.4.2 Multiple linear regression The outputs of interest from the VIC model were those components included in the water balance: surface runoff ( Q d ), baseflow ( Q b ), total runoff ( Q t ), actual evapotranspiration ( AET ) and the soil moisture content of the three soil layers ( SM , SM and SM ). For each component and for each run of the Monte Carlo simulation, the 6 mean value of the simulated series was calculated, and a multiple linear regression model was then adjusted relating the mean values of each component with the sampling parameters: Where y is a column vector with the m mean values of the component, a is the intercept of the hyperplane, a i is the regression coefficient of the parameter i and Θ i is the column of the sampling matrix corresponding to the parameter i . The standardized regression coefficients, β i , are then calculated for each parameter: Here, σ Θ i and σ y p are the standard deviations of Θ i and the predicted values of y , respectively. β i2 represents the relative contribution of the parameter i to the variance of the model output of interest, being and equal to the coefficient of determination r of the adjustment. A threshold of r ≥ β i are valid measures of the sensitivity (Saltelli et al., 2006), although they can be robust and reliable measures even for nonlinear models (Saltelli et al., 2008). β i can take values between -1 and 1. A high absolute value of β i implies that the component is sensitive to the parameter and its sign indicates whether the effect is positive or negative.
4. Results
The values of
NSE and r for the calibration and validation periods are shown in Table 3. Figure 3 depicts the simulated streamflows during both periods together with the 7 observed streamflows for six selected subwatersheds (R-2011, R-2037, R-2038, GS-3005, GS-3089 and GS-3150) located in different parts of the basin. The NSE for the calibration period presents values above 0.75 in 19 out of the 31 subwatersheds and reaches values above 0.85 in 10 subwatersheds, and the corresponding r values are high too. For the validation period, both NSE and r values are predominantly high, and generally lower than the corresponding ones for the calibration process, although minimum values of NSE below 0 were attained for 3 subwatersheds. Note that the results of the calibration and the validation processes are slightly better for the reservoirs, what could indicate a quality difference between the streamflow databases from the reservoirs and from gauging stations. Some of the high r values are obtained for low NSE estimations, which indicates that the model is able to capture the intra-annual variability of the streamflow observations but is not able to reach a good fit for the peaks of streamflow (Table 3, Fig. 3). It is interesting to note that high
NSE values were obtained for subwatersheds with varying sizes, with good fits for both small-sized (e.g. R-2011 and GS-3089) and medium-sized (e.g. R-2038 and GS-3005) subwatersheds, emphasizing the ability of the VIC model to provide accurate predictions of the streamflow across different spatial scales.
Through the application of the SRC method the β coefficients for the five calibration parameters of the water balance components in the VIC model were obtained, and the results are shown in Fig. 4 (a to g) for all the subwatersheds. The r value obtained from the multiple linear regression and the estimation of r as the sum of the squares of β coefficients are also depicted in Fig. 4 (h, i), reflecting very similar values. The results of the sensitivity analysis for the selected components to the parameters are given below 8 in a component-by-component basis providing the necessary explanations when there is a strong dependency between them: - Q d : the values of r are above 0.7 for all the subwatersheds (Fig. 4), fulfilling the criterion of enough linearity for interpreting the results of the sensitivity analysis. The strongest positive effect corresponds to the parameter b i , which means that a higher value of b i leads to more surface runoff. This is clearly evidenced in Eq. 1, where a relation of exponential type between Q d and b i is established. D m produces a negative effect on the surface runoff, suggesting that an increase of the maximum baseflow brings a reduction of the surface component under the assumption of the same meteorological forcings. d also yields a negative effect on Q d , and this effect is related to an increase in AET . - Q b : the values of r mostly range from 0.5 to 0.7, with some values close to 0.8 (Fig. 4). In this case the linearity criterion is hardly reached and therefore it is difficult to interpret the β coefficients. However, it is interesting to note that, with the exception of d , the β coefficients of the parameters are characterized by a low dispersion. This is an indicator of the robustness of the VIC model response, and although the threshold of linearity is not always achieved, the dependency of Q b with respect to these parameters can be accepted. As expected from the previous analysis of the surface component of the runoff, b i reflects a strong negative effect. The positive effects now correspond to D S and D m . This is obvious in the case of D m but not so evident for D S since a higher value of D S only means that the baseflow law tends to be more linear (see Eq. 2). The amplitude of the β coefficients for d is broader than for the rest of the parameters but always negative except for two subwatersheds. - Q t : the total runoff exhibits an additive effect of the previous components for both r and the β coefficients as it is computed through the sum of the surface runoff 9 and the baseflow (Fig. 4). Thus, higher values of b i , D S and D m lead to an increase of Q t and higher values of d produce a negative effect on Q t due to a rise of AET . This component is of particular interest given that it is the component subject to calibration in this work. In order to provide a better understanding on its behavior the spaghetti plots of the Monte Carlo simulation for an example subwatershed (R-2038) are depicted in Fig. 5. As shown there, the time series of the observed streamflow (Fig. 5c) falls into the range of responses of the model for almost all the study period and therefore one or more sets of parameters will afford a god fit with the observations. -
AET : this component and Q t are linked through the law of conservation of mass applied to the system defined by each subwatershed, being the precipitation equal to the sum of Q t , AET and the variation of the storage in the hydrologic system. Moreover, the study period for the sensitivity analysis is long-enough to neglect the last term of the water balance equation, and the precipitation is fixed for each subwatershed as a boundary condition. In consequence, the linearity of both Q t and AET with respect to the parameters must be similar and the β coefficients for AET are essentially identical to the corresponding ones for Q t but with opposite signs (Fig. 4). Figure 5d shows the spaghetti plots of this component together with the potential evapotranspiration ( PET ) profile. The reason of the existence of some values of
AET above the
PET curve responds to the internal handling of the Penman-Monteith equation used in the VIC model because various different approaches are considered when computing the potential evapotranspiration, and the curve presented in Fig. 5d corresponds to the current vegetation parameters. - SM : the values of r are widely scattered and range from 0.35 to values above 0.9 (Fig. 4). The nature of such a scattered distribution may be an outcome of the closeness between the soil moisture profiles in this layer, making it difficult to adjust a 0 multiple linear regression model to its mean values. Similarly to the case of Q b , most of the β coefficients present a relatively low dispersion and subsequently the results of the sensitivity analysis can be interpreted. The negative effects mainly concern to b i and D m , demonstrating that a higher exponent in the surface runoff equation and a higher maximum baseflow are related to lower soil moisture values for the upper soil layer. On the other hand, the positive effects are associated with increasing values of d despite of revealing highly dispersed β coefficients. The spaghetti plots of SM (Fig. 5e) seem to reproduce the PET cycles and this results from the evaporative fluxes themselves as the transpiration process occurs from the roots of the vegetation. - SM : it is the component with the highest linearity with regard to the calibration parameters, proffering values of r very close to 1 for all the subwatersheds (Fig. 4). In this case d dominates the sensitivity of SM with a noticeable positive effect (i.e. values of β near to 1). The PET cycles are also markedly represented in the soil moisture profiles of this layer (Fig. 5f). - SM : even though the values of r lie between 0.6 and 0.7 predominantly, some of them fall below 0.6 with minimum values close to 0.4 (Fig. 4). Once again the dispersion of the β coefficients is relatively low and in this occasion this is also true for d . Baseflow takes place from this layer and this is reflected in the β coefficients corresponding to D S , W S and D m , which present opposite signs and similar absolute values to the calculated ones for Q b . As for the previous soil layers, the PET cycles are present too but here there is a lag in the valleys of the soil moisture profiles due to the delay in the baseflow generation process (Fig. 5g).
1 Equifinality and the efficiency of the calibration algorithm were assessed through the evaluation of the
NSE values for the Monte Carlo simulations of all the subwatersheds by comparing the total runoff of each simulation with the observed streamflow during the calibration period. For this purpose, two counts of the number of simulations satisfying certain criteria were carried out: first, the number of simulations for each subwatershed presenting
NSE values above the
NSE determined during the calibration (
NSE cal ) minus 0.05 was used as indicator of equifinality of the VIC model and the parameter samples; and second, the number of simulations with
NSE values above
NSE cal hinted at the efficiency of the SCE-UA algorithm in finding the optimal set of parameters producing the best fit with the streamflow records. The results of this exercise are expressed in Table 4. It is clear that for the majority of the subwatersheds there are many simulations with
NSE values very close to the optimal model, and in some cases the number of simulations is fairly high (> 3000). This can be also appreciated when the columns of the sampling matrix are plotted against the
NSE of each simulation in a “dotty plot”.
Figure 6 shows the dotty plots of the five parameters of the calibration for two subwatersheds (R-2038 and GS-3089) as an example of this analysis. For both subwatersheds the
NSE values of the simulations are above 0 and the points clouds are concentrated on the top of the diagrams, suggesting that a high number of them are close to the optimal fit (see also Table 4). The shape of the dotty plots supplies useful information about the behavior of the parameter samples as a set. For example, the dotty plot of d for the subwatershed R-2038 reflects a trend to produce high NSE values when the parameter values are near to the upper bound. The optimum was reached for d = 0.8995 m, while the rest of the fitted parameters were located between the fixed limits. Also, most of the NSE values were below
NSE cal when the second count was 2 executed, implying that the SCE-UA algorithm was highly efficient in searching for the optimal set of parameters.
5. Discussion
The results of the calibration and validation processes are comparable to other studies using hydrologic models developed in northern Spain (Morán-Tejeda et al., 2014) and recently in the south of Spain (Pellicer-Martínez and Martínez-Paz, 2018;
Yeste et al., 2018). It is to be expected that the application of a single model structure over an heterogeneous spatial domain, such as the Duero River Basin, does not conduct to a good adjustment of the simulated streamflow with the observations for all the studied subwatersheds. Furthermore, the existence of other potential pressures over the water resources may be responsible for those cases where the calibration and/or validation showed poor results, and therefore further research is required in order to identify the origin of the biases with respect to the observations of the simulated streamflows for these subwatersheds. Nevertheless, the results of the calibration and the validation suggest that the macroscale application of the VIC model carried out in this study performs well for a large number of subwatersheds in the Duero River Basin. Moreover, the routing procedure has proven to be accurate and efficient in this work, permitting its application in other studies using hydrologic models that operate over the grid cell in a similar mode to the VIC model. Concerning the sensitivity analysis, the application of the SRC method allowed a deep understanding of the existing relationships between the components of the water balance in the VIC model and the selected parameters for the calibration as long as the linearity criterion was fulfilled. Even when the coefficient of determination of the fitted model did not satisfy the linearity criterion, the relatively low dispersion of the β
3 coefficients permitted the interpretation of the results. Special attention deserves the component Q t since it is the component that was compared with the streamflow observations during the calibration and validation processes. The sensitivity of this component to the five soil parameters reflected an additive effect of the sensitivity measures of Q d and Q b as Q t is calculated as the sum of the surface and the subsurface components. Q t was mainly sensitive to b i and d , and this is consistent with the sensitivity measures for the simulated streamflow carried out in Demaria et al. (2007) for four studied subwatersheds. At the sight of the results of the equifinality assessment, it is unavoidable accepting that no parameter set leads to a single optimal model, or in other words, that there are many parameter samples with performances as good as the optimum calculated with the calibration algorithm. As in the GLUE method (Beven and Binley, 1991; Beven, 2012), this fact could be the starting point of the calibration process, in which a measure of belief is associated to each parameter set according to the degree of proximity to the optimum. This will be an interesting research line for further investigation in the Duero River Basin. In any case, we consider that the use of a calibration algorithm provides a first-look into the goodness-of-fit response surface of the hydrologic model in a computationally more efficient way than the Monte Carlo experiment, serving as a sign of the goodness-of-fit of the overall parameter samples.
6. Conclusions
The main conclusions of this work can be summarized as follows: [1] The calibration and validation of the VIC model reflected good results for most of the studied subwatersheds in the Duero River Basin with a predominance of high
NSE values. The results were slightly better for the reservoirs than for the gauging 4 stations and this may be a consequence of a quality difference between the streamflow databases. The calibration and validation processes showed poor results in a few subwatersheds. This may be caused by the existence of pressures over the water resources that have not been taken into account in the modeling exercise. However, this is out of the scope of this work since the main interest is placed on the macroscale application of the VIC model, which has shown to perform well for a great part of the Duero River Basin. [2] The β coefficients calculated during the sensitivity analysis allowed to quantify the sensitivity of the water balance components to the selected parameters for the calibration. The surface runoff and the soil moisture content of the soil layer 2 were the components with the highest linearity and were mainly dominated by the values of the infiltration shape parameter and the thickness of soil layer 2, respectively, both with a positive effect. The total runoff presented a combined behavior from the surface runoff and the baseflow components, and the sensitivity analysis yielded similar results to other sensitivity measures previously reported in the literature. The potential evapotranspiration cycles were noticeable in the whole soil profile and more evidently in the upper two soil layers. [3] A final exercise for assessing equifinality and the efficiency of the calibration algorithm was carried out, finding that there are many parameter sets with NSE values as high as the
NSE determined during the calibration. The calibration algorithm was efficient and reached the optimal fit for almost all the studied subwatersheds. The use of a calibration algorithm is also in line with other possible practical applications of the VIC model for studying the impacts of climate change on water resources in the Duero River Basin, where a parameter set must be chosen prior to the simulations using climate change data. 5
Acknowledgements
All the simulations were conducted in the ALHAMBRA cluster (http://alhambra.ugr.es/) of the University of Granada. This work was partially funded by the Spanish Ministry of Economy and Competitiveness projects CGL2013-48539-R and CGL2017-89836-R, with additional support from the European Community Funds (FEDER). The first author was supported by the Ministry of Education, Culture and Sport of Spain (FPU grant FPU17/02098).
References
Bennett, K.E., Urrego Blanco, J.R., Jonko, A., Bohn, T.J., Atchley, A.L., Urban, N.M., Middleton, R.S., 2018. Global Sensitivity of Simulated Water Balance Indicators Under Future Climate Change in the Colorado Basin. Water Resour. Res. 54, 132 – –
36. https://doi.org/10.1016/j.jhydrol.2005.07.007 Beven, K., 2012. Rainfall-Runoff Modelling: The Primer: Second Edition. Rainfall-Runoff Model. Prim. Second Ed. https://doi.org/10.1002/9781119951001 Beven, K., Binley, A., 1992. The future of distributed models: Model calibration and uncertainty prediction. Hydrol. Process. 6, 279 – –
29. https://doi.org/10.1016/S0022-1694(01)00421-8 6 Ceballos-Barbancho, A., Morán-Tejeda, E., Luengo-Ugidos, M.Á., Llorente-Pinto, J.M., 2008. Water resources and environmental change in a Mediterranean environment: The south-west sector of the Duero river basin (Spain). J. Hydrol. 351, 126 – – – ‐ up behavior in the North American Land Data Assimilation System (NLDAS). J. Geophys. Res. Atmos. https://doi.org/10.1029/2002jd003316 Demaria, E.M., Nijssen, B., Wagener, T., 2007. Monte Carlo sensitivity analysis of land surface parameters using the Variable Infiltration Capacity model. J. Geophys. Res. Atmos. https://doi.org/10.1029/2006JD007534 Diffenbaugh, N.S., Giorgi, F., 2012. Climate change hotspots in the CMIP5 global climate model ensemble. Clim. Change 114, 813 – – – – Impacts, Adaptation and Vulnerability: Part B: Regional Aspects: Working Group II Contribution to the IPCC Fifth Assessment Report, vol. 2, Cambridge University Press, Cambridge, UK, https://doi.org/10.1017/CBO9781107415386 Li, H., Luo, L., Wood, E.F., Schaake, J., 2009. The role of initial conditions and forcing uncertainties in seasonal hydrologic forecasting. J. Geophys. Res. Atmos. https://doi.org/10.1029/2008JD010969 Liang, X., Guo, J., Leung, L.R., 2004. Assessment of the effects of spatial resolutionson daily water flux simulations. J. Hydrol. 9 https://doi.org/10.1016/j.jhydrol.2003.07.007 Liang, X., Lettenmaier, D.P., Wood, E.F., Burges, S.J., 1994. A simple hydrologically based model of land surface water and energy fluxes for general circulation models. J. Geophys. Res. https://doi.org/10.1029/94JD00483 Liang, X., Wood, E.F., Lettenmaier, D.P., 1996. Surface soil moisture parameterization of the VIC-2L model: Evaluation and modification. Glob. Planet. Change. https://doi.org/10.1016/0921-8181(95)00046-1 Lorenzo-
Lacruz, J., Mo ŕ an -Tejeda, E., Vicente- Serrano, S.M., Ĺ opez -Moreno, J.I., 2013. Streamflow droughts in the Iberian Peninsula between 1945 and 2005: Spatial and temporal patterns. Hydrol. Earth Syst. Sci. 17, 119 – – – –
49. https://doi.org/10.1016/j.gloplacha.2010.03.003 Morán-Tejeda, E., Ceballos-Barbancho, A., Llorente-Pinto, J.M., López-Moreno, J.I., 2012a. Land-cover changes and recent hydrological evolution in the Duero Basin (Spain). Reg. Environ. Chang. 12, 17 –
33. https://doi.org/10.1007/s10113-011-0236-7 Morán-Tejeda, E., Ignacio, L.M., Antonio, C.B., Sergio M., V.S., 2011a. Evaluating
Duero’s basin (Spain) response to the NAO phases: Spatial and seasonal variability.
Hydrol. Process. 25, 1313 – – – – Nijssen, B., O’donnell, G.M ., Hamlet, A.F., Lettenmaier, D.P., 2001. Hydrologic sensitivity of global rivers to climate change. Clim. Change 50, 143 – Savenije, H.H.G., 2009. HESS opinions: “The art of hydrology.” Hydrol. Earth Syst.
Sci. 13, 157 – – – surface modeling: Meeting a grand challenge for monitoring Earth’s terrestrial water. Water Resour. Res. 47. https://doi.org/10.1029/2010WR010090 Yeste, P., Dorador, J., Martin-Rosales, W., Molero, E., Esteban-Parra, M.J., Rueda, F.J., 2018. Climate-driven trends in the streamflow records of a reference hydrologic network in Southern Spain. J. Hydrol. 566, 55 –
72. https://doi.org/10.1016/j.jhydrol.2018.08.063 Zhao, R.J., Zhuang, Y.L., Fang, L.R., Liu, X.R., Zhang, Q.S., 1980. The Xinanjiang Model, Hydrological Forecasting. Oxford Symp. IAHS129. Figure captions Figure 1.
Duero River Basin and the 31 studied subwatersheds. The prefix “R - ” denotes “Reservoir” and the prefix “GS - ” denotes “Gauging Station”. Figure 2.
VIC model implementation.
Figure 3.
Time series of the observed streamflows along with the simulated ones for the calibration and de validation periods for six example subwatersheds.
Figure 4. (a to g) β coefficients for the 31 subwatersheds. h) r value from the multiple linear regression. i) r estimated from the β coefficients. Figure 5.
Spaghetti plots of the water balance components resulting from the Monte Carlo simulation for the subwatershed R-2038.
Figure 6.
Dotty plots for two subwatersheds: a) R-2038 and b) GS-3089. Red dot corresponds to the calibrated value for the corresponding parameter using the SCE-UA algorithm. -2038 R-2037R-2036R-2032 R-2043R-2042R-2039 R-2030R-2027R-2026 R-2013 R-2012R-2011R-2028 R-2001GS-3818 GS-3150GS-3124GS-3105GS-3104 GS-3089 GS-3057GS-3051 GS-3049GS-3047 GS-3041GS-3028GS-3016 GS-3005R-2014 GS-3035 ° ' N ° N ° ' N ° N ° ' N ° N ° ' N ° N DUERO RIVER BASIN
Geodetic System: ETRS89Units: Degrees
SITUATION
Spain
Gauging PointDrainage networkSubwatershed< 500500 - 750750 - 10001000 - 12501250 - 15001500 - 17501750 - 2000> 2000Elevation (m) ubwatershed boundaryGrid cellDuero River Basin
Model resolution: 0.05º Q ( h m / m o n t h ) GS-3089 Q ( h m / m o n t h ) R-2037 Q ( h m / m o n t h ) R-2011 Q ( h m / m o n t h ) R-2038
ValidationCalibration Q ( h m / m o n t h ) GS-3150 Q ( h m / m o n t h ) GS-3005 Q CAL Q VAL Q OBS i D S W S D m d β -1.0-0.50.00.51.0 b i D S W S D m d β -1.0-0.50.00.51.0 b i D S W S D m d β -1.0-0.50.00.51.0 b i D S W S D m d β -1.0-0.50.00.51.0 b i D S W S D m d β -1.0-0.50.00.51.0 b i D S W S D m d β -1.0-0.50.00.51.0 b i D S W S D m d β Q d Q b Q t A E T S M S M S M r Q d Q b Q t A E T S M S M S M r a) Q d b) Q b c) Q t d) AET e) SM f) SM g) SM h) r i) r bserved streamflowPET Q d ( h m / m o n t h ) Q b ( h m / m o n t h ) Q t ( h m / m o n t h ) A E T ( mm / m o n t h ) S M ( mm ) S M ( mm ) S M ( mm ) PET cycles alibrated valuesNSE = 0.9145 b i = 0.3240 D S = 0.6915W S = 0.5536D m (mm/d) = 24.1871 d (m) = 0.8995 a) R-2038b) GS-3089 Calibrated valuesNSE = 0.9116 b i = 0.2487D S = 0.5722W S = 0.6330D m (mm/d) = 17.1926 d (m) = 0.5914 able 1. Main characteristics of the 31 subwatersheds.
Code Name Area (km ) Mean elevation (m) P an (mm/year) PET an (mm/year) Q an (hm /year) R-2001 CUERDA DEL POZO 546.7 1319 1122 900 188.7 R-2011 ARLANZON 106.7 1440 1293 798 77.6 R-2012 CERVERA DE RUESGA 54.3 1284 1005 859 86.1 R-2013 REQUEJADA, LA 220.7 1378 1022 780 153.2 R-2014 CAMPORREDONDO 229.6 1673 1457 747 223.9 R-2026 BARRIOS DE LUNA 482.9 1496 1468 757 400.1 R-2027 VILLAMECA 45.8 1180 946 941 34.1 R-2028 MONCABRIL (SISTEMA) 62.9 1712 1242 744 93.5 R-2030 PORMA / JUAN BENET 250.3 1412 1251 753 304.5 R-2032 RIAÑO 592.3 1451 1569 748 623.7 R-2036 LINARES DEL ARROYO 761.3 1111 515 1222 52.4 R-2037 BURGOMILLODO 803.1 1097 573 1232 87.0 R-2038 SANTA TERESA 1845.4 1326 882 1235 734.8 R-2039 AGUEDA 788.4 895 1189 1333 392.7 R-2042 CASTRO DE LAS COGOTAS 853.3 1279 527 1231 89.0 R-2043 PONTON ALTO 150.4 1582 990 993 78.3 GS-3005 OSMA 893.1 1090 728 1084 121.0 GS-3016 PAJARES DE PEDRAZA 281.3 1298 634 1196 61.6 GS-3028 SALAS DE LOS INFANTES 353.5 1257 988 926 104.7 GS-3035 OTERO DE GUARDO 69.2 1492 1702 838 29.3 GS-3041 VILLALCAZAR DE SIRGA 307.7 929 664 1011 35.5 GS-3047 MEDIANA DE VOLTOYA 130.4 1347 518 1128 15.4 GS-3049 CABAÑES DE ESGUEVA 270.2 995 658 1046 24.7 GS-3051 ESPINAR, EL 36.7 1610 905 1032 18.8 GS-3057 VILLOVELA DE PIRON 202.0 1183 615 1212 33.0 GS-3089 MORLA DE LA VALDERIA 281.1 1369 1001 945 146.1 GS-3104 VILLAVERDE DE ARCAYOS 371.1 1146 1041 942 135.3 GS-3105 SANTERVAS DE CAMPOS 277.1 900 591 1107 28.8 GS-3124 MEDINA DE RIOSECO 908.4 790 450 1315 37.2 GS-3150 PARDAVE 223.8 1448 1356 787 220.8 GS-3818 RABAL 597.9 678 1328 1146 305.7 able 2.
Selected parameters for the calibration.
Parameter Units Lower bound Upper bound Description b i - 10 -5 D S - 10 -9
1 Fraction of D m where non-linear baseflow starts (see Eq. 2) W S - 10 -9
1 Fraction of the porosity of soil layer 3 where non-linear baseflow starts (see Eq. 2) D m mm/day 10 -9
30 Maximum baseflow (see Eq. 2) d m 0.1 0.9 Thickness of soil layer 2 able 3. Values of
NSE and r for the calibration and validation periods. Code
NSE cal
NSE val r r Code
NSE cal
NSE val r r R-2001 0.8520 0.4813 0.9494 0.8006 GS-3005 0.8307 0.5003 0.9156 0.8161 R-2011 0.9237 0.7606 0.9696 0.8947 GS-3016 0.8350 0.7798 0.9173 0.9126 R-2012 0.2723 0.3758 0.8983 0.8733 GS-3028 0.6290 0.6980 0.8320 0.8465 R-2013 0.8263 0.8141 0.9263 0.9152 GS-3035 0.5047 -2.1931 0.8850 0.7974 R-2014 0.8777 0.8736 0.9370 0.9347 GS-3041 0.7106 0.1784 0.8497 0.8344 R-2026 0.8865 0.7904 0.9645 0.8958 GS-3047 0.6277 0.4813 0.8148 0.7479 R-2027 0.8254 0.7003 0.9477 0.9278 GS-3049 0.7333 -0.0513 0.8738 0.8926 R-2028 0.7475 -0.3382 0.8812 0.8078 GS-3051 0.7743 0.7591 0.8909 0.9226 R-2030 0.6478 0.5780 0.8774 0.8730 GS-3057 0.7230 0.6515 0.8864 0.9176 R-2032 0.9409 0.8354 0.9705 0.9220 GS-3089 0.9116 0.9026 0.9564 0.9561 R-2036 0.6662 0.8070 0.8217 0.9177 GS-3104 0.8137 0.6694 0.9208 0.8985 R-2037 0.8222 0.3201 0.9114 0.7912 GS-3105 0.6204 0.3723 0.8032 0.8010 R-2038 0.9145 0.8301 0.9587 0.9188 GS-3124 0.7397 0.3129 0.8604 0.8246 R-2039 0.9521 0.2494 0.9760 0.8033 GS-3150 0.7909 0.8507 0.9118 0.9436 R-2042 0.8903 0.8309 0.9500 0.9559 GS-3818 0.8724 0.8975 0.9403 0.9720 R-2043 0.8322 0.8960 0.9237 0.9476 able 4.
Behavior of
NSE in the Monte Carlo simulations for assessing equifinality and the efficiency of the calibration algorithm.
Code
NSE cal
Number of simulations with
NSE > NSE cal – 0.05 Number of simulations with
NSE > NSE cal