# Comparison of Stationary and Dynamic Conceptual Models in a Mountainous and Data-Sparse Catchment in the South-Central Chilean Andes

Authors: Camila Toledo; Enrique Muñoz; Mauricio Zambrano-Bigiarini
Journal: Advances in Meteorology
 (2015)

---

## Abstract

In recent years, it has been documented that climatic variability influences hydrological processes; however, these influences, such as hydrologic dynamics, have not yet been incorporated into models, which have been assumed as stationary with regard to climatic conditions. In this study, the temporal variability of hydrological processes and their influence on the water balance of a mountainous and data-sparse catchment in Chile are observed and modeled through the comparison of a stationary (time-invariant parameters) and dynamic (time-variant parameters) model. Since conceptual models are the most adequate option for a data-scarce basin, a conceptual model integrated with the Monte Carlo Analysis Toolbox is used to perform the analyses. Simple analyses aimed at increasing the amount of information obtained from models were used. The General and Dynamic Identifiability Analyses were used to perform stationary and dynamic calibration strategies, respectively. As a result we concluded that the dynamic model is more robust than the stationary one. Additionally, DYNIA helped us to observe the temporal variability of hydrological processes. This analysis contributed to a better understanding of hydrological processes in a data-sparse Andean catchment and thus could potentially help reduce uncertainties in the outputs of hydrological models under scenarios of climate change and/or variability.


---

## Body

## 1. Introduction

Climate change has significant effects on the hydrological cycle; among some of the expected effects of climate change is an increase in the number (frequency) of floods and droughts on the planet [1]. This has a direct impact on the availability of water for agriculture, silviculture, and hydroelectricity, among other activities. Agriculture and hydroelectricity are the most important economic activities in south-central Chile, and they depend directly on the availability of water from Andean basins. Thus, it is necessary to improve the quality of and confidence in hydrological predictions and, at the same time, ascertain and better understand the behavior of hydrological processes in Andean basins. It is here that hydrological modeling comes into play, since it is one of the most used tools for water resources management and planning [2].Currently, there are many tools to support water resources planning and management. In particular, conceptual hydrological models have been widely used by the hydrological community [3, 4] to provide a better understanding of hydrological processes [5], to estimate long-term water availability [6], and to make projections of climate [7–9] and land use change [10, 11], among other uses. However, there is still a need not only to reproduce the past behavior of a basin, but also to evaluate the representativeness of the chosen model and its processes to assess the quality of the simulations for a given basin. This must be done not only in terms of a given measurement or objective function, but also by incorporating a range of parameter values in which the model exhibits physically meaningful behavior. Recently, there has been an increasing concern for addressing equifinality, identifiability, and quantifying the uncertainty associated with model simulations [12].The concept of equifinality [13] recognizes that several parameter sets might result in simulations with similar goodness of fit when compared to observations and therefore a single “best” parameter set cannot be identified [13]. However, under climate change conditions, there are several reasons to expect a change in model parameters as well [12]. Therefore, calibration parameters might change along different time periods or parts of the hydrograph [9].Experience has shown that in rainfall-runoff modeling one set of calibrated parameters (i.e., a stationary model) may not yield equally good approximations of all events or for different parts of the observed hydrograph [14]. Moreover, for long-term simulations, the optimum parameter set may be different for separate periods and also must change concurrently with changes in boundary conditions and process characteristics. Most scientists now acknowledge that the traditional assumption of stationarity in a hydrological model may no longer be adequate [15–17]. Therefore, to properly consider the effects of hydrological variability (nonstationarity), different parameter sets are needed to represent the dynamic behavior of the catchment. In other words, dynamic calibration strategies are required for better confidence in model predictions.Some initial studies attempting to address the aforementioned idea have been published [9, 18, 19]. They have shown through various methods why and how models must change in time and/or according to climatic and catchment properties.Recently, the concept of hydrological dynamics has been much discussed in the hydrological literature to be further incorporated into the conceptual formulation of a model [12, 15, 20]. This means that changes in the response of a basin caused by, for example, climate variability and change, are taken into account in the modeling stage. It is therefore important to find and study ways of detecting this behavior in new models, in order to have more reliable simulations of these changes and increase the efficiency and confidence in the management of water resources and systems.Traditionally, stationary model calibration (with fixed model parameters) is conducted through parameter optimization to maximize the agreement between observed and simulated streamflows. Usually, a single objective function or a combination of objective functions is used over the calibration period for better model representativeness and then for model validation the calibrated model is verified against the rest of the data. But is this method result adequate for long-term predictions? Is this method valid for a basin affected by climate variability such as the El Niño Southern Oscillation Phenomenon (ENSO)? Are these results reliable for predictions? As Luo et al. [21] mentioned, it is essential to carefully design model calibration procedures in order to ensure good performance of the hydrologic model for streamflow forecasting. Therefore, depending on the objective of the study and the predictability expected for the model, different calibration strategies will be needed.Most climate change impact studies run hydrological models calibrated against past observations with future climate inputs to predict future runoff characteristics. This approach is probably adequate for predicting changes in near-term runoff resulting from small changes in the climate inputs. However, these models cannot be extrapolated to predict changes in the future resulting from large changes in climate inputs and hydrological processes [17]. In other words, hydrological models calibrated using past observations may not adequately capture processes that may become more dominant in future and therefore a dynamic calibration strategy would probably be better in such cases.This study aims to improve the representativeness and simulation ability of models by incorporating observed hydrological dynamics in a data-sparse Andean basin highly affected by climate variability. In addition, observation of the process dynamics of the basin will improve our understanding of hydrological processes and deepen our knowledge of the hydrology of an Andean basin. We thus expect to improve the conceptual design of a hydrological model by including dynamic calibrated parameters within the model formulation. The General Identifiability Analysis (GIA) and Dynamic Identifiability Analysis (DYNIA) [22] are carried out to perform a stationary and dynamic model calibration, respectively. Finally, the representativeness of the stationary model is compared against the model calibrated with the dynamic approach.


## 2. Study Area

