Nonstationary Annual Maximum Flood Frequency Analysis Using a Conceptual Hydrologic Model with Time-Varying Parameters

Abstract

Recent evidence of the impact of watershed underlying conditions on hydrological processes have made the assumption of stationarity widely questioned. In this study, the temporal variations of frequency distributions of the annual maximum flood were investigated by continuous hydrological simulation considering nonstationarity for Weihe River Basin (WRB) in northwestern China. To this end, two nonstationary versions of the GR4J model were introduced, where the production storage capacity parameter was regarded as a function of time and watershed conditions (e.g., reservoir storage and soil-water conservation land area), respectively. Then the models were used to generate long-term runoff series to derive flood frequency distributions, with synthetic rainfall series generated by a stochastic rainfall model as input. The results show a better performance of the nonstationary GR4J model in runoff simulation than the stationary version, especially for the annual maximum flow series, with the corresponding NSE metric increasing from 0.721 to 0.808. The application of the nonstationary flood frequency analysis indicates the presence of significant nonstationarity in the flood quantiles and magnitudes, where the flood quantiles for an annual exceedance probability of 0.01 range from 4187 m³/s to 8335 m³/s for the past decades. This study can serve as a reference for flood risk management in WRB and possibly for other basins undergoing drastic changes caused by intense human activities.

1. Introduction

The occurrence of hydrological extremes has always been a challenge for human civilization. In the last decades, the exposure of critical infrastructures to flood has brought the attention of the scientific community on regional and national scales [1,2]. The Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) highlighted the increase in extreme precipitation and the associated increase in the frequency and magnitude of river floods in many regions, including Asia, central Africa, and western Europe [3]. The change in the magnitude and frequency of these extremes has questioned the assumption of stationarity, on which most flood frequency analysis (FFA) studies were based [4,5]. Failure to account for nonstationarity, especially nowadays when the change of flood regimes due to climate change or/and human activities can no longer be ignored, may lead to underestimation or overestimation of the design flood value.

Design flood value associated with a specific return period (e.g., 100-year that is most commonly used in relevant studies) is traditionally estimated based on flood frequency analysis through statistical approaches [6]. Recently, with the concept that stationarity is dead [4], various methodologies using a probabilistic model of flood frequency considering nonstationarity have been introduced [7,8,9,10,11,12,13,14,15]. Khaliq et al. [16] provide a comprehensive review of these approaches, including the incorporation of trends in distribution parameters, trends in statistical moments, etc. Most studies of FFA considering nonstationarity adopted time-varying models to describe the variation of flood regimes [17,18,19,20]. However, applying a time-dependent model is inadequate because it fails to account for the change caused by climate change and/or human activities, which are usually the fundamental drivers. Thus, recent studies [21,22,23,24,25,26] have investigated the possibility of incorporating external forcings, such as climate factors, into models for FFA, and their work reported a better description of variation in flood regimes. However, it is unclear how these methods can be used to project future conditions since they are not physically based.

An alternative to statistical approaches is to derive flood frequency distribution using hydrological modeling, where design flood value for a specific return period was estimated using simulated runoff series from a hydrological model forced by observed or synthetic rainfall series [27,28,29,30,31,32,33,34]. The simulation can be event-based or continuous, with the latter approach increasing in use for its independence on the assumption about equal return periods for the design rainfall event and resulting design flood. Moreover, it becomes less necessary to account for the initial soil moisture conditions within the catchment when applying continuous simulation. A detailed review of applications of continuous simulation in design flood estimation can be found in the work of Boughton and Droop [35]. Through the continuous simulation approach, the design flood value can be determined by derived flood frequency analysis using long-term simulated flows series [28,29,30,31,32,33]. To this end, long-term continuous rainfall series with a sufficient temporal resolution are required, especially for small catchments. Due to the limitation of the resolution and length of available observed data, some researchers turn to synthetic rainfall series generated by stochastic rainfall models [36,37,38,39,40]. Moreover, recent years have witnessed the increasing availability of satellite-based rainfall products (e.g., GPM) and reanalysis rainfall products (e.g., ERA5 and UERRA), which can serve as a good solution for data-scare areas [41,42].

