A large-sample modelling approach towards integrating streamflow and evaporation data for the Spanish catchments

The simultaneous incorporation of streamflow and evaporation data into sensitivity analysis and calibration approaches has great potential to improve the representation of hydrologic processes in modelling frameworks. This work aims to investigate the capabilities of the Variable Infiltration Capacity (VIC) model in a large-sample application focused on the joint integration of streamflow and evaporation data for 189 headwater catchments located in Spain. The study has been articulated into three parts: (1) a regional sensitivity analysis for a total of 20 soil, routing, and vegetation parameters to select the most important parameters conducive to an adequate representation of the streamflow and evaporation dynamics; (2) a two-fold calibration approach against daily streamflow and monthly evaporation data based on the previous parameter selection for VIC; and (3) an evaluation of model performance based on a benchmark comparison against a well-established hydrologic model for the Spanish domain and a cross-validation test using multiple meteorological datasets to assess the generalizability of the calibrated parameters. The regional sensitivity analysis revealed that only two vegetation parameters – namely, the leaf area index and the minimum stomatal resistance – were sufficient to improve the performance of VIC for evaporation. These parameters were added to the soil and routing parameter during the calibration stage. Results from the two calibration experiments suggested that, while the streamflow performance remained close in both cases, the evaporation performance was highly improved if the objectives for streamflow and evaporation were combined into a single composite function during optimization. The VIC model outperformed the reference benchmark, and the independent meteorological datasets yielded a slight to moderate loss in model performance depending on the calibration experiment considered. Results from this investigation provide valuable insights into VIC parameter sensitivities, with a particular focus on large-sample applications, and highlight the importance of integrating multiple datasets into model calibration as a measure to reduce model equifinality.

1 Introduction

Large-sample hydrology (Addor et al., 2020) and large-scale hydrology (e.g. Bierkens, 2015; Wood et al., 2011) aim to promote the transferability of knowledge between regions and to assess the applicability of hydrologic models and theories at regional, continental, and global scales. Large-sample hydrology involves large sets (tens to thousands) of catchments, and its main focus is to provide generalizable knowledge of hydrological processes and models based on a large sample of catchments representing different hydroclimatic conditions (Addor et al., 2020). Similarly, large-scale hydrology relies on simulations from land surface models carried out at the so-called spatial hyper-resolution (> 1 km) to quantify and monitor the terrestrial water cycle at multiple scales (Bierkens, 2015; Bierkens et al., 2015; Wood et al., 2011).

Large-sample and large-scale hydrologic studies play an important role in supporting water resources planning and in quantifying hydrologic changes across scales in the context of a changing climate (Addor et al., 2020; Wood et al., 2011). The gap between both hydrologic fields is becoming increasingly reduced, and both can be considered to be two complementary approaches that attempt to provide a solid understanding of the spatial variability of hydrologic processes and to facilitate the intercomparison of model structures across climates (Addor et al., 2018, 2020). This is also manifested in the greater areas covered by hydrologic models (Beck et al., 2016, 2020), the development of gridded runoff observations (Gudmundsson and Seneviratne, 2016; Ghiggi et al., 2019), the tendency towards finer resolutions in land surface models (Bierkens, 2015; Wood et al., 2011), and the use of macroscale hydrologic models in large-sample studies (e.g. Mizukami et al., 2017; Newman et al., 2017; Rakovec et al., 2016a, b, 2019; Sepúlveda et al., 2022; Yeste et al., 2020, 2021).

From a data perspective, large-sample hydrology encompasses hydrologic studies founded on large-sample datasets of streamflow observations, hydrometeorological data, and hydroclimatic and landscape attributes (Addor et al., 2020; Kratzert et al., 2023). This includes investigations on extreme events (e.g. Blöschl et al., 2017; Do et al., 2017; Gudmundsson et al., 2019), climate change impacts (e.g. Marx et al., 2018; Melsen et al., 2018; Vormoor et al., 2015), variations in terrestrial water storage (e.g. Zhang et al., 2017), model evaluation and benchmarking (e.g. Aerts et al., 2022; Newman et al., 2017; Prieto et al., 2021, 2022; Rakovec et al., 2019; Yeste et al., 2020), data and modelling uncertainties (e.g. Beck et al., 2017; Coxon et al., 2015), parameter estimates during calibration (e.g. Beck et al., 2016, 2020; Mizukami et al., 2017; Rakovec et al., 2016a, b), and the transferability of parameters in space based on parameter regionalization techniques (e.g. Almeida et al., 2016; Beck et al., 2020; Pool et al., 2021; Prieto et al., 2019; Rakovec et al., 2019).

But over and above the extensive hydroclimatic characterization commonly provided in large-sample datasets, streamflow is considered to be a category of its own (Addor et al., 2020). Streamflow datasets are primarily based on individual contributions from national hydrologic services, which constitute the building blocks of continental and global streamflow repositories. The role of national water archives is of capital importance in this respect, and, ultimately, it is the international collaboration among national authorities worldwide which makes it possible to tackle this complex challenge (Addor et al., 2020). In this respect, the Model Parameter Estimation Experiment (MOPEX, Duan et al., 2006) was arguably the first open large-sample dataset, containing hydrologic information for 438 catchments within the CONUS (contiguous United States) domain. The MOPEX dataset was a highly significant contribution to the hydrologic community and constituted the basis for many hydrologic studies during the Prediction in Ungauged Basins (PUB) decade (Andreassian et al., 2007). However, MOPEX covers only up to 2003 and is no longer updated. A similar approach using up-to-date hydrologic data and landscape attributes materialized in the Catchments Attributes and MEteorology for Large-sample Studies (CAMELS) dataset (Addor et al., 2017) for the CONUS domain. The standards in CAMELS were also applied to develop large-sample datasets for other countries such as Chile (CAMELS-CL, Alvarez-Garreton et al., 2018), Great Britain (CAMELS-GB, Coxon et al., 2020), Brazil (CAMELS-BR, Chagas et al., 2020), Australia (CAMELS-AUS, Fowler et al., 2021) and Switzerland (CAMELS-CH, Höge et al., 2023), and most of them have been recently included in the Caravan dataset (Kratzert et al., 2023) as a step towards aggregating existing hydrologic information and producing a global dataset.