The study area (Figure1) corresponds to the upper Chillán catchment, measured at theChillán en Esperanza streamflow station (latitudes 36.76° to 36.87°S and longitudes 71.4 to 71.6°W). The drainage area of this catchment is 210.4 km², with elevations ranging from 435 m a.s.l. at the gauging station up to 3200 m a.s.l. at the Nevados de Chillán volcano in the Andes mountain range. The catchment is located in the Chilean Andes in the piedmont area. It exhibits a complex topography and is a data-sparse area without rainfall and temperature stations in the catchment (Figure 1).Figure 1
Digital elevation model and location of theChillán en Esperanza basin (top), seasonal rainfall distribution (bottom left), and seasonal variation curves of 5, 50, and 95% exceedance probability of theChillán en Esperanza river basin (bottom right).
(a)
(b)
(c)According to the rainfall rate described by Mardones et al. [23] at 36°45′ south, there exists a positive gradient of rainfall with the terrain altitude due to the so-called orographic effect (eastward, from the sea to the Andes mountains). Thus, rainfall increases at the western slopes of the Andean mountains, reaching the values of up to 3000 (mm) per year in the study area [24]. More than 70% of annual precipitation is concentrated in the rainy period (from May to September), and precipitation increases with elevation due to the orographic effect of the Andes mountain range [24, 25]. The annual rainfall is 2200 mm and the average seasonal temperature varies from 3.3°C in winter to 21.9°C in summer.The river has a pluvial regime (see Figure1), since more than 90% of the catchment is located below the annual average snow line of 2500 m a.s.l. [26]. Although the basin contains part of the western slope of the Chillán volcano, it does not present a marked snowmelt influence in spring or early summer (see the seasonal variation curve of 5% of probability of exceedance in Figure 1 where a second peak between October and January is not observed).Because the basin is located on the slope of a volcano, it presents high slopes and highly permeable soil layers. According to SERNAGEOMIN [27], most (over 80%) of the basin area is covered by a PPl3 substrate (from the Pliocene-Pleistocene period). This substrate is composed mainly of volcanic (extrusive) rocks formed when volcanic lava emerges, favoring a high degree of fracturing and permeable soil layers. Moreover, due to the old age of this substrate, weathering processes have started, increasing rock permeability. Therefore, the high slopes and permeable soil layers related to the volcano influence the hydrology of the Chillán basin with regard to rapid surface and groundwater runoff processes and permit a high amount of water to be stored in the basin each year. Moreover, these characteristics produce fast runoff responses in the basin (see Figure 1).TheChillán en Esperanza river basin is covered mostly by native forest with scarce anthropogenic alterations [28], mainly because the basin is located at the slope of the Chillán volcano in a very difficult-to-access area. Therefore, negligible land use and cover change have occurred during recent decades.


## 3. Hydrological Model Description

Because the main characteristics in mountainous environments are the complexity and spatial variability of hydroclimatic patterns, one would be tempted to choose a detailed, distributed, and physically based hydrological model to represent the hydrological behavior of a basin and the hydrological processes which occur in it [4]. But, in the case of a data-sparse basin, there will always be a need for more information (land use, soil map, vegetation map, etc.) than is available. Therefore, conceptual models are the simplest and most adequate option for data-sparse basins such as Andean basins. Thus, considering that conceptual pluvial models do not differ significantly in their conceptualization and capabilities [8], we chose the Muñoz model [29] to perform the analyses due to its simplicity and because it integrates the General and Dynamic Identifiability Analyses used in this study.The model used in this paper is the conceptual water balance model presented in Muñoz [29] and Muñoz et al. [30] (Figure 2). This model simulates the rainfall-runoff and snowmelt-runoff processes separately. Because the pluvial regime predominates in the study area, only the rainfall-runoff module has been used, and it will be briefly described here.Figure 2
Conceptual model.The lumped model considers the watershed as a double storage system: subsurface (SS) and deep storage (US). SS represents the water stored in the unsaturated soil layer as soil moisture, while US is the water covering the saturated soil layer. The main input data for the model are rainfall (PM) and potential evapotranspiration (PET). The model output is the total runoff (ETOT) at the watershed outlet, which includes both base-flow (ES) and surface runoff (EI). For a further description of the model’s characteristics, please refer to Muñoz [29] and Muñoz et al. [30].


## 4. Data

To estimate representative rainfall data of the catchment, 10 rain gauge stations located in the vicinity of the catchment were used as shown in Figure1. It is important to point out that Andean catchments are usually difficult-to-access areas, and therefore most (even all in this case) meteorological stations are located outside of the catchment and far from the highest areas of the catchments under study. Therefore, in studies of these catchments, it becomes a challenge to adequately represent their meteorological characteristics.Monthly data were provided by the National Water Agency (Dirección General de Aguas, or DGA) and were then interpolated using the inverse distance weighting method to obtain representative catchment values. This method was chosen due to its efficiency in interpolating rainfall from point measurements and due to the mountainous nature of the basin and its influence on the spatial variability of precipitation [31].Due to the absence of long-term and reliable temperature data from DGA stations, temperature time-series were obtained from the Center for Climatic Research, University of Delaware [32]. To estimate distributed temperatures in the catchment, the Thiessen polygon method was used. Temperature data were then used to compute potential evapotranspiration using the Thornthwaite method.The morphological characterization of the catchment and its area was defined using Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images of 1 arc-second resolution (30 m approximately) (Figure1).


## 5. Methodology

