Assessment of the effectiveness of a multi-site stochastic weather generator on hydrological modelling in the Red Deer River watershed, Canada

ABSTRACT

To improve the convergence of multiple-site weather generators (SWGs) based on the brute force algorithm (MBFA), a genetic algorithm (GA) is proposed to search the overall optimal correlation matrix. Precipitation series from weather generators are used as input to the hydrological model, the soil and water assessment tool (SWAT), to generate runoff over the Red Deer watershed, Canada for further runoff analysis. The results indicate that the SWAT model using SWG-generated data accurately represents the mean monthly streamflow for most of the months. The multi-site generators were capable of better representing the monthly streamflow variability, which was notably underestimated by the single-site version. In terms of extreme flows, the proposed method reproduced the observed extreme flow with smaller bias than MBFA, while the single-site generator significantly underestimated the annual maximum flows due to its poor capability in addressing partial precipitation correlations.

KEYWORDS:

• multi-site stochastic weather generator
• genetic algorithm
• hydrological model
• spatial correlation
• Red Deer River watershed

Editor A. Castellarin Associate editor A. Requena

1 Introduction

Stochastic weather generators (SWGs) are numerical tools that aim to synthesize long time series of various climate variables, such as precipitation and temperature, which are statistically identical to historical observations (Wilks and Wilby 1999). The Weather Generator (WGEN), developed by Richardson and Wright (1984), and the Long Ashton Research Station-Weather Generator (LARS-WG) (Semenov and Barrow 1997) are two commonly-used single-site SWGs. The former is a parametric-type model, which relies on a first-order two-state Markov chain to generate the occurrence of wet and dry days and a two-parameter gamma distribution to simulate the precipitation amounts on wet days. The latter is a semi-parametric type model, relying on semi-empirical distributions to model precipitation occurrence and amounts. The output from single-site SWGs is spatially independent climate variables, making it difficult to match the actual climate data with strong spatial correlations. In the past, various types of multi-site SWGs have been proposed (Baigorria and Jones 2010), including parametric (Wilks 1998, Qian et al. 2002, Brissette et al. 2007, Chen et al. 2014, Evin et al. 2018), non-parametric (Burton et al. 2008, Leander and Buishand 2009) and hybrid (Palutikof et al. 2002, Breinl et al. 2015, 2017) models.

In view of parametric-type modelling, the general idea of developing a multi-site SWG is to extend the existing single-site version (e.g. WGEN) by driving it with temporally independent but spatially correlated random numbers. Based on the empirically monotonic relationship between the correlations of the random numbers and the observed correlations of precipitation occurrence and amounts, Wilks (1998) constructed the corresponding empirically-derived curves to identify the needed correlation of random numbers. Wilks (1999) and Qian et al. (2002) enhanced the framework of the Wilks (1998) approach by adopting a high-dimensional autoregression technology to produce the non-precipitation variables, including maximum and minimum temperatures as well as daily solar radiation. Even though this approach can take into account the spatial dependence of climatic variables, it has two main drawbacks: (i) the huge computational burden due to its low-efficient Monte Carlo sampling process in construction of empirically-derived curves for all possible pairs of stations, and (ii) the ill-defined correlation matrix problem in estimating correlation on a pairwise basis (Mehrotra et al. 2006, Khalili et al. 2007, Chen et al. 2014). Brissette et al. (2007) proposed a brute force algorithm (BFA) to improve the computational efficiency, and used a spectral decomposition method (Rebonato and Jäckel 1999) to solve the problem of ill-defined correlation matrices in the Wilks (1998) approach. The BFA is a recursion formula, capable of continually updating the needed random number correlation matrix with a varying step size and search direction according to the residual errors of the correlation between observed and generated precipitation occurrences or amounts. Chen et al. (2014) coded the methodology of Brissette et al. (2007) based on MATLAB, and published a MATLAB software package (MulGETS), adding functions of simulating maximum and minimum temperatures. From the perspective of optimization algorithms, the BFA may easily encounter convergence problems as it is generally challenging to obtain an overall minimum solution. Meanwhile, the iteration may be more frequently stopped by the predefined maximum iteration number, instead of the error tolerance, which may increase the risk of obtaining a poor result (Brissette et al. 2007, Chen et al. 2014). In fact, the problem of finding temporally independent but spatially correlated random numbers belongs to an unconstrained nonlinear optimization problem. Many heuristic solvers, such as genetic algorithm (GA), have been reported to be advantageous in dealing with such a problem in the hydrological field (Muleta and Nicklow 2005, Arsenault et al. 2014, Chlumecký et al. 2017, Manfreda et al. 2018). Thus, it is desirable to use a similar approach to replace the BFA in the framework of Brissette et al. (2007) for achieving global solutions.