Large-sample hydrologic studies can benefit strongly from the integration of satellite remote sensing data into modelling frameworks in order to draw more robust conclusions on catchment functioning (Clark et al., 2017; Rakovec et al., 2016a, b, 2019; Yeste et al., 2020, 2021, 2023). In particular, the use of satellite-based algorithms to retrieve evaporation information represents an unprecedented opportunity to monitor the dynamics and the climate-driven changes in evaporative fluxes (Konapala et al., 2020; Koppa et al., 2022). Evaporation represents the second-largest component of the global water balance and is expected to increase as a consequence of global warming (IPCC, 2023). These changes can pose a challenge for future water security and water resources availability from regional to global scales (Lehner et al., 2019; Konapala et al., 2020; Koppa et al., 2022). Therefore, the integration of evaporation data into large-sample modelling approaches is a promising solution to calibrate and evaluate models for more than one hydrologic variable (traditionally streamflow; Dembélé et al., 2020a, b) and thus achieve a more reliable quantification of the water balance (Yeste et al., 2023).

This study aims to develop a hydrologic modelling framework to investigate the streamflow and evaporation dynamics for a large set of Spanish catchments. As part of the Iberian Peninsula, Spain constitutes a region where the effects of climate change are already noticeable and are expected to be much more pronounced by the end of the 21th century (IPCC, 2023). The Iberian Peninsula has been previously identified as a hotspot of climate change (Diffenbaugh and Giorgi, 2012; Vogel et al., 2021) and has manifested recurrent droughts and an increasing tendency towards arid conditions over the last decades, with a similar pattern for the Spanish catchments and a clear latitudinal gradient indicating greater aridity for the southern catchments (García-Valdecasas Ojeda et al., 2021a, b; Páscoa et al., 2017). From a hydrologic perspective, the Spanish catchments have undergone dramatic streamflow decreases during the last decades (Lorenzo-Lacruz et al., 2012, 2013), and evaporative fluxes play a dominant role in the water balance for the entire region (García-Valdecasas Ojeda et al., 2020a; Vicente-Serrano et al., 2014). These changes are expected to be exacerbated under climate change (García-Valdecasas Ojeda et al., 2020b, 2021a, b) and can pose an important threat to the future water planning and management in the country.

The large-sample modelling approach followed in this work will be focused on the Variable Infiltration Capacity (VIC) model (Liang et al., 1994, 1996), one of the most widely used hydrologic models in hydrologic studies (Addor and Melsen, 2019). The VIC model has been successfully implemented in many previous large-sample studies and large-scale applications (e.g. Melsen et al., 2018; Mizukami et al., 2017; Rakovec et al., 2019; Sepúlveda et al., 2022), including in Spain (Yeste et al., 2020, 2021), which makes it an excellent choice for the purpose of this investigation. The capabilities of VIC to integrate streamflow observations and satellite-based evaporation data will be thoroughly examined, and its performance will be compared against current modelling efforts providing the basis for water resource planning and management in Spain.

2 Study area and data

2.1 The Spanish catchments and streamflow dataset

This work is focused on a set of 189 headwater catchments defined for 94 reservoirs and 95 gauging stations belonging to the main river basin districts in Spain (Fig. 1a, Table 1). These catchments are representative of the hydroclimatic variability within the country, with mean annual precipitation ranging from 279 to 2183 mm yr⁻¹, mean annual runoff from 6 to 1821 mm yr⁻¹, mean annual potential evaporation from 770 to 2067 mm yr⁻¹, and runoff ratios from 0.02 to 0.98. Their physiography comprises areas ranging from 9 to 3825 km², with mean elevations from 147 to 1982 m and mean slopes from 4 to 100 m km⁻¹ (see also Fig. 3 and Table 3, which will be introduced later in Sect. 3.2). Streamflow observations for the Spanish catchments are monitored using automatic hydrological information systems (SAIHs, Sistemas Automáticos de Información Hidrológica) and by the official network of gauging stations (ROEA, Red Oficial de Estaciones de Aforo) and are estimated via a daily water balance of water storages and releases for reservoirs and using rating curves for gauging stations.

Figure 1 (a) Topographic boundaries of the 189 headwater catchments and the main river basin districts in Spain. The term “northern districts” refers to the districts of the Galicia coast, Miño-Sil, and Cantabria, which, for the sake of simplicity, were grouped using a common label. (b) Spatial distribution of the runoff ratio ($Q/P$). (c) Spatial distribution of the sum of the runoff and the evaporation ratios to precipitation ($(Q+E)/P$). Values in (b) and (c) were calculated for the study period and the hydroclimatic datasets considered in this work.

Table 1 Number of headwater reservoirs and gauging stations per river basin district included in this study.