The potential advantage of deriving flood frequency analysis using a hydrological model over a statistical one is being able to account for the flood-related physical processes explicitly. However, this method may be invalid under changing environments, such as alteration in land use or the construction of new flood protection projects, whose effect cannot be reflected by the observed historical flood records. Thus, some recent studies try to diagnose and interpret hydrological nonstationarity, aiming to improve the structure of hydrological models and corresponding simulation/forecasting performance under such conditions [43,44,45]. The current state of the art is the application of conceptual hydrological models with time-varying parameters, which are usually assumed to represent changing watershed conditions. For example, Wallner and Haberlandt [44] reported a significantly improved streamflow simulation performance of the HBV model when the model parameter was made dependent on the climate indices. Deng et al. [45] linked the parameters of the two-parameter monthly water balance model with catchment properties. Their case study in southern China shows that the incorporation of time-varying parameters can lead to a more reliable prediction of monthly streamflow under changing environments.

However, most recent studies focus on deriving flood frequency analysis using a statistical model with time-varying parameters. Although the method of deriving flood frequency using continuous hydrological simulation and hydrological modeling with time-varying parameters have been validated in previous studies, few studies link them together to derive flood frequency analysis under nonstationary conditions. Thus, this study aims to derive flood frequency using a hydrological model with time-varying parameters under changing watershed conditions. The flood quantiles in different sub-period were derived using continuous simulation results by the stationary and nonstationary GR4J model, with the Weihe River basin (WRB) as the study area. Then the differences in estimated flood quantiles and magnitudes by the stationary and nonstationary models were presented to emphasize the importance of applying nonstationary modeling in flood designing.

2. Study Area and Data

The Weihe River basin (WRB) in northwestern China was chosen as the study area for its high nonstationarity, with runoff series from Huaxian station (near the outlet of the basin) on a declining trend over the past decades [46,47,48,49,50,51,52,53,54,55], mainly due to the intensifying human activities associated with water extraction and diversion [49,50]. The WRB is located between 104°00′ E–110°20′ E and 33°50′ N–37°18′ N, with a length of 818 km and a drainage area of 134,800 km², as shown in Figure 1. The elevation within WRB ranges from 340 to 3671 m, with a decrease from northwest to east.

The daily runoff series from 1951 to 2011 were obtained from the Huaxian station. The daily rainfall and mean temperature series were obtained from the 21 meteorological stations (see Figure 1) within or around WRB for the same period as the runoff data. The daily potential evaporation of the corresponding period was estimated by the modified Hargreaves equation [56]. To account for the possible impact of human activities on the underlying conditions and flow regimes within WRB, the socio-economic factors, including population (POP), gross domestic product (GDP), irrigated land area (IA), cultivated land area (CA), reservoir storage (RS) and soil-water conservation land area (SWC) were considered as the candidates of external covariates for time-varying parameters. These factors can serve as indicators for the alterations of watershed conditions (e.g., land use and land cover) and have been used to assess the impact of human activities on flow regimes in previous work [57,58]. Considering the majority of cities, population, and farmed land area within WRB located in the Shaanxi Province, the covariates mentioned above for the period of 1951–2011 were represented by data of this province [59,60], as shown in Figure 2.

3. Methodology

3.1. GR4J Model: Stationary and Nonstationary

The hydrological model applied in this study is the lumped conceptual GR4J model [61]. The GR4J model was chosen for its parsimony in model parameters while still maintaining adequate process representation, including fast and slow responses together with a soil moisture component. The GR4J model has 4 parameters that are required to be calibrated, as presented in

Table 1. Description of the parameters within the GR4J model.

Parameter	Description	Unit	Feasible Range
${\theta }_1$	production storage capacity	mm	20–1200
${\theta }_2$	groundwater exchange coefficient	mm	−5–3
${\theta }_3$	one day ahead maximum capacity of the routing store	mm	20–500
${\theta }_4$	time base of unit hydrograph	days	1–5

.

The stationary and two nonstationary versions of the GR4J model were applied for continuous simulation purposes in the study, i.e., the stationary model (model zero), where parameter values keep constant over time; the time-dependent model (model one), where the parameters vary as a function of time; and the time-varying model (model two), where the parameters vary as a function of covariate related to watershed underlying conditions. For both models one and two, the parameter ${\theta }_1$ was allowed to vary over time because it stands for the production storage capacity of the watershed and has proved to be the most sensitive parameter of GR4J in previous studies [62,63,64].