To compare a model which includes hydrological dynamics with one that does not, dynamic and stationary calibration strategies were used. Both calibration strategies were executed based on the regional sensitivity analysis (RSA) described by Hornberger and Spear [33]. For the stationary calibration, an RSA was performed over the whole period of analysis, while for the dynamic calibration the same analysis based on RSA was performed, but with a moving window.To study the sensitivity and temporal variability of the model parameters the Monte Carlo Analysis Toolbox (MCAT) [22] was combined with the Muñoz et al. [30] hydrological model.MCAT starts by running several simulations with different parameter sets sampled from physically meaningful ranges previously defined by the user. Then, it stores the parameter sets, the model outputs, and the values of the objective function(s) associated with each simulation, in order to evaluate model behavior afterwards [30]. The main analyses that were carried out in this study were (i) the General Identifiability Analysis (GIA) and (ii) Dynamic Identifiability Analysis (DYNIA) [22, 34]. Both analyses are based on the RSA and consist of a graphical consolidation of the results obtained from the MCAT simulations.On one hand, GIA gives an overview of the probability distribution function (pdf) and the cumulative pdf (cdf) of individual model parameters for the best 10% of the simulations based on the objective function (OF) computed over the whole simulation period. GIA relates the influence that each parameter has on the model outputs according to a predefined OF. This is achieved through a graphical analysis that compares the cumulative distribution curve of each parameter versus the value of each parameter for different intervals within the parameters’ range. Thus, a greater slope in the cumulative distribution curve indicates that the analyzed parameter has a strong influence on the model outputs. In the same way a constant slope indicates that the model is not sensitive to the parameter.On the other hand, DYNIA carries out a similar sensitivity analysis, but instead of using the entire simulation period, it considers a moving window of lengthn defined by the user. This window is moved through the available time-series in a stepwise fashion; that is, at time step t, it computes the residuals between t
-
n and t
+
n. The size of n is selected depending on the time length in which each parameter is influential and on the period that is to be studied. The measure of model performance is then used to condition the marginal parameter distribution at that particular time step. A threshold is applied to consider only the best performing parameter sets (the best 10% in this study) according to an appropriate support measure.In this study, windows of 3, 7, and 13 months were used, aiming at (i) obtaining the greatest variability associated with each parameter and objective function (3 months) and (ii) obtaining interannual variability in parameters (7 and 13 months).It should be noted that DYNIA processes a great amount of information that varies with modeling time, objective function, and number of parameters. Thus, its graphical representation consists of a grayscale graphic where the darkest shading at any time step represents the steepest slope of the cumulative distribution curve, that is, greater identifiability. Similarly, lighter shading represents a gentler slope of the cumulative distribution curve and therefore lower identifiability (see Figure3). A further description of GIA and DYNIA can be found in Wagener et al. [34, 35] and Wagener and Kollat [22].Figure 3
Dynamic Identifiability Analysis plots for the RMSE OF and for 3, 7, and 13 months and for parametersA (plots a, b, and c, resp.), C
max
⁡ (plots d, e, and f, resp.), and C
k (plots g, h, and I, resp.). The dashed lines show the 5 and 95% confidence limits, the black lines show the relative time-series of the streamflow data, and the blue lines show the time-series of highest identifiability for each time-window analyzed.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)In this study, GIA was used to constrain model parameters according to their identifiability and therefore to calibrate the model as stationary with time-invariant model parameters. DYNIA was used to identify the time variability of the most identifiable parameters and to calibrate a dynamic model. Then, further analyses were carried out based on the comparison of both calibration strategies.


### 5.1. Stationary Parameter Calibration

For the stationary parameter calibration, the method described by Muñoz et al. [30] was used as summarized below.(1)
Estimation of the adequate number of simulations: because the MCAT is based on random trials, it normally requires a large number of simulations to cover the widest spectrum of possible simulations. The number of Monte Carlo simulations was estimated via trial and error, where the stop criterion was met when the correlation (according to the Pearson correlation coefficient) between uncertainty bands (calculated as a linear correlation between the time-series of the upper and lower limits of the bands of uncertainty) of two different trials, but with same number of simulations, was equal to or greater than 0.999. Under this criterion, it was determined that the adequate number of simulations for this study is 10,000.
(2)
Identifiability of the input modification factors: because it is possible for the inputs to not be representative of the hydrometeorological characteristics of the basin, the model considers the possibility of modifying them in terms of amounts in order to ensure the closure of the long-term mass balance. To do this, factorsA and B serve the function of modifying the precipitation and potential evapotranspiration series, respectively. Therefore, the identifiability analysis was first performed on A and B. These factors were constrained according to their positive identifiability and through an iterative process of constraining and then repeating the 10,000 Monte Carlo simulations. The stop criterion was met when A and B did not show further identifiability and then the mean value of the parameter range was used as a calibrated factor.
(3)
Identifiability of the model parameters: due to the parameter interactions, various Monte Carlo simulations were carried out, in which the ranges of each were constrained in accordance with its observed identifiability from GIA. If the parameter showed a range which has a positive influence on the model outputs (i.e., positive identifiability), the parameter was constrained to that range. From this analysis the new range of an identifiable parameter was defined and then a new iteration (Monte Carlo simulations) was carried out.
A steeper slope in GIA indicates that the parameter is more observed (or repeated) in that range and therefore is identifiable in that range. Otherwise, a gentler slope indicates that the model parameter is not identifiable in such range and therefore was used to discard a nonidentifiable range of a model parameter. In this way, the range assigned to a particular parameter was constrained. The experiment was repeated in order to observe identifiability in the remaining parameters related to the aforementioned one. The exercise stopped when identifiability was not further observed within the range defined for each parameter.


### 5.2. Dynamic Parameter Calibration

An approach similar to the GIA calibration was used for the dynamic calibration, but instead of time-fixed, that is, stationary, parameter values, time-variant parameters were used, following the variability observed in the DYNIA. The calibration procedure is summarized below.(1)
To define the adequate number of random trials for the Monte Carlo analysis, the same analysis as for the stationary calibration was used.
(2)
Dynamic identifiability of the input modification factors: as described above, factorsA and B serve in modifying the precipitation and potential evapotranspiration series, respectively. Therefore, the Dynamic Identifiability Analysis was first performed on A and B starting with the one that showed the highest identifiability. The criterion to estimate the factor with higher identifiability was the highest slope in the cdf plot for the best 10% of the models (highest slope in the GIA plots).
Then, to obtain the time-varying factor identifiability, the same criterion as above was used but along a moving window. Therefore, the darkest shading for every time step was used to define the most identifiable value of the factor (or parameter) along time.
(3)
After ascertaining the time-series of the identifiable parameter, it was processed, and its seasonality (based on the mean monthly value of parameter) was recorded. Then, a seasonal parameter was incorporated into the model for model calibration instead of a fixed (stationary) parameter.
(4)
The same approach was then repeated for the remaining input modification factor, and then the same method was used for the dynamic calibration of the most identifiable model parameters.
It is important to point out that (i) the procedure consists of an iterative procedure where dynamic calibration is performed for each parameter one at a time and (ii) the parameters that did not show seasonality after DYNIA were calibrated as fixed according to the GIA.The analysis was performed using the full range of each parameter aimed at covering the full spectrum of variations that can present the different parameters along the moving window.