From a practical application point of view, SWGs are widely used to generate long climate series, which are used as the input to hydrological models for runoff analysis. Theoretically, the climate data generated by single-site SWGs would lose the information of spatial dependence, leading to potential modelling errors, especially for simulating extreme hydrological events (Caron et al. 2008). Khalili et al. (2011) used a spatial autocorrelation-based multi-site SWG and a single-site WGEN, respectively, to drive a physically based distributed hydrological model (HydroTel; Fortin et al. 2001) for streamflow simulation in the Chute du Diable River Watershed (approx. 9700 km²), Quebec, Canada. The study suggested that it may be preferable to use multi-site SWGs, coupled with distributed hydrological models, for hydrological simulation in basins where precipitation varies spatially in a significant manner. More recently, Li et al. (2017) evaluated the multi-site SWG (i.e. MulGETS) and the single-site SWG (i.e. WGEN) in hydrological modelling in the Jing River watershed (approx. 45 421 km²), China. Alodah and Seidou (2019) examined the effects of a single-site SWG (WeaGETS inspired from the WGEN approach), the multi-site MulGETS and a non-parametric k-nearest-neighbour weather generator, on the hydrological output of the soil and water assessment tool (SWAT) for the South Nation watershed in the Eastern Ontario, Canada (approx. 4000 km²). The above-mentioned three studies concluded that distributed hydrological models driven by multi-site SWGs performed better in simulating extreme streamflow, which was consistent with the findings from Khalili et al. (2011). Even when using a lumped hydrological model, a multi-site SWG can be advantageous in helping evaluate the impact of rainfall–runoff modelling on hydrological output, as reported by Breinl (2016), in two Alpine study areas, the Salzach basin on the border of Germany and Austria (4637 km²) and the Ubaye basin in France (548 km²).

However, there are also some studies that have come up with different results in evaluating the effect of multi-site SWGs on hydrological modelling. Watson et al. (2005) assessed the response of a distributed hydrological model (SWAT) to the weather data generated by a single-site SWG (Srikanthan 2005) and a multi-site one (Wilks 1998) over the Woady Yaloak River catchment (306 km²) in Australia. Little difference was found between the two models, probably due to the small size of the watershed and the flat topography, which limited the spatial variability of weather variables. Chen et al. (2016) compared two methods for hydrological modelling over the large Lac-Saint-Jean watershed (45 432 km²), Québec, Canada. One method was based on the combination of a single-site SWG (i.e. WGEN) using a lumped approach and a conceptually lumped hydrological model (HSAMI); the other was based on the combination of a multi-site SWG (i.e. MulGETS) and a physically-based distributed hydrological model (CEQUEAU; Ayadi and Bargaoui 1998). The two methods both performed well in simulating the mean and standard deviation of monthly average and extreme flows, but the advantage of the distributed hydrological model using the multi-site SWG was not obvious, possibly due to the strongly snowmelt-dominated hydrographs in the Lac-Saint-Jean watershed. These studies demonstrated that the performance of multi-site SWGs may be affected by many factors, such as watershed size and topography, climate variability and hydrological type of watershed. More quantitative assessment of the impact of using multi-site SWGs on hydrological modelling is still highly demanded.

Therefore, the objective of this research is twofold. Firstly, GA is used to search a global optimal correlation matrix for generating the spatially correlated multi-site precipitation series by replacing the BFA of the framework of Brissette et al. (2007). Its searching performance related to different error functions and effectiveness in simulating the spatial correlation of precipitation occurrence and amounts is carefully evaluated. Secondly, the flow predictions from a hydrological model (i.e. SWAT) driven by both single and multiple weather generators are evaluated for the Red Deer River watershed, Alberta, Canada. The statistics of precipitation from single and multiple weather generators are analysed to give information on their possible effects on hydrological modelling. The hydrological responses of SWAT to single and multiple weather generators will be compared to give suggestions for the choice of a weather generator in this basin.

2 Data and methodology

2.1 Study area