For model one, the parameter ${\theta }_1$ was treated as a liner function of time, as follows

$${\theta }_1^{1,t}=a+bt$$

(1)

where ${\theta }_1^{1,t}$ stands for ${\theta }_1$ in model one for the year $t$. The coefficients {$a$, $b$} are unknown and are therefore calibrated together with the parameter ${\theta }_2$, ${\theta }_3$ and ${\theta }_4$ in model one.

To reveal how the parameter ${\theta }_1$ would vary with the underlying conditions variation caused by human activities, ${\theta }_1$ in model two was treated as a function of multiple covariates, as follows

$${\theta }_{1,t}^2={\beta }_0+{\displaystyle \sum _{i=1}^m{\beta }_i{h}_i(x_{i,t})}$$

(2)

where ${\theta }_{1,t}^2$ represents the value of ${\theta }_1$ in model two for the year $t$, $x_{i,t}(i=1,2,...,m)$ are $m$ explanatory variables for the year $t$, ${\beta }_i(i=0,1,...,m)$ are $m+1$ regression coefficients. Note that function $h_i(\cdot )$ stands for the dependence of the parameter ${\theta }_{1,t}$ on explanatory variables $x_{i,t}$, which can be identical, logarithmic, or exponential, so as to represent possible linear or nonlinear relationships of the parameter ${\theta }_1$ to changing watershed conditions.

3.2. Parameter Estimation and Model Selection

For model two, the selection of proper covariates for the parameter ${\theta }_1$ from several candidate explanatory variables is essential. In total, 6 covariates (POP, GDP, IA, CA, RS, and SWC) were considered in this study, as mentioned above. In order to investigate the potential for the nonstationary GR4J model in runoff simulation, all combinations of the 6 external covariates (in total 63 covariate combinations, $C_6^1+C_6^2+C_6^3+C_6^4+C_6^5+C_6^6=63$) were tested, and their corresponding performance was evaluated in terms of the Nash-Sutcliffe efficiency coefficient (NSE) metric. The combination of covariates corresponding to the highest NSE metric was regarded as optimal.

Parameters of all 3 models were calibrated using the SCEM [65] method, which has been widely used for practical assessment of parameter uncertainty in hydrological modeling [65,66,67]. The most likely parameter set in terms of NSE metric can be identified through the SCEM method, as well as their underlying posterior distribution [65,66,67,68,69]. To account for uncertainty in parameters, 1000 parameter sets were then sampled from the posterior probability distribution. Of the 1000 sampled parameter sets, the one corresponding to the highest NSE metric was regarded as the approximation of the global optimum and used for model evaluation and selection. Considering the limitation of model comparison and selection based on a single metric, additional criteria were adopted, including the Kling-Gupta efficiency coefficient (KGE), the relative error (RE) between simulated and observed annual total volume, and the daily-scale flow-duration curves, which enables comparison between the probability distribution of simulated and observed runoff and may contribute to figuring out the potential biases (e.g., possible underestimation of high flows).

The period of 1951–1990 was used for parameter calibration (for all models) and model selection (for model two only). The period of 1991–2011 was used for model validation. In order to reduce the possible impact of initial conditions, a 1-year warm-up period was applied for both the calibration and validation period, and the simulation within the warm-up period was excluded before the model evaluation procedure.

Note that WRB is highly seasonal, typically with low flow during winter. In order to avoid the possible impact of low flow on the analysis, a peak-over-threshold calibration strategy was adopted. This strategy excludes flows below a threshold of 100 m³/s for the model evaluation procedure, thus placing greater emphasis on the simulation performance of high flow [70]. This censoring threshold was chosen because roughly 40% days from 1951 to 2011 witnessed a flow below 100 m³/s for WRB.

3.3. Deriving Flood Frequency Distribution through Monte Carlo Simulation

For flood frequency analysis under the stationarity assumption, the underlying watershed conditions were assumed to be time-invariant; thus, a constant value of design flood for a specific return period could be determined. However, under nonstationary conditions, the design flood value for a certain return period may change over time due to the evolving watershed conditions. In order to achieve flood frequency analysis under nonstationary conditions, the historical period was divided into several sub-period. Based on the parameter from each sub-period, the Monte Carlo simulation method using the GR4J model was adopted to generate long-term runoff series; subsequently, the flood frequency distribution for certain sub-period could be determined.