### 5.3. Objective Functions

For the analysis, three objective functions were used to analyze the influence of the hydrological processes in the different parts of the hydrograph used. The Root-Mean-Square Error (RMSE), transformed root-mean-square error (TRMSE), and runoff coefficient error (ROCE) were used (see (1)). These functions were selected because they are focused on the analysis of high flows, low flows, and the mass balance, respectively [36, 37]: (1)
RMSE
=
1
n
∑
j
=
1
n
S
j
-
O
j
2
,
TRMSE
=
1
n
∑
j
=
1
n
1
+
S
j
λ
-
1
λ
-
1
+
O
j
λ
-
1
λ
2
,
ROCE
=
abs
Q
s
-
P
-
-
Q
o
-
P
-
,where O
j and S
j are the observed and simulated flows, respectively, λ is related to the Box-Cox transformation for the data series with a value of 0.3, and Q
/
P represents the mean annual runoff coefficient and is calculated by the rate between the mean annual flow Q
- and the mean annual rainfall P
- for the simulated (s) or observed (o) flows.


### 5.4. Model Implementation

The hydrological model was implemented on a monthly time step. Due to the data availability, 17 years of hydrometeorological records were used for calibration (1961–1977) and validation (1978–1994). Additionally 10,000 simulations by iteration were carried out, drawing random values from each parameter according to their initial or constrained range. The initial range considered is shown in Table1.Table 1
Description of the model parameters and adjustment factors range for the rainfall-runoff module.
Parameter
Units
Description
Influence over
Parameter range
C
max
⁡
[—]
Maximum runoff coefficient when the subsurface layer is saturated
EI
0.1
0.7
P
Lim
[mm]
Limit of rainfall over which PPD exists
PPD
50
1000
D
[—]
Percentage of rainfall overP
Lim transformed into PPD
PPD
1
100
H
max
⁡
[mm]
Maximum capacity of the soil layer to retain water
C
max
⁡ and ER
200
700
PORC
[—]
Fraction ofH
max
⁡ that defines the soil water content restricting the evaporation processes
H
crit and ER
1
100
C
k
[—]
Subterranean runoff coefficient
ES
0.1
0.7
A
[—]
Adjustment factor of the precipitation data
PM
0.5
2.5
B
[—]
Adjustment factor of the potential evapotranspiration data
PET and ER
0.5
2.5
EI: direct runoff; PPD: direct deep percolation; ER; real evapotranspiration; ES: subterranean runoff; PM: precipitation; PET: potential evapotranspiration.Parameter values, simulated flows, and related OF values were stored for each iteration. Then, using this information, GIA and DYNIA were carried out.It is important to point out that the streamflow station was destroyed by a flood in 1995 and thus from then on no new records are available at this location.


## 5.1. Stationary Parameter Calibration

For the stationary parameter calibration, the method described by Muñoz et al. [30] was used as summarized below.(1)
Estimation of the adequate number of simulations: because the MCAT is based on random trials, it normally requires a large number of simulations to cover the widest spectrum of possible simulations. The number of Monte Carlo simulations was estimated via trial and error, where the stop criterion was met when the correlation (according to the Pearson correlation coefficient) between uncertainty bands (calculated as a linear correlation between the time-series of the upper and lower limits of the bands of uncertainty) of two different trials, but with same number of simulations, was equal to or greater than 0.999. Under this criterion, it was determined that the adequate number of simulations for this study is 10,000.
(2)
Identifiability of the input modification factors: because it is possible for the inputs to not be representative of the hydrometeorological characteristics of the basin, the model considers the possibility of modifying them in terms of amounts in order to ensure the closure of the long-term mass balance. To do this, factorsA and B serve the function of modifying the precipitation and potential evapotranspiration series, respectively. Therefore, the identifiability analysis was first performed on A and B. These factors were constrained according to their positive identifiability and through an iterative process of constraining and then repeating the 10,000 Monte Carlo simulations. The stop criterion was met when A and B did not show further identifiability and then the mean value of the parameter range was used as a calibrated factor.
(3)
Identifiability of the model parameters: due to the parameter interactions, various Monte Carlo simulations were carried out, in which the ranges of each were constrained in accordance with its observed identifiability from GIA. If the parameter showed a range which has a positive influence on the model outputs (i.e., positive identifiability), the parameter was constrained to that range. From this analysis the new range of an identifiable parameter was defined and then a new iteration (Monte Carlo simulations) was carried out.
A steeper slope in GIA indicates that the parameter is more observed (or repeated) in that range and therefore is identifiable in that range. Otherwise, a gentler slope indicates that the model parameter is not identifiable in such range and therefore was used to discard a nonidentifiable range of a model parameter. In this way, the range assigned to a particular parameter was constrained. The experiment was repeated in order to observe identifiability in the remaining parameters related to the aforementioned one. The exercise stopped when identifiability was not further observed within the range defined for each parameter.


## 5.2. Dynamic Parameter Calibration

An approach similar to the GIA calibration was used for the dynamic calibration, but instead of time-fixed, that is, stationary, parameter values, time-variant parameters were used, following the variability observed in the DYNIA. The calibration procedure is summarized below.(1)
To define the adequate number of random trials for the Monte Carlo analysis, the same analysis as for the stationary calibration was used.
(2)
Dynamic identifiability of the input modification factors: as described above, factorsA and B serve in modifying the precipitation and potential evapotranspiration series, respectively. Therefore, the Dynamic Identifiability Analysis was first performed on A and B starting with the one that showed the highest identifiability. The criterion to estimate the factor with higher identifiability was the highest slope in the cdf plot for the best 10% of the models (highest slope in the GIA plots).
Then, to obtain the time-varying factor identifiability, the same criterion as above was used but along a moving window. Therefore, the darkest shading for every time step was used to define the most identifiable value of the factor (or parameter) along time.
(3)
After ascertaining the time-series of the identifiable parameter, it was processed, and its seasonality (based on the mean monthly value of parameter) was recorded. Then, a seasonal parameter was incorporated into the model for model calibration instead of a fixed (stationary) parameter.
(4)
The same approach was then repeated for the remaining input modification factor, and then the same method was used for the dynamic calibration of the most identifiable model parameters.
It is important to point out that (i) the procedure consists of an iterative procedure where dynamic calibration is performed for each parameter one at a time and (ii) the parameters that did not show seasonality after DYNIA were calibrated as fixed according to the GIA.The analysis was performed using the full range of each parameter aimed at covering the full spectrum of variations that can present the different parameters along the moving window.