The 189 study catchments were selected from the Integrated Network of Gauging Stations (SAIH-ROEA) dataset (https://www.miteco.gob.es/en/cartografia-y-sig/ide/descargas/agua/anuario-de-aforos.aspx, last access: 11 December 2024), a national archive of streamflow observations maintained and annually updated by the Spanish National Public Works Research Centre (CEDEX, Centro de Estudios y Experimentación de Obras Públicas). Similarly to Yeste et al. (2018, 2020, 2023), the study catchments were selected considering a maximum percentage of missing values in the streamflow series of 10 % for the period October 1990–September 2010, which was chosen as the study period in this work. Among those catchments, 171 (84 reservoirs + 87 gauging stations) presented less than 5 % of missing values, and 144 (80 reservoirs + 64 gauging stations) were below 1 %.

An exploratory data analysis of negative values during the study period was subsequently conducted and revealed that 47 reservoirs presented up to 5 % of negative estimates of daily streamflow and that 46 reservoirs presented more than 5 % of negative estimates, whereas all 95 gauging stations and only 1 reservoir did not present values below 0 (Fig. 2a). In addition, the percentage ratio of negative to positive values was calculated for each reservoir to quantify their relative importance, suggesting that negative values are close to 0 for reservoirs with less than 5 % of negative records and become more visible above 5 % (the median percentage ratio of negative to positive values for reservoirs with more than 5 % of negative records is 2.9 %).

Figure 2 Exploratory data analysis of the negative values in the daily time series of streamflow gathered from the SAIH-ROEA dataset for the 189 studied catchments during the study period. (a) Scatterplot representation of the percentage ratio of negative to positive values (y axis) against the percentage ratio of the number of negative values to the number of records (x axis) for each catchment. “R” denotes “reservoir”, and “GS” denotes “gauging station”. (b) Monthly distribution of the percentage ratio of the number of negative values per month to the total number of negative values. The blue line corresponds to the average hydrologic year of streamflow expressed as the monthly mean percentage of annual streamflow calculated over all the catchments.

One feasible explanation for the presence of negative values in the streamflow series of reservoirs in the SAIH-ROEA dataset is that inflow data are calculated by applying a daily water balance exclusively to water storages and releases. As opposed to gauging stations, where streamflow is derived from rating curves, the daily water balance in reservoirs can produce negative estimates of streamflow (i.e. inflow) when the variation in the storage is negative and when its magnitude is greater than the water releases. Given the prominent role of evaporative fluxes in the Spanish catchments, the presence of negative values is likely to happen on the warmest days of the hydrologic year, when streamflow is at its minimum or when null and open-water evaporation in reservoirs become relevant. To further test this hypothesis, the monthly distribution of negative values for all the catchments was calculated and compared against the average hydrologic year of streamflow (Fig. 2b). Results confirm that most negative values occur in summer months, emerging in late spring and extending to the beginning of the hydrologic year. The hypothesis is also supported by the location of reservoirs with more than 5 % of negative records as they are mostly concentrated in southern Spain (Fig. 2a) and are thus characterized by a warmer climate.

Unfortunately, incorporating the effect of evaporative fluxes into the water balance requires additional data that are not provided in the SAIH-ROEA dataset and that are neither publicly disclosed nor publicly available for all the reservoirs, such as pan-evaporation measurements and elevation area capacity curves. Hence, this course of action could not be adopted and was left out of the scope of this work. The effect of other potential driving factors on the negative records, such as seepage losses, is supposed to be minor in comparison to open-water evaporation as negative estimates tend to happen in summer for the southernmost reservoirs. Therefore, on the basis of this initial exploratory analysis, negative values were considered to be null, and all 189 headwater catchments were included in the modelling framework. The validity of such an assumption will be further discussed in the light of results from the modelling exercise.

2.2 Meteorological forcings and evaporation dataset

Daily precipitation and temperature data were collected from the Spanish PREcipitation At Daily scale (SPREAD, Serrano-Notivoli et al., 2017) and the Spanish TEmperature At Daily scale (STEAD, Serrano-Notivoli et al., 2019) datasets, two gridded products at ∼ 5 km resolution constructed by interpolating daily observations from a dense network of meteorological stations distributed across Spain. Monthly evaporation data were collected at 0.25° resolution from the Global Land Evaporation Amsterdam Model (GLEAM) version 3.5a (Martens et al., 2017; Miralles et al., 2011) and were remapped to the study catchments following a first-order conservative approach. GLEAM has shown less uncertainty compared to other satellite-based evaporation products (Xu et al., 2019) and has seen extensive use within hydrologic studies for the calibration and evaluation of hydrologic models (e.g. Bouaziz et al., 2021; Dembélé et al., 2020a, b; Koppa et al., 2019; Mei et al., 2023), particularly in studies involving data-scarce areas (e.g. Dembélé et al., 2020a, b; López López et al., 2017) and/or regions where evaporative fluxes are dominant (e.g. Dembélé et al., 2020a, b; Yeste et al., 2020, 2021, 2023).

Fig. 1b and c show the computed values of the runoff ratio ($Q/P$) and the sum of the runoff and evaporation ratio to precipitation ($(Q+E)/P$) for the study period and the hydroclimatic datasets described in this and the previous section. Two-thirds of the catchments manifested $Q/P$ estimates below 0.4 and were predominantly located in the southeastern sector of the country (Fig. 1b). Approximately 50 % of the catchments produced $(Q+E)/P$ estimates between 0.9 and 1.1, with negative imbalances ($(Q+E)/P<\mathrm{1}$) mostly corresponding to catchments with a low runoff ratio and with positive imbalances ($(Q+E)/P>\mathrm{1}$) towards the northwest (Fig. 1c). The effect of these imbalances will be thoroughly examined in light of the results.

3 Methods

3.1 The VIC model

The Variable Infiltration Capacity (VIC) model (Liang et al., 1994, 1996) version 4.2.d is a semi-distributed macroscale hydrologic model implementing water and energy balances within a gridded domain at daily and subdaily time steps. The model utilizes a three-layer soil profile to conceptualize runoff generation. Surface runoff, based on the Xinanjiang formulation (Zhao et al., 1980), occurs in the first two soil layers, while baseflow is generated in the bottom layer using the Arno equation (Franchini and Pacciani, 1991). The VIC model considers the subgrid variability in land uses through vegetation tiles. Evaporation in each grid cell is calculated as the sum of evaporation from bare soil, evaporation from the canopy interception, and transpiration and is constrained by the atmospheric demand for water vapour according to the Penman–Monteith equation. The model structure includes a snow model for accumulation and melting processes, employing snow bands to consider subgrid variability in topography, land uses, and precipitation, making it applicable across diverse geographical domains. While originally designed as a land surface scheme for Earth system models, VIC has seen extensive use globally as a hydrologic model and stands out for its widespread usage within the hydrologic community (Addor and Melsen, 2019).

The VIC model was applied with a gridded configuration at 0.05° resolution (∼ 5 km) and with a spin-up period of 10 years preceding the study period. Meteorological forcings were interpolated to the model resolution through a nearest-neighbour assignment. Notably, the VIC model lacks a consideration of the horizontal fluxes between adjacent grid cells, typically addressed by coupling a routing model. A gamma function was selected in this study to post-process the runoff simulations and account for the delay between runoff generation and catchment discharge (i.e. streamflow). The required soil and vegetation parameters to run VIC were collected at 1 km resolution from the following datasets: bulk density and soil textural classes were collected from SoilGrids1km (Hengl et al., 2014); porosity, saturated hydraulic conductivity, field capacity, and wilting point were collected from EU-SoilHydroGrids ver1.0 (Tóth et al., 2017); and land uses were collected from the UMD Global Land Cover Classification (Hansen et al., 2000), with associated vegetation parameters aligned with Global Land Data Assimilation System (GLDAS) specifications for VIC (Rodell et al., 2004). Soil parameters were regridded to the model resolution using a first-order conservative remapping, whereas land uses were kept at their original resolution as the subgrid variability in land uses is handled statistically within VIC.

3.2 Regional sensitivity analysis

Parameter sensitivities were analysed using the implementation of the regional sensitivity analysis (RSA) method of Hornberger and Spear (1981) in the SAFE Toolbox (Pianosi et al., 2015). RSA is based on a classification of model simulations into behavioural and non-behavioural according to one or more performance metrics and evaluates differences between parameter distributions corresponding to both classes. The RSA sensitivity index for a given parameter represents the maximum vertical distance between the cumulative distribution functions (CDFs) corresponding to the behavioural and non-behavioural classes, which is equivalent to the Kolmogorov–Smirnov distance statistic computed in the Kolmogorov–Smirnov test. Hence, the RSA sensitivity index ranges from 0 to 1, with values closer to 1 indicating a greater parameter sensitivity.

For each of the 189 catchments, the parametric space was explored by conducting a Monte Carlo simulation for 10 000 Latin hypercube samples (Iman and Conover, 1982) extracted from the parameter ranges of the 20 soil, vegetation, and routing parameters analysed in Yeste et al. (2023) and described in Table 2. Model performance during the Monte Carlo experiment was evaluated using the Nash–Sutcliffe efficiency (NSE), which can be defined as follows (e.g. Knoben et al., 2019):

$$\begin{array}{}\text{(1)}& \mathrm{NSE}=\mathrm{2}\mathit{\alpha }r-{\mathit{\alpha }}^{\mathrm{2}}-{\displaystyle \frac{(\mathit{\beta }-\mathrm{1}{)}^{\mathrm{2}}}{{\mathrm{CV}}_{\mathrm{obs}}^{\mathrm{2}}}},\end{array}$$

where r is the correlation coefficient between simulations and observations, α is the ratio of the standard deviations (i.e. $\mathit{\alpha }={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{sim}}}{{\mathit{\sigma }}_{\mathrm{obs}}}}$), β is the ratio for mean values (i.e. $\mathit{\beta }={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{sim}}}{{\mathit{\mu }}_{\mathrm{obs}}}}$), and CV_obs is the coefficient of variation of the observations. r, α, and β represent the dynamics, variability, and bias components for the comparison of simulations and observations, respectively, and CV_obs is model-independent.