As for model input for the Monte Carlo simulations, long-term rainfall series were generated by the weather generator BCC/RCG-WG [71], which was developed by Beijing Climate Center of China Meteorological Administration and Regional Climate Group at the University of Gothenburg. The stochastic model generates continuous rainfall series whose arrival time, duration, and intensity are all random, governed by a gamma distribution with parameters assumed to change monthly [72]. The calibration and validation of monthly parameters of the stochastic model have been conducted for 672 meteorological stations all over China in previous works [72,73]. Details on the model and its calibration/validation procedure can be found in the work of Liao et al. [73].

In practice, maximum annual flows are mainly relevant to annual maximum rainfall events and their antecedent rainfalls. In the case of WRB, Xiong, and Jiang [25] reveal a considerable correlation between annual maximum daily runoff and its 15-day antecedent rainfall, which shows a negligible variation. This study assumed that the characteristics of storms (annual maximum rainfall events) were stationary over the study period, thus enabling investigations on how the variation of watershed conditions affects the design flood value. Note that when applied in future scenarios, a perturbation with a climate change signal to the weather generator is necessary to generate rainfall series corresponding to the possible climate scenarios.

The GR4J model and the weather generator were used in the Monte Carlo simulation to generate daily runoff series of 100 years. Ten rainfall realizations generated by the stochastic model were used to consider uncertainty in rainfall (i.e., input for the GR4J model). As mentioned above, a total of 1000 parameter combinations were obtained from the posterior parameter distribution to account for uncertainty in the parameters of GR4J. Thus, in total 10,000 model runs were applied to derive the flood frequency distributions and corresponding uncertainty for certain sub-period.

4. Results

4.1. Simulation Performance of the Stationary and Nonstationary GR4J Models

For model two, 63 covariate combinations were tested, and their corresponding performance was evaluated in terms of the NSE metric. For the sake of brevity, only 4 cases are presented in

Table 2. Model evaluation and selection for the nonstationary GR4J model (model 2).

Case	Covariate		Equation for the Time-Varying Parameter ${\mathit{\theta }}_{1,\mathit{t}}$	NSE [-]	KGE [-]	RE [%]
	RS	SWC
C0			${\theta }_{1,t}=542.1$	0.711	0.801	11.5
C1	✓		${\theta }_{1,t}=290.3+195.6\frac{R{S}_t}{\overline{RS}}$	0.753	0.845	8.9
C2		✓	${\theta }_{1,t}=296.8+152.4\frac{SW{C}_t}{\overline{SWC}}$	0.748	0.839	9.4
C3	✓	✓	${\theta }_{1,t}=304.5+102.1\frac{R{S}_t}{\overline{RS}}+69.2\frac{SW{C}_t}{\overline{SWC}}$	0.765	0.853	8.5

for illustration purposes. Note that for each covariate combination, 1000 parameter sets were obtained from the posterior probability distribution, and only the highest-NSE one was presented here for model selection purposes. Case C0 was designed to represent the stationary GR4J model. Under cases C1–C3, the external covariate RS and SWC were introduced in different schemes to describe the temporal variation of the production storage capacity parameter ${\theta }_1$. It can be found that the runoff simulation performance improved when the parameter ${\theta }_1$ was allowed to change over time. For example, the NSE value is 0.753 under case C1, where only the covariate RS was considered when compared to 0.711 for the stationary GR4J model. The highest NSE value (0.765) was achieved under case C3, where both RS and SWC were incorporated into the time-varying parameter. Thus, the time-varying model under C3 was chosen as the optimal model two. Table 2 also gives the equation for the time-varying parameter ${\theta }_{1,t}$ under 4 cases, where $\overline{RS}$ and $\overline{SWC}$ represent the mean value of covariate RS and SWC for the period from 1951 to 2011.

Table 3. Simulation performances of the stationary and nonstationary GR4J model for WRB.

Model	Covariates	Calibration Period (1951–1990)			Validation Period (1991–2011)
		NSE [-]	KGE [-]	RE [%]	NSE [-]	KGE [-]	RE [%]