## 5.3. Objective Functions

For the analysis, three objective functions were used to analyze the influence of the hydrological processes in the different parts of the hydrograph used. The Root-Mean-Square Error (RMSE), transformed root-mean-square error (TRMSE), and runoff coefficient error (ROCE) were used (see (1)). These functions were selected because they are focused on the analysis of high flows, low flows, and the mass balance, respectively [36, 37]: (1)
RMSE
=
1
n
∑
j
=
1
n
S
j
-
O
j
2
,
TRMSE
=
1
n
∑
j
=
1
n
1
+
S
j
λ
-
1
λ
-
1
+
O
j
λ
-
1
λ
2
,
ROCE
=
abs
Q
s
-
P
-
-
Q
o
-
P
-
,where O
j and S
j are the observed and simulated flows, respectively, λ is related to the Box-Cox transformation for the data series with a value of 0.3, and Q
/
P represents the mean annual runoff coefficient and is calculated by the rate between the mean annual flow Q
- and the mean annual rainfall P
- for the simulated (s) or observed (o) flows.


## 5.4. Model Implementation

The hydrological model was implemented on a monthly time step. Due to the data availability, 17 years of hydrometeorological records were used for calibration (1961–1977) and validation (1978–1994). Additionally 10,000 simulations by iteration were carried out, drawing random values from each parameter according to their initial or constrained range. The initial range considered is shown in Table1.Table 1
Description of the model parameters and adjustment factors range for the rainfall-runoff module.
Parameter
Units
Description
Influence over
Parameter range
C
max
⁡
[—]
Maximum runoff coefficient when the subsurface layer is saturated
EI
0.1
0.7
P
Lim
[mm]
Limit of rainfall over which PPD exists
PPD
50
1000
D
[—]
Percentage of rainfall overP
Lim transformed into PPD
PPD
1
100
H
max
⁡
[mm]
Maximum capacity of the soil layer to retain water
C
max
⁡ and ER
200
700
PORC
[—]
Fraction ofH
max
⁡ that defines the soil water content restricting the evaporation processes
H
crit and ER
1
100
C
k
[—]
Subterranean runoff coefficient
ES
0.1
0.7
A
[—]
Adjustment factor of the precipitation data
PM
0.5
2.5
B
[—]
Adjustment factor of the potential evapotranspiration data
PET and ER
0.5
2.5
EI: direct runoff; PPD: direct deep percolation; ER; real evapotranspiration; ES: subterranean runoff; PM: precipitation; PET: potential evapotranspiration.Parameter values, simulated flows, and related OF values were stored for each iteration. Then, using this information, GIA and DYNIA were carried out.It is important to point out that the streamflow station was destroyed by a flood in 1995 and thus from then on no new records are available at this location.


## 6. Results and Discussion

### 6.1. Hydrological Processes Representation