Table 2 Parameters included in the RSA sensitivity analysis. As in Yeste et al. (2023), the subscript “f” corresponds to the VIC vegetation parameters that were modified using dimensionless multiplication factors.

RSA was applied to the NSE for the Monte Carlo simulations of daily streamflow (NSE(Q_d)) and monthly evaporation (NSE(E_m)) calculated for the study period, choosing the median NSE(Q_d) and median NSE(E_m) to classify behavioural and non-behavioural simulations. Similarly to Sepúlveda et al. (2022), the Spearman correlation coefficient (r_S) between the RSA sensitivity indices for NSE(Q_d) and NSE(E_m) and the physiographic and hydroclimatic characteristics defined in Table 3 and depicted in Fig. 3 was calculated. r_S measures how similar the spatial patterns of the parameter sensitivities and the selected attributes are, with their sign being indicative of a matching pattern (i.e. positive sign) or an opposite pattern (i.e. negative sign). These attributes were initially selected based on their ease of access for the study catchments and allowed for further investigation of the potential drivers of parameter sensitivities.

(Serrano-Notivoli et al., 2017)(Serrano-Notivoli et al., 2019)(Serrano-Notivoli et al., 2017)(Serrano-Notivoli et al., 2017)(Tóth et al., 2017)

Table 3 Definition of the physiographic and hydroclimatic characteristics analysed in this work as potential drivers of parameter sensitivities.

Figure 3 Spatial distribution of the physiographic and hydroclimatic characteristics analysed in this work as potential drivers of parameter sensitivities. These attributes are defined in Table 3.

The two most influential vegetation parameters in relation to any of the two performance metrics were finally selected for each catchment and were incorporated, together with the five soil parameters and the two routing parameters (SR parameterization hereafter), into the calibration stage. This selection will be referred to as the SR parameterization hereafter. Adding two extra VIC vegetation parameters to the SR parameterization is sufficient to improve the joint performance against streamflow and evaporation data according to previous research in Yeste et al. (2023).

3.3 Calibration and evaluation approach

A split-sample test (SST, Klemeš, 1986) was applied to calibrate and evaluate VIC considering two independent periods of equal duration belonging to the study period: a calibration period from October 1990 to September 2000 and an evaluation period from October 2000 to September 2010. A spin-up simulation of 10 hydrologic years preceding both the calibration and evaluation periods was conducted to provide initial states of model storages free from the effect of initial conditions (this strategy was also applied to the Monte Carlo simulation described in the previous section). The performance of the VIC model was evaluated through NSE(Q_d) and NSE(E_m) and their decomposition into correlation r, variability α, and bias β.

The calibration was performed using the shuffled-complex-evolution algorithm (SCE-UA) of Duan et al. (1994) and following a single-objective optimization approach for the nine selected parameters (five soil parameters, two routing parameters, and two vegetation parameters) to minimize a composite function that aggregates the performance metrics for streamflow and evaporation:

$$\begin{array}{}\text{(2)}& \mathrm{Minimize}\sqrt{w_Q\cdot (\mathrm{1}-\mathrm{NSE}(Q_{\mathrm{d}}){)}^{\mathrm{2}}+w_E\cdot (\mathrm{1}-\mathrm{NSE}\left(E_{\mathrm{m}}\right){)}^{\mathrm{2}}}.\end{array}$$

This problem minimizes the two-dimensional weighted Euclidean in relation to the ideal vector (1, 1) and belongs to the more general weighted-metric method used to minimize distances (see Yeste et al., 2023, for a detailed implementation of this method to integrate streamflow and evaporation data). In this work, two calibration experiments for different weights combinations in Eq. (2) were applied to the VIC model: firstly, a streamflow-only calibration (Q-only calibration hereafter) was performed by choosing w_Q=1 and w_E=0. Secondly, the model was calibrated for a weighted Euclidean distance (Q-E calibration hereafter) selecting two equal weights $w_Q=w_E=\mathrm{0.5}$. The case of two equal weights is equivalent to minimizing the pure Euclidean distance, that is, $w_Q=w_E=\mathrm{1}$, as equal weights do not affect the optimization problem in Eq. (2).