Model zero	none	0.711	0.801	11.5	0.692	0.791	13.5
Model one	time	0.743	0.832	9.8	0.714	0.812	10.1
Model two	RS, SWC	0.765	0.853	8.5	0.716	0.818	9.8

gives a summary of statistics of simulation performance using the stationary GR4J model (model zero) and the modified versions (models one and two). Note that there are numerous versions of model two, and only the highest-NSE one is present here. The NSE, KGE, and RE metrics in Table 3 indicate that all three models achieved acceptable simulation performance for both calibration and validation periods. Model two, which incorporates RS and SWC as external variables, achieved the highest NSE value for both the calibration and validation period, followed by model one and model zero. This indicates that the nonstationary GR4J model can achieve a better simulation performance than the stationary one.

Figure 3 shows the comparison between the observed and simulated flow-duration curves (FDC) from models zero, one, and two with the parameter set corresponding to the highest NSE metric. All models achieved a good simulation performance for flow with an exceedance probability of less than 20%. Model two apparently outperforms model zero and one in terms of high flows. The flow-duration curves were also presented for the autumn season when floods are most likely to occur (Figure 3b). The results indicate good performance for models zero, one, and two in simulating high flows and poor performance in reproducing low flows. This is expected given NSE is the objective function and the peak-over-threshold calibration strategy was applied.

Predictions of the annual maximum flow were obtained from the simulated runoff series using the GR4J models (model zero, one, and two) with the 1000 parameter sets. Figure 4 compares the observed and simulated annual maximum flows for the period of 1951–2011. The observed annual maximum flows are represented by blue dots. The 95% prediction uncertainty bounds associated with the parameter are represented by the dark-shaded area, and the 95% prediction interval related to the model residuals corresponds to the light-shaded area. The NSE metric for the annual maximum flood series is denoted by NSE_AMF in Figure 4.

The results from the stationary model with a constant parameter set (model zero) and the nonstationary model assuming temporal dependence only (model one) can adequately describe the variation in maximum annual flows (see Figure 4a,b), with most observed floods lying within the uncertainty bounds. For WRB, model two produced the most reliable estimation of maximum annual flows under nonstationary conditions (NSE_AMF = 0.808) when compared to model zero (NSE_AMF = 0.721) and model one (NSE_AMF = 0.748), together with comparative narrower uncertainty bounds. This also proves the benefit and advantage of incorporating reservoir storage RS and soil-water conservation project SWC that are related to local regulation strategies as covariates.

Note that the incorporation of the time-varying parameter can be regarded as a modification of the model structure since one or more parameters require calibration are incorporated. Thus, whether a significant improvement in simulation performance can be achieved highly depends on the rationality of the modified model structure. For models with time-varying parameters, the selection of appropriate covariate to describe the dynamic of watershed conditions is of great importance. In this study, the combination of covariate RS and SWC leads to the highest NSE metric. Other covariates related to watershed conditions (e.g., vegetation cover) may achieve a more significant improvement in simulation performance. This is an interesting topic that deserves further research.

4.2. Dynamic Variation of Time-Varying Parameter

Figure 5 presents the dynamic variation of the parameter ${\theta }_1$ for models one and two over the whole period with a 5-year spacing. The production storage capacity determines the catchment responsiveness and the allocation of rainfall input into interception, evapotranspiration, groundwater, and runoff. A larger production storage capacity means stronger resistance effects of rainfall input on runoff generation, i.e., weaker responsiveness of the watershed to rainfall variability.

For both models one and two, the increasing trend in the parameter ${\theta }_1$ is apparent. This can be partially explained by the aforementioned external variables RS (reservoir storage) and SWC (soil-water conservation land area). The construction of large reservoirs in WRB sprang up in the 1960s and reached an upsurge in the late 1970s, which may lead to a significant shift in flood regime and, therefore, a noticeable change in flood magnitude at the outlet of the basin. The soil-water conservation projects (including the construction of terraces, afforestation, etc.) that started in the late 1950s and still persist today also contributed to larger storage for the rainfall during runoff generation, thus causing a much lower flow given the same rainfall input.

4.3. Performance of the Stochastic Rainfall Model