Results of the GIA (not shown here) revealed that for RMSE, TRMSE, and ROCE the most identifiable parameters wereA, C
max
⁡, and C
k. Therefore, the DYNIA was carried out only on these three parameters.Figure3 shows, for example, the dynamic identifiability plots according to the RMSE. It is observed that the shorter the time-window used for the DYNIA is, the more variability the time-series shows (e.g., Figures 3(a), 3(d), and 3(g)), and therefore wider areas are covered by the 90% confidence limits (area between the dashed lines) and the parameter is less identifiable. These results are probably related to uncertainties and noise in data (rainfall, evapotranspiration, and streamflow data), which are reduced and averaged out with larger time-windows. But also with larger time-windows, less variability can be observed in parameter series and aggregated responses are observed mainly because with a larger time-window the DYNIA tends to approximate a GIA. Figures 4, 5, and 6 summarize the information from the DYNIA plots for A, C
max
⁡, and C
k and according to RMSE, TRMSE, and ROCE OFs, respectively. The boxes show the 25th and 75th percentiles (bottom and top blue lines) for each parameter, while the monthly average is shown by a red plus sign. Figures 5 and 6 do not show the boxes for the ROCE, because C
max
⁡ and C
k did not show identifiability on this OF.Figure 4
Seasonal variation of parameterA based on the DYNIA for the RMSE (left column), TRMSE (central column), and ROCE (right column) OFs and for 3 (top), 7 (middle), and 13 (bottom) months.Figure 5
Seasonal variation of parameterC
max
⁡ based on the DYNIA for the RMSE (left column), TRMSE (central column), and ROCE (right column) OFs and for 3 (top), 7 (middle), and 13 (bottom) months.Figure 6
Seasonal variation of parameterC
k based on the DYNIA for the RMSE (left column), TRMSE (central column), and ROCE (right column) OFs and for 3 (top), 7 (middle), and 13 (bottom) months.Figure3 shows that parameter A achieves values greater than 1 for the three OFs studied. A fulfills the function of correcting for orographic precipitation amounts. This parameter allows the long-term mass balance to be closed. In the cordillera and precordillera in south-central Chile, it is necessary to carry out this correction due to the absence of meteorological stations in the high zones of the Andes, due to which the increase in orographic precipitation is not well measured by the available instruments [25]. Similar analyses have been carried out by Muñoz et al. [30, 38, 39], and Zúñiga et al. [40].Figure3(a) shows that parameter A exhibits greater identifiability during high-flow periods (darker shaded areas and narrower ranges in high-flow periods). Figure 4 shows the seasonal variation of each parameter using boxes with the 25th and 75th percentiles (bottom and top blue lines) and the average (red plus sign) of the monthly time-series of each parameter. The range between the 25th and 75th percentiles is narrower for the rainy season (May–October), indicating that A is more identifiable in winter than in summer (December–March) (wider box) and suggesting that A is highly related to orography. Moreover, these results suggest that, for better representativeness of the rainfall amounts in the basin, the measured and interpolated data need to be amplified about 1.25 to 1.45 times in the rainy period.ParameterA is not identifiable (wider box) in summer because there is no rainfall that needs to be corrected. Similar results were obtained by Muñoz [4] in the Polcura River basin, close to the basin under study, in south-central Chile. These results suggest that the documented phenomenon of orography described by previous studies (e.g., [23, 25]) must be considered as part of the main hydrometeorological processes of the basin.ParameterC
max
⁡ represents the maximum runoff coefficient that is reached when the subsurface storage system becomes saturated, and it is related to the surface runoff processes.C
max
⁡ exhibits seasonal variability for RMSE (time-window of 3 and 7 months) and TRMSE (time-window of 3 months) (Figure 5). Comparing the dry and wet season, low parameter values are observed during summer (an average of about 0.2), which are related to less surface runoff, while the highest parameter values are observed during the rainy season (May–October), mostly with mean values between 0.3 and 0.4. However, the results suggest not only an influence of rainfall variability, but also the fact that geomorphology has a strong influence on surface runoff values and variability. More than 80% of the study area presents rocky and fractured formations [41]. This type of geology along with steep hill-slopes influences the maximum runoff coefficient, which reaches only 0.2 during the dry period and up to 0.4-0.5 in rainy periods. Despite the high rainfall amounts that can occur in the Chillán basin throughout the year and its high slopes, the maximum runoff coefficient is not as high as would be expected. This is probably related to the volcanic nature of the basin and the high permeability of the soil layers [27]. Therefore, these results suggest that surface runoff tends to be reduced by infiltration and percolation processes.In addition, the observed seasonality has a direct connection to the rainfall patterns of the zone, in which there are greaterC
max
⁡ values in periods with increased precipitation and therefore a greater degree of saturation in the surface layers of the basin.Regarding variability ofC
max
⁡ observed with the OF ROCE, a clear trend is not observed. Lower identifiability was observed in the dynamic identifiability plots (summarized by Figure 5) than in the RMSE and TRMSE, probably because ROCE is focused on the mass balance, and thus the runoff rates will be, relative to the other functions, much less identifiable than, for example, rainfall adjustment factor A. This suggests that the OF ROCE is only adequate to study the mass balance of the basin; after achieving it with the model, no further processes can be clearly identified.Figure6 shows the seasonal variation of parameter C
k, which represents the proportion of the subterranean storage system that becomes the base-flow. It is observed that C
k shows higher identifiability during low-flow periods (darker shadings in Figure 3(g)), indicating that it is more influential on model outputs during summer, where the only observed discharge is the base-flow. Similar results can be found in Figure 6 where the boxes become narrower during summer.Additionally,C
k tends to be larger during winter periods than during summer. This behavior is probably due to the type of geology described in the study area by SERNAGEOMIN [27] (fractured volcanic rocks), which favors the rapid movement of groundwater. Thus, as the basin recharges, the proportion of the deep storage that is discharged as base-flow also increases (in winter and early spring) increasing the C
k values, and then when the catchment starts to get dry the deep storage contribution decreases and therefore the proportion of groundwater storage that becomes the base-flow, that is, C
k, decreases.For the TRMSE, it is observed thatC
k values tend to be more bounded than the rest of the analyzed OFs, suggesting that the TRMSE is adequate for studying base-flow dynamics.These results are in line with the expected conceptual behavior of a volcanic basin with high slopes and highly permeable and fractured soil layers. The DYNIA is an adequate method to observe the main hydrological process dynamics and to perform a dynamic calibration strategy because it recovers the natural variability of identifiable processes.


### 6.2. Stationary against Dynamic Calibration

To compare both strategies of calibration, OF values were compared (Table2). Also the Nash–Sutcliffe efficiency index (NSE) related to each calibration strategy (stationary or dynamic) and related to each OF was included to better discuss the model results.Table 2
OF and NSE-related values for the stationary and dynamic strategies of calibration.
Stationary
Dynamic
Calibration
Validation
Calibration (1961–1977)
Validation (1978–1994)
1961–1977
1978–1994
3 months
7 months
13 months
3 months
7 months
13 months
RMSE
4.327
4.694
5.343
4.464
4.475
6.173
4.754
4.455
TRMSE
0.483
0.606
0.571
0.503
0.489
0.676
0.605
0.585
ROCE
0.001
0.006
0.002
0.000
0.000
0.010
0.007
0.008
NSE-related
RMSE
0.840
0.860
0.774
0.841
0.843
0.785
0.867
0.872
TRMSE
0.840
0.860
0.740
0.827
0.843
0.779
0.867
0.869
ROCE
0.840
0.850
0.728
0.770
0.749
0.734
0.757
0.758Table2 shows the OF values for the stationary and dynamic calibration of 3, 7, and 13 months. For the RMSE and TRMSE, the stationary calibration shows slightly better values than the dynamic approach, which is counterintuitive. However, for the 13-month time-window, the dynamic calibration shows a better RMSE value during the validation period and also shows a slightly smaller value than that for calibration. Something similar is observed for the TRMSE. The calibration-validation varies from 0.483 to 0.606 with a stationary calibration while it varies from 0.503 to 0.605 for a time-window of 7 months and from 0.489 to 0.585 for a time-window of 13 months. This indicates that although with dynamic calibration better values are not achieved, the validation tends to be better (smaller). In other words, the difference between calibration and validation values tends to be smaller in a dynamic than in a stationary calibration.A stationary model calibration seeks to maximize the agreement between observed and simulated streamflows in the evaluation period; therefore, all parameter values are defined to maximize a measure of evaluation (i.e., an objective function). On the other hand, dynamic calibration is performed to maximize the measure of evaluation as well as hydrological process representativeness. These processes (represented by parameters) have exhibited interannual variability and seasonality. The observed seasonality was realistic for processes such as orography, surface runoff, and groundwater runoff in a mountainous Chillán river basin.As many authors have discussed (e.g., [12, 17, 18]), hydrological processes vary in time, and therefore a dynamic model should better simulate basin behavior and its processes. The narrower range observed for the calibration-validation analysis suggests that a dynamic model tends to be better (more robust) than a stationary model and tends to better approximate realistic basin behavior. Moreover, the dynamic calibration also contributed to a better understanding of the system under study and was a very useful option in poorly understood basins such as Chilean Andean basins.When we compare the ROCE no benefits are observed in the dynamic calibration. This is probably due to the OF which is only focused on the mass balance and therefore does not take into account the hydrological processes of the basin and the exact shape of the hydrograph while a good water balance is achieved. Further studies using ROCE for studying the mass balance combined with a secondary OF for comparing the shape of the hydrograph can be recommended for the future.The NSE-related values show the same behavior as discussed above. Validation NSE values tend to be slightly better with the dynamic calibration (7 and 13 months), suggesting that the model representativeness is improved by including the dynamic variation of parameters in the model conceptualization and its equations.This study has shown for an Andean basin that, using a time-varying approach for model conceptualization, a better (more robust and representative) model can be obtained, which is in full agreement with the results of Oudin et al. [18].