The Red Deer River watershed (50º29′–52º59′N, 110º01′–116º06′W) is located in southern Alberta, Canada (Fig. 1) and covers approximately 47 954 km² (Mishra and Coulibaly 2010). The Red Deer River originates from the eastern slopes of the Rocky Mountains in the Banff National Park near Lake Louise, with a total length of 724 km (Tanzeeba and Gan 2012). The region that this river traverses is mountainous in the west, prairie in the east and cropland in the middle. With a continental climate, the region has an average annual temperature of 4°C and a median annual precipitation of approximately 393 mm (Kerr and Cooke 2017). We used the precipitation data from four climate stations along the mainstream of the Red Deer River. Table 1 lists the background information of the four climate stations. The observed period of all the climate stations covers 24 years from 1979 to 2005 (excluding 1981, 1982 and 1990 due to lack of data). Further information regarding the average monthly precipitation and inter-station statistics can be found in the Supplementary material (Fig. S1). The majority of annual precipitation (about 70%) occurs from May to September and the annual maximum value of the inter-station correlation coefficient is 0.54. According to the Red Deer River State of the Watershed Report (https://www.rdrwa.ca), the multi-year average water flow rate of the Red Deer River watershed is about 70 m³/s at the watershed outlet. There have been several high streamflow and flood watch/warning advisories in the upstream of the Red Deer River since 2005, which resulted from snowmelt early in the year and high precipitation events throughout the summer months.

Table 1. Information on the climate stations in this study.

Name	Climate ID	Longitude	Latitude	Elevation	Period
Oyen Cappon	CA003024961	–110.52	51.17	792.5	1974–2013
Brooks North^a	CA003030862	–111.8833	50.6167	759	1960–2008
Drumheller Andrew^b	CA003022136	–112.8667	51.4667	719	1960–2008
Sundre Garrington^c	CA003026259	–114.45	51.9667	991	1978–2008

^aMissing data was supplemented by averaging the data from the neighbouring stations Brooks 1 (CA003030840) and Brooks (CA003030QLP);

^bMissing data was supplemented by averaging the data from the neighbouring stations Drumheller (CA003022120), Drumheller Institution (CA003022137), Drumheller City (CA003022139) and Drumheller East (CA0030221LG);

^cMissing data was supplemented by averaging the data from the neighbouring stations Sundre (CA003026258) and Sundre Ranger (CA003026260).

Figure 1. Red Deer River watershed and the used hydro-meteorological stations.

2.2 Hydrological modelling

The well-known SWAT model (Tripathi et al. 2004) is used in this study to simulate hydrological processes. The purpose, characteristics and sources of the primary data used for the SWAT modelling are presented in the Supplementary material (Table S1). The related digital elevation model (DEM), land-use and soil data were derived from the NASA (US National Aeronautics and Space Administration) Shuttle Radar Topographic Mission (SRTM) project (Butt and Bilal 2011), the US Geological Survey Global Land Cover Characterization (GLCC) database (Brown et al. 1999), and the Food and Agriculture Organization of the United Nations (FAO) (Nachtergaele et al. 2009), respectively. The precipitation data at the four stations (Table 1) and other meteorological data (i.e. air temperature, relative humidity, wind speed and solar radiation) were collected from the Government of Canada (http://climate.weather.gc.ca) and the National Centers for Environmental Prediction Climate Forecast System Reanalysis (CFSR) (Dile and Srinivasan 2014). The daily streamflow data from 1979 to 2005 at the station 05CK004 (shown in Fig. 1) was derived from Environment Canada (http://www.ec.gc.ca/rhc-wsc).

Based on the DEM and defined outlets, the Red Deer River watershed was divided into 33 sub-basins (Fig. 1), with areas ranging from 6 to 4535 km² (average: 1398 km²). Each sub-basin was further divided into multiple hydrologic response units (HRU) based on the information of slope, land use and soil type, resulting in 97 HRUs. The SWAT model was calibrated by using the sequential uncertainty fitting algorithm version-2 (SUFI-2) tool (Mousavi et al. 2012, Abbaspour et al. 2015) based on the observed daily streamflow data from the hydrometric station (05CK004). Table 2 presents the description of the main hydrological parameters for calibration and their calibrated values. The Nash-Sutcliffe efficiency ($E_{\mathrm{N}\mathrm{S}}$), is used as the index for evaluating model performance (Confesor and Whittaker 2007):

(1) $E_{\mathrm{N}\mathrm{S}}=1-\frac{{\sum }_i{\left(O_i-S_i\right)}^2}{{\sum }_i{\left(O_i-\bar{O}\right)}^2}$(1)

Table 2. SWAT parameters and the best calibrated values.

			Calibration range^c
Hydrological process	Parameter name^a	Calibration method^b	Min.	Max.	Best calibrated value
Snow	SFTMP	v	–1.5	1	–0.93
	SMTMP	v	0	3	1.30
	SMFMX	v	1.4	6.9	5.65
	SMFMN	v	1.4	6.9	5.12
	TIMP	v	0	1	0.14
Groundwater	ALPHA_BF	v	0.001	1	0.92
	REVAPMN	a	1	10	7.82
	GW_DELAY	v	30	450	168.07
	GW_REVAP	a	0	0.2	0.13
	GWQMN	a	0	2	0.89
	RCHRG_DP	v	0.001	1	0.96
Soil	SOL_AWC	r	–0.2	0.2	–0.01
	SOL_K	r	–0.8	0.8	0.08
Evaporation	EPCO	v	0.001	1	0.22
	ESCO	v	0.001	1	0.93
Interflow	OV_N	r	–0.2	0	−0.18
	SLSUBBSN	r	0	0.2	0.01
Runoff	CN2	r	–0.2	0.2	0.03
Channel	CH_N2	v	0.001	0.3	0.28
	CH_K2	v	50	130	107.42

^aFor definition of these hydrological parameters see the Supplementary material, Table S2 (Faramarzi et al. 2015);

^bThe calibration methods, a, v and r, denote an absolute increase, a replacement and a relative change to the default parameter value, respectively;

^cThe ranges of snow-related parameters is based on Wang and Melesse (2005) and the ranges of the remaining parameters on Rouhani et al. (2007) and Her et al. (2015).

where $O_i$ and $S_i$ are, respectively, the observed and simulated streamflow on the ith day and $\bar{O}$ is the average value of $O_i$.

2.3 Stochastic weather generator

2.3.1 Multi-site precipitation generation based on brute force algorithm (MBFA)

The WGEN uses a first-order two-state Markov chain to determine the daily precipitation occurrence and a two-parameter gamma distribution to generate the daily precipitation amount after a wet day is confirmed (Richardson 1981, Richardson and Wright 1984, Wilks 1998). It is meant for generating precipitation series at a single location and is incapable of addressing spatial correlation for multiple sites. The BFA (Brissette et al. 2007) is used to search for a suitable correlation matrix to produce spatially correlated random numbers and then drive the precipitation generator. The procedure (Fig. 2) can be briefly described as follows.

Figure 2. Methodology framework of the overall processes in this study. Note: the dotted, dashed, dash-dot and solid arrow lines denote the operation procedures relevant to WGEN, MBFA, MGA and all, respectively.

Firstly, we assume $C_k$ is a candidate correlation matrix with size of $n\times n$ (where $n$ is the number of stations) at the $k\mathrm{t}\mathrm{h}$ iteration $(k=0,1,2,\dots )$, and initialized with the observed correlation matrix of precipitation occurrence or amounts ($C_{\mathrm{o}\mathrm{b}\mathrm{s}}$), i.e. $C_0=C_{obs}$. If $C_k$ is a non-positive-definite matrix, it needs to be modified as a positive-definite one by using the spectral decomposition method (Rebonato and Jäckel 1999). Next, we can produce the spatially correlated uniform distributed matrix $R_c^U$ by using the Cholesky factorization method (Chen et al. 2014) for the positive-definite matrix $C_k$. Then, $R_c^U$ is imported into the WGEN to generate spatially correlated precipitation occurrences or amounts. Meanwhile, we can calculate the generated correlation matrix of precipitation occurrence or amount ($C_{\mathrm{g}\mathrm{e}\mathrm{n}}$). In order to evaluate the overall differences between the observed and generated correlation matrixes, an error function ($z$) can be defined as (Brissette et al. 2007):

(2) $z=\left\{f\left(c_{i,j}^{\mathrm{\Delta }}\right)|c_{ij}^{\mathrm{\Delta }}\in \mathrm{\Delta }C,\mathrm{\Delta }C=C_{\mathrm{o}\mathrm{b}\mathrm{s}}-C_{\mathrm{g}\mathrm{e}\mathrm{n}}\right\}$(2)

where $i,j=1,2,\dots n$,; and the definition of error function, $f\left(c_{i,j}^{\mathrm{\Delta }}\right)$, is defined in Section 2.4. If $z$ is lower than a predefined error tolerance (i.e. $\rho$) or the number of iteration ($k$) is more than the predefined maximum iteration number (i.e. $\tau$), the iteration process should stop. Otherwise, $C_k$ is replaced by $C_{k+1}$ and the iteration process will restart. The relationship between $C_{k+1}$ and $C_k$ is (Brissette et al. 2007):

(3) $C_{k+1}=C_k+\eta \cdot \left(C_{\mathrm{o}\mathrm{b}\mathrm{s}}-C_{\mathrm{s}\mathrm{y}\mathrm{n}}\right)$(3)

where $\eta$ is a convergence criterion that affects the convergence speed. Moreover, Brissette et al. (2007) linked the mean precipitation amounts to an occurrence index to generate the daily precipitation amounts for solving the spatial intermittence problem mentioned by Wilks (1998). Brissette et al. (2007) also described the physical meaning and calculation method for the occurrence index, and divided the occurrence index into nine approximately uniform classes for fitting mean precipitation amounts. To avoid the data scattering problem, this study determines the number of classes of the occurrence index based on the number of precipitation events within each class. The maximum number of classes is set to 10; if the number of precipitation events in any class is lower than 30, this class will be combined with neighbouring classes until all classes contain at least 30 events. Moreover, this study also used the continuity ratio method, as proposed by Wilks (1998), to judge whether it can deal appropriately with the spatial intermittence problem.

2.3.2 Multi-site precipitation generation based on genetic algorithm (MGA)

The BFA is generally difficult to converge and a poor solution could be obtained if the iteration process is compulsorily terminated by using a predefined iteration number (Chen et al. 2014). In practice, searching an optimized correlation matrix for generating the spatially correlated multi-site precipitation series is an unconstrained nonlinear optimization problem. The objective is to minimize the error function derived from the difference between the observed and generated correlation matrix of precipitation series. The corresponding optimization model can be formulated as:

(4a) $\mathrm{M}\mathrm{i}\mathrm{n}f\left(c_{i,j}^{\mathrm{\Delta }}\right)$(4a)

subject to

(4b) $C_x=U\left(x_1,x_2,\dots ,x_{n\left(n-1\right)/2}\right)+L\left(x_1,x_2,\dots ,x_{n\left(n-1\right)/2}\right)+I_n$(4b)

(4c) $C_{\mathrm{g}\mathrm{e}\mathrm{n}}^x=\mathrm{B}\left(C_x\right)$(4c)

(4d) $\mathrm{\Delta }C=C_{\mathrm{o}\mathrm{b}\mathrm{s}}-C_{\mathrm{g}\mathrm{e}\mathrm{n}}^x$(4d)

(4e) $c_{i,j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\le c_{i,j}^x\le 1,\mathrm{\forall }i,j$(4e)

(4f) $c_{i,j}^{\mathrm{\Delta }}\in \mathrm{\Delta }C$(4f)

(4g) $c_{i,j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\in C_{\mathrm{o}\mathrm{b}\mathrm{s}}$(4g)

(4h) $c_{i,j}^x\in C_x$(4h)

where $x_1,x_2,\dots ,x_{n\left(n-1\right)/2}$ are the decision variables, $U$ and $L$ are the strictly upper and lower triangular matrix established by decision variables, respectively; $I_n$ is an n-level unit matrix; $C_x$ is n-level candidate correlation matrix; $C_{\mathrm{g}\mathrm{e}\mathrm{n}}^x$ is the generated correlation matrix of precipitation occurrence or amount based on the candidate correlation matrix $C_x$; $\mathrm{B}\left(\cdot \right)$ denotes the method from Brissette et al. (2007) to estimate the $C_{\mathrm{g}\mathrm{e}\mathrm{n}}^x$, which needs three main steps: (i) producing spatially correlated random numbers with uniform distribution based on the $C_x$, (ii) importing it into a single-site weather generator to generate spatially correlated precipitation occurrence and amounts, and (iii) calculating the correlation matrix of synthetic precipitation occurrence or amounts (i.e. $C_{\mathrm{g}\mathrm{e}\mathrm{n}}^x$). Constraint (4c) implies that $C_{\mathrm{g}\mathrm{e}\mathrm{n}}^x$ can be generated with the inputs of $C_x$. Constraint (4d) indicates that $\mathrm{\Delta }C$ is the error matrix between $C_{\mathrm{o}\mathrm{b}\mathrm{s}}$ and $C_{\mathrm{g}\mathrm{e}\mathrm{n}}^x$. Constraint (4e) provides the boundary conditions for decision variables, which indicates that the elements of $C_x$ should be more than those of $C_{\mathrm{o}\mathrm{b}\mathrm{s}}$ and less than 1. If $n=4$, $C_x$ can be written as:

(11) $C_x=\left[\begin{array}{c}\begin{array}{ccc}1& x_1& \begin{array}{cc}{x}_2& x_3\end{array}\\ x_1& 1& \begin{array}{cc}{x}_4& x_5\end{array}\\ x_2& x_4& \begin{array}{cc}1& x_6\end{array}\end{array}\\ \begin{array}{ccc}{x}_3& x_5& \begin{array}{cc}{x}_6& 1\end{array}\end{array}\end{array}\right],\mathrm{i}\mathrm{f}n=4$(11)

This unconstrained nonlinear optimization problem can be solved by GA. In a GA, the population of candidate solutions to an optimization problem evolves toward better solutions by mimicking the natural selection of biological evolution (Whitley 1994). GA is one of the adaptive heuristic search algorithms, and requires adjustment of a number of parameters, such as population size, mutation rate and crossover rate (Karahan et al. 2007). Figure 2 also presents the searching process of GA. Firstly, a set of individuals is generated randomly within the bounds of decision variables (i.e. constraint (4e)) to initialize the population. Each individual is then transformed into a matrix via constraint (4b). Afterwards, each population contains various combinations of candidate correlation matrixes ($C_x$), and it will be fed into constraints (4c) and (4d) in sequence to calculate the value of the objective function and evaluate the difference between the observed and synthetic correlation matrixes of the precipitation series. As a result, the fitness score (the same as the value of the objective function) can be obtained. Such process will repeat through three main genetic operations, i.e. selection, crossover and mutation, until the termination criteria is satisfied. The relationship between selection, crossover and mutation operations, as well as the mutation and crossover rates, are described in Rajesh et al. (2010). When the GA stops, a close-to-optimal spatially correlated multi-site precipitation series will be identified.

2.4 Approach design

In this study, we adopt five types of error function: (i) maximum absolute error (MAE), (ii) sum of squared error (SSE), (iii) average absolute error (AAE), (iv) root mean square error (RMSE), and (v) sum of mean and standard-deviation of absolute errors (SMSAE). The specific equations are listed as follows:

(5a) $\mathrm{M}\mathrm{A}\mathrm{E}:f\left(c_{i,j}^{\mathrm{\Delta }}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|c_{i,j}^{\mathrm{\Delta }}\right|\right\}$(5a)

(5b) $\text{SSE}:f\left(c_{i,j}^{\Delta }\right)={\displaystyle \sum _n^{i=1}{\displaystyle \sum _n^{j=1}{c}_{i,j}^{\Delta }{ }^2}}$(5b)

(5c) $\mathrm{A}\mathrm{A}\mathrm{E}:f\left(c_{i,j}^{\mathrm{\Delta }}\right)=E\left[\left|c_{i,j}^{\mathrm{\Delta }}\right|\right]$(5c)

(5d) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}:f\left(c_{i,j}^{\mathrm{\Delta }}\right)=\sqrt{\frac{1}{n^2}\sum _n^{i=1}\sum _n^{j=1}c{{ }_{i,j}^{\mathrm{\Delta }}}^2}$(5d)

(5e) $\mathrm{S}\mathrm{M}\mathrm{S}\mathrm{A}\mathrm{E}:f\left(c_{i,j}^{\mathrm{\Delta }}\right)=E\left[\left|c_{i,j}^{\mathrm{\Delta }}\right|\right]+\mathrm{S}\mathrm{t}\mathrm{d}\left[\left|c_{i,j}^{\mathrm{\Delta }}\right|\right]$(5e)

where $\mathrm{E}\left[\cdot \right]$ and $\mathrm{S}\mathrm{t}\mathrm{d}\left[\cdot \right]$ denote computations of statistical mean and standard deviation, respectively. In this study, three types of stochastic weather generator are to be compared. They are: (i) WGEN, (ii) MBFA, which drives the WGEN with spatially correlated random numbers identified based on the BFA, and (iii) MGA which drives the WGEN with spatially correlated random numbers identified based on GA. Corresponding to Equations (5a)–(5e), MGA is further divided into five types, namely MGA1 (based on MAE), MGA2 (based on SSE), MGA3 (based on AAE), MGA4 (based on RMSE), and MGA5 (based on SMSAE), Table 3 shows the differences among the various stochastic weather generators. The outputs from these stochastic weather generators will be fed into the SWAT model for evaluating the effects of synthetic precipitation on the simulated streamflow (Fig. 2).

Table 3. Features of various stochastic weather generators.

Feature of SWG	WGEN	MBFA	MGAs
			MGA1	MGA2	MGA3	MGA4	MGA5
Is a single-site SWG	√	×	×	×	×	×	×
Is a multiple-site SWG	×	√	√	√	√	√	√
Uses Markov chain to generate the precipitation occurrence	√	√	√	√	√	√	√
Uses gamma distribution to generate the precipitation amount	√	√	√	√	√	√	√
Uses occurrence index to solve the spatial intermittence problem	×	√	√	√	√	√	√
Method to search a satisfied correlation matrix	○	BFA	GA	GA	GA	GA	GA
Objective function of GA	○	○	MAE	SSE	AAE	RMSE	SMSAE

√, × and ○ denote yes, no and none, respectively.

3 Results and discussion

3.1 SWAT model calibration and validation

Table 2 shows the calibrated parameters of the SWAT model. We set the 11-year dataset from 1980 to 1993 (excluding 1981, 1982 and 1990 due to lack of part of meteorological data) as the calibration period, and the 11-year dataset from 1994 to 2004 as the validation period. The observed and simulated daily streamflow for the calibration and validation periods are shown in Fig. 3(a) and (b), respectively. It is indicated that the SWAT model would somewhat underestimate the runoff of some peak-flood events, although most simulated runoff closely matches the observed runoff for both calibration and validation periods. It may be the case that (i) the SWAT model did not consider the effect of the Gleniffer Reservoir in the upstream due to a lack of reservoir operation records, and (ii) the calibration efforts place too much emphasis on the overall accuracy instead of the peak flows. Overall, the values of $E_{\mathrm{N}\mathrm{S}}$ in the calibration and validation periods reached 0.52 and 0.57, respectively. This performance is considered satisfactory for this study, as our focus is on evaluating the effect of the spatial correlation of rainfall on hydrological flows.

Figure 3. Observed vs simulated daily streamflow for (a) the calibration period, and (b) the validation period.

3.2 Identification of the best error functions as the objective of MGAs

Figure 4 shows the Taylor diagrams for monthly inter-station correlations calculated by observed and six multi-site SWGs. Figure 4(a) and (b) shows the precipitation occurrence and amounts, respectively. The Taylor diagram plots three statistics to quantify how realistically each SWG reproduces the observations with similar inter-station correlations: (i) the correlation coefficient (Taylor-CC) related to the azimuthal angle; (ii) the centred root mean square difference (Taylor-CRD) proportional to the distance between each multi-site SWG represented by different marks on the diagram and the point labeled “observed”; and (iii) the standard deviation (Taylor-SD) proportional to the radial distance from the origin point.

Figure 4. Taylor diagrams displaying a statistical comparison of the monthly inter-station correlations calculated by the observed and generated precipitation: (a) occurrence and (b) amounts.

It can be seen from Fig. 4(a) that the values of Taylor-SD of all the multi-site SWGs are higher than that of the observed record. For MBFA, the values of Taylor-CC, Taylor-CRD and Taylor-SD are approx. 0.68, 0.059 and 0.078, respectively. All MGAs have a higher Taylor-CC and a lower Taylor-CRD than the MBFA, implying a better performance of MGAs over MBFA in generating precipitation occurrence. The MGA2, MGA3 and MGA5 models are found to have slightly higher Taylor-CC values, and MGA1 and MGA4 have low Taylor-SD values. In Fig. 4(b), the values of Taylor-CC, Taylor-CRD and Taylor-SD for MBFA are approx. 0.52, 0.106 and 0.125, respectively. The MGA3 and MGA5 models have the highest Taylor-CC values, but MGA3 shows higher variations (i.e. Taylor-SD) and greater distance from the observed point (i.e. Taylor-CRD) than MGA5. It is also noted that MGA3 has similar Taylor-SD and Taylor-CRD levels as compared to MBFA, but a slightly higher Taylor-CC value.

Fig. 5 shows boxplots of errors (both residual and absolute error) for the monthly inter-station correlations of the generated precipitation occurrence and amounts. Figure 5(a) and (c) are based on precipitation occurrence, and Fig. 5(b) and (d) are based on precipitation amounts. In the boxplots, the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively, the difference between 75th and 25th percentiles are defined as interquartile range (IQR), and the central black line within the box indicates the median. The whiskers above and below the box extend to the most extreme data without considering outliers, and the outliers are marked as black points above and/or below the whiskers. It was found that MGA2 and MGA4 in Fig. 5(a) and MGA1 and MGA2 in Fig. 5(b) may have an overestimation problem, as the majority of their boxes are higher than zero. Moreover, MBFA, MGA1, MGA3 and MGA5 in Fig. 5(a) have similar ranges and locations of boxes and whiskers. From Fig. 5(b), MGA3 is identified to be superior to the other three SWGs since it has the smallest IQR with a close-to-zero median. In terms of absolute error for rainfall occurrence (Fig. 5(c)), the five MGAs perform at relatively similar levels but they are generally superior to MBFA due to their lower median and narrower IQR and whiskers. From Fig. 5(d), MGA3 is also seen to be the best alternative, judging from its lowest median, percentile and whisker levels. The selection of the fitness function is crucial for the search performance of the GA (Chlumecký et al. 2017). The reason that the MGA3 shows the best performance might be that the average absolute error is the most suitable fitness function of GA to search the optimized correlation matrix for generating the spatially correlated multi-site precipitation series in this basin.

Figure 5. Boxplots of residual errors for the monthly inter-station correlations of the generated precipitation: (a) occurrence and (b) amounts, and absolute errors for the correlation coefficient of the generated precipitation: (c) occurrence and (d) amounts.

3.3 Performance of MGA3, MBFA and WGEN

The results shown in the Taylor diagrams and boxplots demonstrate that MGA3 (based on AAE error function) is the best option among all MGAs in serving as the objection function for MGAs to produce spatially correlated precipitation series. In the following sections, we compare the differences among MGA3, MBFA and WGEN. Figure 6 presents scatterplots of the observed vs generated inter-station correlations of precipitation occurrence and amounts for all station pairs over the Red Deer River watershed. The correlation of precipitation occurrence is computed based on a monthly scale, whereas the correlation of precipitation amounts is calculated at a seasonal scale. This is because (i) MGA3, MBFA and WGEN work on a monthly basis to estimate the precipitation transition probabilities of Markov chains to generate the occurrence, and (ii) for MGA3 and MBFA, a seasonal link between the occurrence index and precipitation mean value can be used for avoiding the spatial intermittence problem in the precipitation amount generation process. In Fig. 6, most inter-station correlations of occurrence and amounts generated by WGEN would be close to zero indicating WGEN’s incapability in reproducing inter-station correlations. Both MBFA and MGA3 can reasonably reproduce the observed monthly inter-station correlations of both precipitation occurrence and amounts, but MGA3 seems to have a closer scattering of points along the 1:1 benchmark line than MBFA, which is consistent to the results in Fig. 5.

Figure 6. Observed vs generated inter-station correlations for the precipitation: (a) occurrence and (b) amounts.

Figure 7. Converging processes of MBAF and MGA3 in generating precipitation: (a) occurrence and (b) amounts.

The basic statistics of the observed against generated precipitation series for all seasons and station pairs are shown in the Supplementary material (Fig. S2(a)–(c)). It is indicated that the three SWGs have similar performance in terms of mean, standard deviation and maximum daily precipitation. This is mainly because all generators use similar techniques in both precipitation occurrence and amount generation. It is also found that there is overestimation of the larger values of mean precipitation. In this study, we used the default gamma distribution to simulate the precipitation amounts. For a comparison of the fitting performance between the gamma distribution and the other three distributions (i.e. the exponential, Weibull and the generalized Pareto distributions) the reader is referred to Figs. S3–S6 in the Supplementary material, while Fig. S7 shows the basic statistics of the observed against generated precipitation series using these four distributions. We found that the performance of the four distributions in fitting the observed precipitation data is acceptable, but there is still bias in the mean and maximum values of the generated precipitation data. The performance of the exponential, Weibull and generalized Pareto distributions is not notably superior to that of the gamma distribution. The problem of overestimation of the larger values of mean precipitation may be caused by the stochastic nature of precipitation generation of this basin. The continuity ratios of the observed against generated precipitation series for all seasons and all combinations of station pairs are presented in Fig. S2(d). The results indicate that both MBFA and MGA3 could well reproduce the observed continuity ratios implying their capability of addressing spatial intermittence. But MGA3 seems to have relatively close scatter plots to the 1:1 line than MBFA. WGEN obviously performs the poorest in reproducing the observed continuity ratios as it does not take the spatial intermittence into consideration.

For the MBAF method, we set the convergence criterion, $\eta$, to 0.1. Two types of termination criteria are adopted, either when the maximum iteration number (i.e. $\tau =1000$) is reached, or when a predefined error tolerance (i.e. $\rho ={10}^{-3}$) is exceeded. For MGA3, the population size, elite count, crossover rate and mutation rate for GA operations are configured to 2000, 200, 0.85 and 0.01, respectively; the corresponding termination criteria is based on the maximum generation number (set to 1000). Figure 7 shows the search processes of MBFA and MGA3 in generating the precipitation occurrence and amounts. For further comparison of MBFA to other MGAs (i.e. MGA1, MGA2, MGA4, and MGA5) in searching, the reader is referred to the Supplementary material (Fig. S8). For both MBFA and MGA3, the overall trends of AAE would decrease and converge, with rapid fluctuations throughout the processes. Obviously, the magnitude of such fluctuation is more significant for MBFA, especially regarding the precipitation amounts. This implies that the searching process of MBFA exhibits a higher level of uncertainty and easily leads to a poorer solution when terminated by the maximum number of iterations. MGA3, on the other hand, could notably alleviate such a problem. Generally, Fig. 7 demonstrates that MGA3 may show a slower convergence speed, but could reach a solution with better objective value. This is consistent with the findings from Figs. 5 and 6, where the inter-station correlations could be better regenerated through MGA3. In terms of hardware and software setting, MBAF and MGA3 (on MATLAB platform) were run on a workstation with 96 GB RAM and AMD Opteron 2.6 GHz 16 Core 32 logical processor. We used the Parallel Computing Toolbox of MATLAB and executed the MGA3 on a parallel pool of 30 works. Table 4 shows the computational costs of MBFA and MGA3. It took a total of 24.76 h for MGA3, including 2.43 h generating the precipitation occurrence and 22.33 h generating the precipitation amount. Compared with MBFA, MGA3 is more computationally intensive obviously. Although the 24.67 h calculation time is accepted in this study and such computing time could be further shortened if more CPU cores are available, the large computational burden is the conspicuous drawback of MGAs.

Table 4. Computational costs of MBFA and MGA3.

	MBFA (h)	MGA3 (h)
Generating the precipitation occurrence	0.01	2.43
Generating the precipitation amount	0.08	22.33

MBAF and MGA3 (on MATLAB platform) were run on a workstation with 96 GB RAM and AMD Opteron 2.6 GHz 16 Core 32 Logical Processor. We used the Parallel Computing Toolbox of MATLAB and executed the MGA3 on a parallel pool of 30 works.

3.4 Hydrological responses to different SWGs

To evaluate the effect of SWGs on hydrological processes, the three SWGs (i.e. WGEN, MBFA and MGA3) are adopted to generate 200-year precipitation series that are used as weather input to drive the SWAT model for simulating runoff. Rather than selecting the observed streamflow directly from the hydrometric station (05CK004), we use the simulated runoff from running SWAT model using 24-year observed precipitation data as the benchmark. This is for the purpose of reducing the uncertainty/error from hydrological modelling and carrying out a fairer comparison among generators. Figure 8 shows the mean and standard deviation of monthly streamflow simulated by SWAT using the observed and generated precipitation series. The mean monthly streamflow would be accurately represented by the three SWGs data during January to June and September to December but overestimated during July and August (Fig. 8(a)). This is probably because of the overestimation of the larger precipitation values during the summer period (see Fig. S2). In terms of standard deviation, MGA3 is found to perform well in most of the months except for May, while MBFA seems to have slight overestimation from April to December. WGEN considerably underestimates the monthly flow variability especially from May to July; the main reason is that WGEN generates precipitation series without considering any spatial dependence, which may attenuate the peak precipitation values.

Figure 8. Statistics of (a) mean and (b) standard deviation of monthly streamflow simulated by SWAT using observed and generated precipitation series. OBS: observed precipitation series.

Figure 9 shows the frequency curves of annual maximum flows simulated by SWAT using observed 24-year and generated 200-year precipitation series. The frequency curves are constructed by using the Pearson type III distribution (Koutrouvelis and Canavos 1999), as it has the most satisfactory fitting performance compared to the extreme value distribution, Weibull distribution and exponential distribution. For detailed fitting performance, the reader is referred to the Supplementary material (Figs. S9 and S10). It is found that WGEN significantly underestimates the extreme flows, while MBFA and MGA3 have slight underestimation. This is somewhat consistent with Fig. 8(b) where high flows in one location of the watershed may be offset by low flows in other locations due to independent precipitation series generated from WGEN. As the inter-station correlations are reasonably reproduced for the observed rainfall time series, the underestimation of rainfall extremes using the gamma distribution most likely leads to the underestimation of the extreme flows by both MBFA and MGA3. In addition, compared with MBFA, the frequency curve of MGA3 is closer to the observed one, indicating that MGA3 could reproduce the observed extreme flows with a smaller bias. Overall, the multiple-site SWGs are obviously more advantageous than single-site ones in simulating extreme flows.

Figure 9. Annual maximum flows over different return periods using observed and generated precipitation series. OBS: observed precipitation series.

4 Conclusion

In this research work, we compared the performance of GA with the BFA in the framework of Brissette et al. (2007) in regard to finding suitable correlation matrices and random numbers for multi-site SWG. We designed five types of error function, MAE, SSE, AAE, RMSE and SMSAE, as the alternative objectives of GA, and found that MGAs with the AAE (i.e. MGA3) could reproduce the observed precipitation series with the best spatial correlations. It should be noted that the selection of error functions would affect the fitness score and a suitable error function should be tested for specific watersheds. We compared the performance of MGA3, MBFA and WGEN in the simulation of precipitation series of four stations distributed over the Red Deer watershed, Alberta, Canada. Because the three SWGs adopted similar techniques, namely a two-state Markov chain coupled with a gamma distribution in generating precipitation series, the results exhibit similar basic statistics, i.e. mean, standard deviation and maximum value. However, all the three generators somewhat overestimated the larger values of mean precipitation, due to statistical fitting and stochastic sampling. Furthermore, by linking average precipitation to the occurrence index as suggested by Brissette et al. (2007), both MBFA and MGA3 reproduced the observed continuity ratio well. As multiple-site SWGs, MBFA and MGA3 both well reproduced the observed monthly inter-station correlations of precipitation occurrence and amounts. Through scatterplots and Taylor diagram analysis, the MGA3-generated inter-station correlations of precipitation occurrence and amounts turned out to be closer to the observed records than MBFA-generated ones. This is due to the better searching capability of Gas for the correlation matrix. It was also found that the larger magnitude of the fluctuations of MBFA compared with MGA3 would lead to a higher level of uncertainty in searching suitable solutions.

After validation of the SWAT model for the Red Deer watershed, we assessed the hydrological responses to the precipitation series generated from MGA3, MBFA and WGEN. Using the data from three SWGs, the SWAT model accurately represented the monthly mean streamflow for most of the months except for some overestimation in July and August. For the standard deviation, MGA3 performed well in almost all months except for May and MBFA slightly overestimated the monthly streamflow variability from April to December. However, WGEN considerably underestimated the variability of monthly flows especially from May to July. Compared with MBFA, MGA3 reproduced the observed extreme flow with smaller bias. WGEN significantly underestimated the annual maximum flows due to its poor capability in addressing precipitation correlations. In general, spatial correlation was found essential in generating precipitation series for multiple sites in the studied Red Deer watershed, for a better representation of flow characteristics and flood peaks.

The main contributions of this research work are: (i) comparison of the performance of GA with BFA in the framework of the Brissette et al. (2007) multiple-site SWG to search for a global optimal solution for generating the spatial correlated precipitation series; and (ii) the assessment of the effects of the multiple-site SWG on hydrological responses in the Red Deer River watershed, which is characterized by large size and humid continental climate. The proposed method has a more stable convergence than the framework proposed by Brissette et al. (2007) and could be used as a viable alternative in regenerating the spatial correlation, basic statistics and extreme data of weather variables. However, the method was compromised by a higher computational need due to GA iterations. Based on our test, GA could work well with less than 10 stations with manageable amount of time; but when the number of stations increases, the computational burden may unavoidably increase. Setting up a MATLAB distributed computing server to run MGA3 on computer clusters, clouds, and grids might be a potential solution.

Acknowledgements

The authors are grateful to the associate editor and the reviewers for their insightful comments and suggestions, which have greatly helped to improve the paper.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary Material

Supplemental data for this article can be accessed here.

Additional information

Funding

This project was supported by Start-Up Grant [M4081327.030] from School of Civil and Environmental Engineering, Nanyang Technological University, Singapore.