In order to validate the reliability of the stochastic rainfall model mentioned above, 10 realizations of long-term synthetic rainfall series on a daily basis (each 100-year in length) were generated for 21 meteorological stations within WRB. For brevity, only the result of the meteorological station Tianshui is shown here.

Table 4. Event characteristics of the observed rainfall series (OBS) and the synthetic rainfall series (STOCH) for the meteorological station Tianshui.

Event Characteristics	Unit	OBS	STOCH
Number of rainfall events per year	-	54	48
Mean of event volume	mm	9.54	10.08
Standard deviation of event volume	mm	14.76	13.58
Skewness of event volume	-	3.24	3.52

shows comparisons of the event characteristics of the 61-year observed rainfall series and the synthetic rainfall series generated by the stochastic rainfall model. The statistics of the generated rainfall series show good agreement with those of the observed historical records, with a slight overestimation of the mean volume of a rainfall event. The standard deviation and skewness of event volume were also well reproduced.

A frequency analysis was conducted on the annual maximum rainfall series for both the observed and synthetic rainfall series. Note that rainfall series can be generated for long periods using the stochastic model but only can be compared with the shorter observation statistic. Selected results are presented in Figure 6, where the observed values mainly lie within the boundary of the corresponding synthetic rainfall records. This indicates a good reproduction of rainfall series using the stochastic model, and the synthetic rainfall series is adequate for the purpose of Monte Carlo simulation.

4.4. Flood Frequency Derivation Using the Stationary and Nonstationary Model

Figure 7 and Figure 8 show the comparison of flood quantiles for an exceedance probability of 0.01 using the stationary model (model zero) and nonstationary models (models one and two). The 95% uncertainty bounds were estimated from simulation results from the 10,000 model runs, as mentioned above.

For the stationary model (model zero), the flood quantile for an exceedance probability of 0.01 (i.e., return period of 100-year) was a constant value of 6521 m³/s (see Figure 7). For model one, the corresponding flood magnitudes were continuously on a downward trend. This is expected since the parameter ${\theta }_1$ of model one was time-dependent only. An FFA of the Huanxian station using model two (see Figure 8) shows that the flood quantiles for an annual exceedance probability of 0.01 range from 4187 m³/s (in 2011) to 8335 m³/s (in 1960). A noticeable decline occurred in the late 1970s when the reservoir storage and soil-water conservation area were undergoing drastic changes, indicating a more reasonable interpretation of the variation of flood magnitudes over time compared to model one. The variation in flood quantiles is in line with the work of Xiong et al. [25] and Hesarkazzazi et al. [74].

The results also stress the problems when assuming stationarity in flood frequency analysis and design. The application of the nonstationary models indicates that flood quantiles and magnitudes for WRB underwent significant variability over the past decades. In such cases, the assumption of stationarity may lead to an apparent misestimate of design flood for the future.

5. Discussion

Applying time-varying parameters in hydrological modeling has proved to be a useful approach to diagnosing hydrological nonstationarity under changing underlying watershed conditions. The key issue remains the selection of proper external covariates that are related to watershed conditions variation. This study chooses six external covariates and searches for the optimal one in a brute-force approach. A more comprehensive approach to identifying the most likely covariate combination requires further investigation.

The introduction of external covariates helps the models grasp the variations in flood frequency during the calibration and validation periods. This is in line with the work of Xiong et al. [25] and Jiang et al. [55]. However, before utilizing such a model as a predictive tool for future scenarios, the predictions of climate and anthropogenic indices are necessary.

6. Conclusions

In this study, flood frequency analysis considering nonstationarity by continuous simulation using the GR4J model with time-varying parameters was conducted for the Weihe River basin, and the main findings are as follows:

• Improved model performance can be achieved when the parameter ${\theta }_1$ (representing the production storage capacity) is treated as time-dependent or as a function of external variables that reflect changes in the hydrological responses within WRB;
• Incorporating the watershed conditions as covariates for the model parameter (model two) can better describe nonstationarity in the flood frequency and magnitude in WRB than the over-simplified time-dependent model (model one);
• The nonstationary model can achieve a more rational description of variations in the frequency and magnitude of floods over time, which reveals the deficiency when applying a stationary model under changing watershed conditions.