## 6.1. Hydrological Processes Representation

Results of the GIA (not shown here) revealed that for RMSE, TRMSE, and ROCE the most identifiable parameters wereA, C
max
⁡, and C
k. Therefore, the DYNIA was carried out only on these three parameters.Figure3 shows, for example, the dynamic identifiability plots according to the RMSE. It is observed that the shorter the time-window used for the DYNIA is, the more variability the time-series shows (e.g., Figures 3(a), 3(d), and 3(g)), and therefore wider areas are covered by the 90% confidence limits (area between the dashed lines) and the parameter is less identifiable. These results are probably related to uncertainties and noise in data (rainfall, evapotranspiration, and streamflow data), which are reduced and averaged out with larger time-windows. But also with larger time-windows, less variability can be observed in parameter series and aggregated responses are observed mainly because with a larger time-window the DYNIA tends to approximate a GIA. Figures 4, 5, and 6 summarize the information from the DYNIA plots for A, C
max
⁡, and C
k and according to RMSE, TRMSE, and ROCE OFs, respectively. The boxes show the 25th and 75th percentiles (bottom and top blue lines) for each parameter, while the monthly average is shown by a red plus sign. Figures 5 and 6 do not show the boxes for the ROCE, because C
max
⁡ and C
k did not show identifiability on this OF.Figure 4
Seasonal variation of parameterA based on the DYNIA for the RMSE (left column), TRMSE (central column), and ROCE (right column) OFs and for 3 (top), 7 (middle), and 13 (bottom) months.Figure 5
Seasonal variation of parameterC
max
⁡ based on the DYNIA for the RMSE (left column), TRMSE (central column), and ROCE (right column) OFs and for 3 (top), 7 (middle), and 13 (bottom) months.Figure 6
Seasonal variation of parameterC
k based on the DYNIA for the RMSE (left column), TRMSE (central column), and ROCE (right column) OFs and for 3 (top), 7 (middle), and 13 (bottom) months.Figure3 shows that parameter A achieves values greater than 1 for the three OFs studied. A fulfills the function of correcting for orographic precipitation amounts. This parameter allows the long-term mass balance to be closed. In the cordillera and precordillera in south-central Chile, it is necessary to carry out this correction due to the absence of meteorological stations in the high zones of the Andes, due to which the increase in orographic precipitation is not well measured by the available instruments [25]. Similar analyses have been carried out by Muñoz et al. [30, 38, 39], and Zúñiga et al. [40].Figure3(a) shows that parameter A exhibits greater identifiability during high-flow periods (darker shaded areas and narrower ranges in high-flow periods). Figure 4 shows the seasonal variation of each parameter using boxes with the 25th and 75th percentiles (bottom and top blue lines) and the average (red plus sign) of the monthly time-series of each parameter. The range between the 25th and 75th percentiles is narrower for the rainy season (May–October), indicating that A is more identifiable in winter than in summer (December–March) (wider box) and suggesting that A is highly related to orography. Moreover, these results suggest that, for better representativeness of the rainfall amounts in the basin, the measured and interpolated data need to be amplified about 1.25 to 1.45 times in the rainy period.ParameterA is not identifiable (wider box) in summer because there is no rainfall that needs to be corrected. Similar results were obtained by Muñoz [4] in the Polcura River basin, close to the basin under study, in south-central Chile. These results suggest that the documented phenomenon of orography described by previous studies (e.g., [23, 25]) must be considered as part of the main hydrometeorological processes of the basin.ParameterC
max
⁡ represents the maximum runoff coefficient that is reached when the subsurface storage system becomes saturated, and it is related to the surface runoff processes.C
max
⁡ exhibits seasonal variability for RMSE (time-window of 3 and 7 months) and TRMSE (time-window of 3 months) (Figure 5). Comparing the dry and wet season, low parameter values are observed during summer (an average of about 0.2), which are related to less surface runoff, while the highest parameter values are observed during the rainy season (May–October), mostly with mean values between 0.3 and 0.4. However, the results suggest not only an influence of rainfall variability, but also the fact that geomorphology has a strong influence on surface runoff values and variability. More than 80% of the study area presents rocky and fractured formations [41]. This type of geology along with steep hill-slopes influences the maximum runoff coefficient, which reaches only 0.2 during the dry period and up to 0.4-0.5 in rainy periods. Despite the high rainfall amounts that can occur in the Chillán basin throughout the year and its high slopes, the maximum runoff coefficient is not as high as would be expected. This is probably related to the volcanic nature of the basin and the high permeability of the soil layers [27]. Therefore, these results suggest that surface runoff tends to be reduced by infiltration and percolation processes.In addition, the observed seasonality has a direct connection to the rainfall patterns of the zone, in which there are greaterC
max
⁡ values in periods with increased precipitation and therefore a greater degree of saturation in the surface layers of the basin.Regarding variability ofC
max
⁡ observed with the OF ROCE, a clear trend is not observed. Lower identifiability was observed in the dynamic identifiability plots (summarized by Figure 5) than in the RMSE and TRMSE, probably because ROCE is focused on the mass balance, and thus the runoff rates will be, relative to the other functions, much less identifiable than, for example, rainfall adjustment factor A. This suggests that the OF ROCE is only adequate to study the mass balance of the basin; after achieving it with the model, no further processes can be clearly identified.Figure6 shows the seasonal variation of parameter C
k, which represents the proportion of the subterranean storage system that becomes the base-flow. It is observed that C
k shows higher identifiability during low-flow periods (darker shadings in Figure 3(g)), indicating that it is more influential on model outputs during summer, where the only observed discharge is the base-flow. Similar results can be found in Figure 6 where the boxes become narrower during summer.Additionally,C
k tends to be larger during winter periods than during summer. This behavior is probably due to the type of geology described in the study area by SERNAGEOMIN [27] (fractured volcanic rocks), which favors the rapid movement of groundwater. Thus, as the basin recharges, the proportion of the deep storage that is discharged as base-flow also increases (in winter and early spring) increasing the C
k values, and then when the catchment starts to get dry the deep storage contribution decreases and therefore the proportion of groundwater storage that becomes the base-flow, that is, C
k, decreases.For the TRMSE, it is observed thatC
k values tend to be more bounded than the rest of the analyzed OFs, suggesting that the TRMSE is adequate for studying base-flow dynamics.These results are in line with the expected conceptual behavior of a volcanic basin with high slopes and highly permeable and fractured soil layers. The DYNIA is an adequate method to observe the main hydrological process dynamics and to perform a dynamic calibration strategy because it recovers the natural variability of identifiable processes.