The results of the split-sample test were benchmarked against the streamflow and evaporation outputs from the Integrated System for Rainfall-Runoff Modeling (SIMPA) model (Estrela and Quintas, 1996; Alvarez et al., 2004) for the study catchments. SIMPA is a well-established hydrologic model calibrated for the Spanish catchments and updated yearly; it provides the foundation for water planning and management in the country (see https://www.miteco.gob.es/en/agua/temas/evaluacion-de-los-recursos-hidricos/evaluacion-recursos-hidricos-regimen-natural.html, last access: 11 December 2024) for additional information). SIMPA constitutes a semi-distributed implementation of the model proposed by Témez (1977), a four-parameter lumped conceptual model that has been extensively used in hydrologic studies in Spain (e.g. Bejarano et al., 2010; Marcos-Garcia et al., 2017; Yeste et al., 2018), and has evolved to include a snow routine and different hydrogeological modules. SIMPA simulations are run at a monthly time step and at 500 m spatial resolution and are commonly used for benchmarking purposes in the region (e.g. Pellicer-Martínez and Martínez-Paz, 2018; Suárez-Almiñana et al., 2020; Yeste et al., 2020). Therefore, given the absence of daily simulations for SIMPA, the performances of VIC and SIMPA were compared for the Nash–Sutcliffe efficiency of monthly streamflow (NSE(Q_m)) and NSE(E_m).

Furthermore, the effect of considering the negative values in the streamflow series of reservoirs to be null (see Sect. 2.2) was tested in an additional implementation of the Q-only calibration experiment for the 47 catchments with less than 5 % of negative records and the 46 catchments with more than 5 % of negative records (Fig. 2) after considering the negative values to be gaps. This strategy made it possible to evaluate the extent to which the performance of VIC was affected by the data-processing approach followed for the negative records.

Finally, a cross-validation test using multiple meteorological datasets was carried out to assess the generalizability of the calibrated parameters. Thus, the performance of the VIC model during the study period was further evaluated for the Q-only and Q-E calibration experiments using precipitation and temperature data gathered from a gridded dataset provided by the Spanish Meteorological Agency (AEMET; see https://www.aemet.es/en/serviciosclimaticos/cambio_climat/datos_diarios?w=2, last access: 11 December 2024) and E-OBS (Cornes et al., 2018).

4 Results

4.1 RSA sensitivity analysis

The RSA sensitivity indices for NSE(Q_d) and NSE(E_m) are depicted in Figs. 4 and 5, respectively. NSE(Q_d) sensitivities were mostly related to the five soil parameters and the two routing parameters (i.e. SR parameterization), with little or no influence from the vegetation parameters and no clear spatial pattern for the RSA indices (Fig. 4). Among these parameters, the highest sensitivities corresponded to d₂, rout₁, and rout₂, although b_i, D_S, W_S, and D_m were also influential in relation to the streamflow metric according to several local estimates.

Figure 4 Spatial distribution of the RSA sensitivity indices calculated for NSE(Q_d).

Figure 5 Spatial distribution of the RSA sensitivity indices calculated for NSE(E_m).

Contrarily to NSE(Q_d), NSE(E_m) scores were greatly influenced by the vegetation parameters (Fig. 5). In this case, the highest sensitivities corresponded to rmin_f and LAI_f and manifested a latitudinal gradient, with minimum sensitivities occurring for the northern catchments. This pattern is also noticeable, but to a lesser extent, for depth_1f, rarc_f, albedo_f, rough_f, and RGL_f. d₂ was revealed to be the most important soil parameter to NSE(E_m), and, as expected from the VIC model implementation, the routing parameters had a null effect due to the routing scheme being exclusively applied to post-processing the runoff simulations.

Figure 6 Spearman correlation coefficient (r_S) between the RSA sensitivity indices calculated for (a) NSE(Q_d) and (b) NSE(E_m) and the physiographic and hydroclimatic characteristics defined in Table 3 and depicted in Fig. 3. Full-size circles indicate statistically significant r_S estimates at the 95 % confidence level.

Figure 6 shows the Spearman correlation coefficient (r_S) between the RSA sensitivity indices calculated for NSE(Q_d) and NSE(E_m) and the physiographic and hydroclimatic characteristics in Fig. 3. The NSE(Q_d) sensitivities for the soil parameters presented an opposite pattern (i.e. negative correlations) compared to the mean annual precipitation; the aridity index; the normalized difference vegetation index (NDVI); and, to a lesser extent, the slope. A similar behaviour can be observed for various vegetation parameters such as rmin_f and LAI_f, although they were not identified as important to NSE(Q_d) (Fig. 4). Conversely, both routing parameters exhibited a matching pattern (i.e. positive correlations) compared to the previous four attributes, but their magnitude was higher for rout₁. Concerning NSE(E_m) sensitivities, the soil parameters produced positive correlations with respect to those characteristics, whereas the vegetation parameters still reflected an opposite pattern. The correlations for mean temperature and saturated hydraulic conductivity (K_S) became noticeable for the NSE(E_m) sensitivities and revealed an opposite pattern compared to the soil parameters and a matching pattern compared to most of the vegetation parameters.

The two most influential vegetation parameters in relation to any of the two performance metrics under study were, lastly, selected according to the values of the RSA index and were added to the SR parameterization during the calibration stage. Figure 7 indicates that LAI_f and rmin_f were the two most influential parameters for the vast majority of the catchments, with little influence from other vegetation parameters.

Figure 7 Frequency of the first-most influential and second-most influential vegetation parameters according to the RSA sensitivity indices calculated for NSE(Q_d) and NSE(E_m).

4.2 Split-sample test: calibration and evaluation

Figure 8 shows the spatial distributions of NSE(Q_d) and NSE(E_m) corresponding to the Q-only (Fig. 8a, c) and Q-E (Fig. 8b, d) calibration experiments for the calibration period. The relative gain or loss in model performance suggests that, while NSE(Q_d) remained similar for both calibrations (Fig. 8a, b), NSE(E_m) was highly improved after calibrating VIC against streamflow and evaporation data simultaneously (Fig. 8c, d).

Figure 8 Spatial distribution of NSE(Q_d) and NSE(E_m) for (a, c) the Q-only calibration and (b, d) the Q-E calibration during the calibration period.

Figure 9 CDFs of the NSE and its decomposition calculated for (a–d) daily streamflow, (e–h) monthly streamflow, and (i–l) monthly evaporation. Blue lines and green lines indicate the performance of VIC for the Q-only and Q-E calibration experiments, respectively, while grey lines correspond to the performance of SIMPA (note that the SIMPA simulations are only available at a monthly time step). Results for the calibration (evaluation) period are represented with solid (dashed) lines.

This can also be appreciated in the CDFs of NSE(Q_d) and NSE(E_m) depicted in Fig. 9a and i. The median NSE(Q_d) was close to 0.6 for both experiments during the calibration period (Fig. 9a), although the streamflow performance deteriorated slightly for the Q-E calibration. On the other hand, the median NSE(E_m) for the Q-E calibration was 0.67 during the calibration period (Fig. 9i), while the median NSE(E_m) for the Q-only calibration did not exceed 0. As for the evaluation period, the slight to moderate loss in model performance for NSE(Q_d) and NSE(E_m) was indicative of an acceptable implementation and an adequate predictive capability.

The decomposition of NSE(Q_d) revealed similar $r_{Q_{\mathrm{d}}}$ estimates for both calibration experiments (Fig. 9b) and ${\mathit{\alpha }}_{Q_{\mathrm{d}}}$ values generally below 1 (Fig. 9c). The ${\mathit{\beta }}_{Q_{\mathrm{d}}}$ distribution is approximately symmetric around the median for both calibrations but reflects a steeper CDF closer to 1 for the Q-only calibration (Fig. 9d). The NSE(E_m) improvement for the Q-E calibration is also evinced in its decomposition, with $r_{E_{\mathrm{m}}}$ being the component subject to the greatest enhancement (Fig. 9j). ${\mathit{\alpha }}_{E_{\mathrm{m}}}$ and ${\mathit{\beta }}_{E_{\mathrm{m}}}$ estimates are comparable for both experiments, with values slightly closer to 1 corresponding to the Q-E calibration (Fig. 9k, l), and point to a generalized overestimation of the variability and a slight underestimation of the bias, respectively.

The results of the split-sample test were subsequently benchmarked against the performance of the SIMPA model for monthly streamflow (Fig. 9e–h) and monthly evaporation (Fig. 9i–l). Both calibration experiments outperformed SIMPA in terms of NSE(Q_m) and its decomposition, and although the poor performance of SIMPA for monthly evaporation was comparable to that of the Q-only calibration, the Q-E calibration produced much higher NSE(E_m) estimates.

The effect of the imbalances in the hydroclimatic datasets on model performance during the split-sample test was evaluated for NSE(Q_d), NSE(E_m), and their decomposition for the complete study period considering the $(Q+E)/P$ ratio (Fig. 1c) as a signature of how far or close the water balance is from being closed. Results for the Q-only and Q-E calibration experiments are depicted in Fig. 10a–d and in Fig. 10e–h, respectively. Both NSE(Q_d) and NSE(E_m) estimates are higher for both experiments for catchments with $(Q+E)/P$ close to 1 (Fig. 10a, e). This effect is particularly noticeable for NSE(E_m) and the Q-only calibration as the model was not calibrated using evaporation data (Fig. 10a). As shown in Fig. 10d and h, both ${\mathit{\beta }}_{Q_{\mathrm{d}}}$ and ${\mathit{\beta }}_{E_{\mathrm{m}}}$ are closer to 1 for catchments with $(Q+E)/P$ close to 1, with a wider distribution of ${\mathit{\beta }}_{Q_{\mathrm{d}}}$ for the Q-E calibration due to the imbalances in the hydroclimatic datasets. These imbalances, however, do not have a marked effect on the dynamics (Fig. 10b, f) and the variability (Fig. 10c, g).

Figure 10 NSE and its decomposition calculated for daily streamflow and monthly evaporation for (a–d) the Q-only experiment and (e–h) the Q-E calibration considering the complete study period and using the $(Q+E)/P$ ratio as a signature of how far or close the water balance is from being closed.

Figure 11 Equifinality assessment based on the 1 %, 2 %, and 5 % best-performing simulations from the Monte Carlo experiment according to (1) NSE(Q_d) and (2) the Euclidean distance for NSE(Q_d) and NSE(E_m). The y axis represents the mean absolute deviation of NSE(Q_d) from the maximum NSE(Q_d) estimated during the Monte Carlo experiment.

In order to assess the role of incorporating evaporation data into model calibration in reducing equifinality (i.e. behavioural parameter combinations leading to similar high-performance estimates), the Monte Carlo simulation performed as part of the sensitivity analysis was leveraged, and the mean absolute deviation of NSE(Q_d) from the maximum NSE(Q_d) was calculated for every catchment considering the 1 %, 2 %, and 5 % best-performing simulations from the Monte Carlo experiment according to two criteria: (1) NSE(Q_d) itself and (2) the Euclidean distance for NSE(Q_d) and NSE(E_m). Figure 11 reflects that the deviations from the maximum NSE(Q_d) (1) become higher as the percentage of best performers considered increases and (2) are more pronounced when the best-performing criterion is based on the Euclidean distance. The first effect is a straightforward consequence of considering an increasing number of simulations to calculate the mean absolute deviation from the maximum NSE(Q_d). The second effect is an indicator of less equifinality as there are fewer parameter combinations yielding a performance close to the maximum NSE(Q_d).

Finally, the effect of handling the negative records in the streamflow series of 93 reservoirs (Fig. 2) was evaluated after considering them to be gaps and re-implementing the Q-only calibration experiment. Figure 12a–d shows almost identical distributions for NSE(Q_d) and its decomposition for the 47 reservoirs with less than 5 % of negative records, suggesting that the data-processing strategy applied to the negative values had a minimum impact on the performance of VIC, thus corroborating the validity of considering them to be null. This is also observable for the 46 reservoirs with more than 5 % of negative records (Fig. 12e–h), although slight differences became apparent for the bias component (i.e. ${\mathit{\beta }}_{Q_{\mathrm{d}}}$) as the number of records modified was greater.

Figure 12 CDFs of NSE(Q_d) and its decomposition for the Q-only calibration experiment in reservoirs presenting negative records of daily streamflow: (a–d) 47 reservoirs with less than 5 % of negative records and (e–h) 46 reservoirs with more than 5 % of negative records (see also Fig. 2). Dark-blue lines and light-blue lines indicate the performance of VIC after considering the negative values to be 0 and gaps, respectively. Results for the calibration (evaluation) period are represented with solid (dashed) lines.

4.3 Cross-validation test using multiple meteorological datasets

To cross-validate the results from the split-sample test, the $Q/P$ bias was firstly calculated as the $Q/P$ ratio difference between the calibrated VIC and the observations using SPREAD and STEAD, AEMET, and E-OBS data as the meteorological forcings of VIC for the complete study period (Fig. 13). $Q/P$ biases corresponding to the Q-only calibrated parameters were broadly in the range of ±0.1 for all the datasets (Fig. 13a–c), while the Q-E calibrated parameters produced increased deviations (Fig. 13d–f).

Figure 13 Spatial distribution of the $Q/P$ bias calculated as the $Q/P$ ratio difference between the calibrated model and the observations (see Fig. 1b) using meteorological data from SPREAD and STEAD, AEMET, and E-OBS to force VIC for the complete study period. (a–c) $Q/P$ bias corresponding to the Q-only calibrated parameters. (d–f) $Q/P$ bias corresponding to the Q-E calibrated parameters.

Results for SPREAD and STEAD and the Q-only calibration suggest that negative biases tend to be associated with higher $Q/P$ values and vice versa (see Fig. 13a and compare to Fig. 1b), whereas the $Q/P$ biases corresponding to the Q-E calibration display an opposite spatial distribution compared to that observed for the $(Q+E)/P$ values (see Fig. 13d and compare to Fig. 1c) and exhibit a high negative correlation ($r=-\mathrm{0.91}$). There is a predominance of negative $Q/P$ biases for both AEMET (Fig. 13b, e) and E-OBS (Fig. 13c, f), even though these differences became exacerbated when forcing VIC with E-OBS and reached values below −0.3 in many of the northern catchments for the Q-E calibration.

Figure 14 Distributions of NSE and its decomposition calculated for (a–d) daily streamflow and (e–h) monthly evaporation using the SPREAD and STEAD (SPST), AEMET, and E-OBS datasets for the complete study period. Blue, red, and purple boxplots (i.e. boxplots to the left in each pair group) correspond to the Q-only calibration experiment. Green, orange, and pink boxplots (i.e. boxplots to the right in each pair group) correspond to the Q-E calibration experiment.

The distributions of NSE(Q_d) and NSE(E_m), as well as their decomposition for each meteorological dataset, are depicted in Fig. 14. The performance attained for NSE(Q_d) using SPREAD and STEAD was closely followed by that obtained with AEMET but reflected a moderate deterioration in the case of E-OBS (Fig. 14a). The daily streamflow dynamics were generally well captured with the three datasets (Fig. 14b), particularly with SPREAD and STEAD and AEMET, although an enhanced underestimation of the variability, together with an underestimation of the bias, became noticeable for AEMET and E-OBS (Fig. 14c, d).

The performance of VIC for NSE(E_m), in turn, was similar for all the datasets and clearly demonstrated the effect of both calibration experiments (Fig. 14e). The monthly evaporation dynamics were well reproduced for the Q-E calibration in all cases (Fig. 14f), with estimates of the variability and the bias close to 1 (Fig. 14g, h).

5 Discussion

5.1 Parameter sensitivities

This study expands on previous investigations using VIC for the Duero River basin in Yeste et al. (2020, 2021, 2023) by involving the main river basin districts in Spain. The large-sample approach followed in this work allowed for drawing more robust conclusions on model realism and the hydrologic functioning of a wide range catchments representing the hydroclimatic variability within the country. Parameter sensitivities were quantified according to the RSA sensitivity analysis method for the 189 study catchments and the various soil, routing, and vegetation parameters indicated in Table 2. As in Yeste et al. (2023), sensitivities were calculated with respect to NSE(Q_d) and NSE(E_m), which were the performance metrics selected to evaluate the goodness of fit of VIC for streamflow and evaporation, respectively.

d₂ and the two routing parameters governing the gamma distribution function (i.e. rout₁ and rout₂) were identified as the most important parameters in relation to NSE(Q_d) (Fig. 4), highlighting the importance of applying a routing procedure to improve model performance for daily streamflow. The rest of the soil parameters were also identified as important to NSE(Q_d) and yielded comparable sensitivities to those uncovered in previous studies (e.g. Gou et al., 2020; Lilhare et al., 2020; Melsen and Guse, 2019; Mendoza et al., 2015; Yeste et al., 2020, 2023). The influence of the vegetation parameters on NSE(Q_d), however, was negligible, resembling the findings of Sepúlveda et al. (2022) for a large-sample application of VIC in Chile. The strong dependencies of the NSE(Q_d) sensitivities for the soil and vegetation parameters on mean annual precipitation, aridity index, and NDVI were manifested as either highly positive (i.e. a matching pattern) or highly negative (i.e. an opposite pattern) Spearman correlations (Fig. 6a), thus corroborating the interdependency between parameter sensitivities and climate variables found in Sepúlveda et al. (2022) for the Chilean catchments.

From a hydrologic perspective, the five soil parameters are responsible for the runoff generation process in VIC, and the negative correlations with respect to mean annual precipitation, aridity index, and NDVI in Fig. 6a indicate that they are more important for catchments characterized by a more arid climate. As the precipitation volume to be transformed into runoff is lower for such catchments, the role of these parameters becomes critical in modulating the runoff generation, whereas their effect is less relevant for catchments belonging to a more humid climate given the higher water availability. The generated runoff volume is subsequently routed to the catchment outlet according to a gamma-based unit hydrograph in a post-processing phase. The two routing parameters control the delay between runoff generation and catchment discharge (i.e. streamflow), and the positive correlations in Fig. 6a suggest that both parameters are important for the humid catchments as a consequence of the higher runoff volumes to be routed.

On the contrary, NSE(E_m) was found to be most sensitive to the vegetation parameters, with LAI_f and rmin_f being the most important vegetation parameters according to the RSA sensitivity indices (Figs. 5, 7). This is in line with the parameter sensitivities reported in Sepúlveda et al. (2022) and Yeste et al. (2023), suggesting that the VIC vegetation parameters have a significant potential to improve the representation of evaporative processes if included in model calibration. Among the soil parameters, d₂ was the most important parameter in relation to NSE(E_m).

In this case, the negative correlations reflected by LAI_f and rmin_f with respect to mean annual precipitation, aridity index, and NDVI in Fig. 6b denote a greater impact for arid catchments that is likely to be associated with the limiting effect on the evaporative processes entailed by a lower vegetation density. As for the soil parameters, the high NSE(E_m) sensitivities for d₂ could be related to the water uptake by vegetation in the root zone as it is directly affected by the thickness of the VIC soil layers. The positive correlations associated with the five soil parameters with respect to the previous characteristics are likely to be connected to the implementation of the closed water balance equation in VIC and manifest an opposite behaviour compared to that observed for NSE(Q_d). This effect was also appreciated in Yeste et al. (2020) for the sensitivities of the VIC soil parameters.

5.2 Model performance during the split-sample test and the cross-validation test

The large-sample application of the VIC model provided valuable insights into its performance for streamflow and evaporation in the 189 study catchments. The capability of VIC to produce satisfactory estimates of NSE(Q_d) and NSE(E_m) simultaneously was tested through a split-sample test encompassing two calibration experiments based on the weighted Euclidean distance definition for two objectives (Eq. 2), namely Q-only and Q-E calibration. While the Q-only calibration led to NSE(E_m) scores below 0 for more than half of the catchments, the Q-E calibration substantially increased the performance for NSE(E_m) and concomitantly produced NSE(Q_d) values closed to those corresponding to the Q-only calibration (Figs. 8, 9). Similar conclusions were reached in Yeste et al. (2023) according to five single-objective calibration experiments carried out for one test catchment located in the Guadalquivir River basin.

The benchmark comparison of VIC against the performance of SIMPA for NSE(Q_m) and NSE(E_m) clearly indicated an increased performance for both metrics when using VIC (Fig. 9). The monthly simulations from SIMPA stand as the greatest modelling effort available for the Spanish domain and play a fundamental role in supporting water resource planning at the national and river basin scales. Therefore, the implementation of VIC developed in this work constitutes an important leap forward in comparison with the SIMPA simulations as VIC was run at a finer temporal resolution and improved the individual and joint representation of streamflow and evaporation. Moreover, the performance of VIC for daily streamflow and monthly evaporation was similar to that reflected in large-sample applications using VIC over the CONUS domain in Mizukami et al. (2017) and Rakovec et al. (2019). The streamflow performance was also comparable to other modelling efforts involving the Duero River basin (Morán-Tejeda et al., 2014; Yeste et al., 2020, 2023), Tajo (Pellicer-Martínez and Martínez-Paz, 2018; Pellicer-Martínez et al., 2021), Guadalquivir (Yeste et al., 2018), Segura (Pellicer-Martinez and Martínez-Paz, 2015; Pellicer-Martínez et al., 2015), Júcar (Marcos-Garcia et al., 2017; Suárez-Almiñana et al., 2020), and the northern districts (Prieto et al., 2019, 2021, 2022).

An important role of incorporating multiple datasets into model calibration is to reduce equifinality as some of the parameter combinations that perform equally well when only one observational dataset is considered are rejected when more datasets are employed. This would be ideally addressed following a Pareto optimization approach for both NSE(Q_d) and NSE(E_m), as in Yeste et al. (2023). Pareto optimization is known for reducing equifinality as the model is calibrated for two or more objective functions simultaneously, resulting in a lower number of behavioural parameter sets (Efstratiadis and Koutsoyiannis, 2010). Although this work does not include a full Pareto optimization, it was still possible to address this issue for the Q-only and the Q-E calibration experiments as they represent two important solutions belonging to the Pareto front: the corner corresponding to the maximum NSE(Q_d) performance and the compromise solution for NSE(Q_d) and NSE(E_m) according to the Euclidean distance definition for both. The equifinality assessment was carried out for the best-performing simulations from the Monte Carlo experiment according to NSE(Q_d) and the Euclidean distance for NSE(Q_d) and NSE(E_m) and revealed a reduction in model equifinality for the Euclidean distance formulation (Fig. 11).

The amount of evaporation represents a key piece of hydrologic information that, together with streamflow and precipitation, allows for identifying gaining and losing catchments as they constitute an indirect measure of the inter-catchment groundwater flow (Liu et al., 2020). Gaining and losing catchments are therefore characterized by an unclosed water balance, which can potentially lead to an unrealistic representation of the partitioning of precipitation into streamflow and evaporation in the case of significant imbalances when using a (closed) water balance hydrologic model as VIC. The effect of the existing imbalances in the hydroclimatic datasets on model calibration was thoroughly examined in this work by analysing the relative gain or loss in model performance (Fig. 9), as well as the NSE decomposition into r, α, and β (Fig. 10) for the Q-only and Q-E calibration experiments. Although NSE(Q_d) values were slightly lower for the Q-E calibration, they remained similar for both calibration experiments, suggesting that the imbalances did not yield a significant deterioration in model performance for the study catchments. As indicated in Yeste et al. (2023), the β component is of capital importance from a water balance perspective when streamflow and evaporation data are integrated together and was revealed to be the component most sensitive to such imbalances, with a lesser influence on r and α.

As described in Sect. 2.2, the presence of negative values in the streamflow series of 93 reservoirs (Figs. 2, 12) is likely to be related to the indirect estimation of inflow data through a daily water balance of water storages and releases without considering the evaporative fluxes from the reservoir. In this respect, the initial exploratory data analysis for the negative records represents a call for action for future releases of the SAIH-ROEA dataset as it is not feasible to handle this issue on the basis of the current hydrologic information provided in this dataset. The effect of considering the negative values to be null was evaluated during the split-sample test for the Q-only calibration experiment to quantify their relative significance in terms of model performance. The distributions of NSE(Q_d) and its decomposition were virtually identical after considering the negative values to be null and gaps, suggesting that the simulated streamflow is also null or low during the warmest part of the hydrologic year.

Finally, the generalizability of the calibrated parameters was evaluated by means of a cross-validation test using meteorological information gathered from AEMET and E-OBS. As demonstrated in Yeste et al. (2023), the integration of streamflow and evaporation data into model calibration is ultimately subject to the law of conservation of mass and the magnitude of the imbalance stemming from merging three independent datasets of precipitation, streamflow, and evaporation. This limitation was thoroughly checked for the Q-only and Q-E calibration experiments and all the meteorological datasets, with results pointing to a higher $Q/P$ bias (Fig. 13) and a slight to moderate loss in model performance for the Q-E calibration (Fig. 14) as a consequence of calibrating VIC against streamflow and evaporation data simultaneously. The potential of the calibrated parameters, as well as the trade-off in model performance arising during the Q-E calibration experiment, will be further explored in future implementations of VIC for the Spanish catchments to produce seamless distributed parameter maps and Spain-wide simulations based on a fully gridded implementation. In addition, an assessment of model structure error can be potentially performed using the Framework for Understanding Structural Errors (FUSE, Clark et al., 2008) or the Structure for Unifying Multiple Modeling Alternatives (SUMMA, Clark et al., 2015) as a way to overcome the limitations of using one single model structure, as in this study.

6 Conclusions

In this work, a large-sample application of the VIC model was carried out for 189 headwater catchments belonging to the main river basin districts in Spain. The potential of combining streamflow and evaporation data into the hydrologic modelling exercise was explored for the sensitivity analysis stage, the calibration of the VIC parameters, and the evaluation of its performance for the streamflow and evaporation simulations. The key findings of this study can be summarized as follows:

• A regional sensitivity analysis allowed for the identification of the parameter sensitivities with respect to the selected metrics to evaluate the performance of VIC against daily streamflow and monthly evaporation data for all the study catchments. The soil and routing parameters were revealed to be the most important parameters in relation to the streamflow performance, while the influence from the vegetation parameters was negligible. The performance of VIC for evaporation was mostly controlled by the soil parameters and two of the vegetation parameters. The two routing parameters were identified to be important for the humid catchments as a consequence of the higher runoff volumes to be routed.

• The calibration of the VIC model was performed by following two single-objective calibration experiments: a calibration against daily streamflow data exclusively and a calibration against daily streamflow and monthly evaporation data simultaneously. The performance of VIC was assessed for two independent periods, suggesting that it is possible to achieve satisfactory adjustment to both hydrologic variables at the same time if their performance metrics are combined into a composite function based on a weighted Euclidean distance definition.

• A benchmark comparison was made between the performance of VIC and the monthly simulations from the SIMPA model, with the latter constituting the greatest modelling effort available to date for the Spanish domain. The VIC model led to an increased performance for both streamflow and evaporation compared to SIMPA, thus indicating promising potential for a fully gridded implementation of VIC in the future to carry out Spain-wide simulations.

• An additional evaluation of the performance of the VIC model was performed using meteorological observations from two independent gridded datasets in order to assess the generalizability of the calibrated parameters. The slight to moderate loss in model performance at this stage was subject to the calibration experiment under study, with a greater imbalance and a trade-off in model performance becoming apparent for the calibration against streamflow and evaporation data simultaneously.