## 6.2. Stationary against Dynamic Calibration

To compare both strategies of calibration, OF values were compared (Table2). Also the Nash–Sutcliffe efficiency index (NSE) related to each calibration strategy (stationary or dynamic) and related to each OF was included to better discuss the model results.Table 2
OF and NSE-related values for the stationary and dynamic strategies of calibration.
Stationary
Dynamic
Calibration
Validation
Calibration (1961–1977)
Validation (1978–1994)
1961–1977
1978–1994
3 months
7 months
13 months
3 months
7 months
13 months
RMSE
4.327
4.694
5.343
4.464
4.475
6.173
4.754
4.455
TRMSE
0.483
0.606
0.571
0.503
0.489
0.676
0.605
0.585
ROCE
0.001
0.006
0.002
0.000
0.000
0.010
0.007
0.008
NSE-related
RMSE
0.840
0.860
0.774
0.841
0.843
0.785
0.867
0.872
TRMSE
0.840
0.860
0.740
0.827
0.843
0.779
0.867
0.869
ROCE
0.840
0.850
0.728
0.770
0.749
0.734
0.757
0.758Table2 shows the OF values for the stationary and dynamic calibration of 3, 7, and 13 months. For the RMSE and TRMSE, the stationary calibration shows slightly better values than the dynamic approach, which is counterintuitive. However, for the 13-month time-window, the dynamic calibration shows a better RMSE value during the validation period and also shows a slightly smaller value than that for calibration. Something similar is observed for the TRMSE. The calibration-validation varies from 0.483 to 0.606 with a stationary calibration while it varies from 0.503 to 0.605 for a time-window of 7 months and from 0.489 to 0.585 for a time-window of 13 months. This indicates that although with dynamic calibration better values are not achieved, the validation tends to be better (smaller). In other words, the difference between calibration and validation values tends to be smaller in a dynamic than in a stationary calibration.A stationary model calibration seeks to maximize the agreement between observed and simulated streamflows in the evaluation period; therefore, all parameter values are defined to maximize a measure of evaluation (i.e., an objective function). On the other hand, dynamic calibration is performed to maximize the measure of evaluation as well as hydrological process representativeness. These processes (represented by parameters) have exhibited interannual variability and seasonality. The observed seasonality was realistic for processes such as orography, surface runoff, and groundwater runoff in a mountainous Chillán river basin.As many authors have discussed (e.g., [12, 17, 18]), hydrological processes vary in time, and therefore a dynamic model should better simulate basin behavior and its processes. The narrower range observed for the calibration-validation analysis suggests that a dynamic model tends to be better (more robust) than a stationary model and tends to better approximate realistic basin behavior. Moreover, the dynamic calibration also contributed to a better understanding of the system under study and was a very useful option in poorly understood basins such as Chilean Andean basins.When we compare the ROCE no benefits are observed in the dynamic calibration. This is probably due to the OF which is only focused on the mass balance and therefore does not take into account the hydrological processes of the basin and the exact shape of the hydrograph while a good water balance is achieved. Further studies using ROCE for studying the mass balance combined with a secondary OF for comparing the shape of the hydrograph can be recommended for the future.The NSE-related values show the same behavior as discussed above. Validation NSE values tend to be slightly better with the dynamic calibration (7 and 13 months), suggesting that the model representativeness is improved by including the dynamic variation of parameters in the model conceptualization and its equations.This study has shown for an Andean basin that, using a time-varying approach for model conceptualization, a better (more robust and representative) model can be obtained, which is in full agreement with the results of Oudin et al. [18].


## 7. Conclusions

Stationary (with fixed model parameters) and dynamic (time-varying parameters) calibration strategies were performed based on the General and Dynamic Identifiability Analyses, respectively.The General Identifiability Analysis was performed to maximize the agreement between simulated and observed streamflows, while the Dynamic Identifiability Analysis also aimed to better represent the variability of basin processes.The DYNIA proved to be a very useful option for studying data-sparse basins such as theChillán en Esperanza river basin. DYNIA, combined with a conceptual model, was an adequate option for obtaining realistic information about hydrological process variability in the basin. Therefore, DYNIA contributed to a better understanding of a data-sparse mountainous basin.Time variability of processes such as the orography, surface runoff, and base-flow was observed and moreover modeled.Dynamic calibration proved to be better than stationary. Although stationary calibration resulted in better OFs in calibration, the dynamic model presented better validation values. Additionally, the differences between calibration and validation were lower with the dynamic calibration strategy, suggesting that the dynamic conceptualization and formulation of a hydrological model are better (more robust) and more realistic than in a stationary model.The results shown in the present study suggest that the incorporation of hydrological dynamics in models must be considered in the modeling stage, not only for a better measurement of the OFs but also because it is the most adequate option for representing the hydrological processes that predominate in the basin and the variability that they could have in time.This study helped to improve our understanding of the hydrological processes related to groundwater runoff, surface runoff, and the orographic effect produced in a data-sparse Andean catchment. Regarding the seasonality of these processes, it was observed that they could be directly related to fluctuations in model inputs and are probably related to climate variability and change.Further studies which include larger scales variability such as climate change trends and variability, caused by, for example, the El Niño Southern Oscillation phenomenon, are recommended.It is expected that including these dynamics in modeling will help make forecasts and predictions less uncertain and more representative of the hydrological processes in the study area.



---

*Source: 526158-2015-11-02.xml*