Multisite Spatiotemporal Streamflow Simulation-With an Application to Irrigation Water Shortage Risk Assessment

Regional water resources management generally requires knowledge of multisite streamflows which exhibit random, yet spatially and temporally correlated, variabilities. The complexity of such correlated randomness makes decision-making for water resources management a difficult task. With presence of uncertainties in space and time, risk-based decision making using stochastic models is sought after. In this study we propose a spatiotemporal stochastic simulation model for multisite streamflow simulation. The model is composed of three components: (1) stochastic simulation of bivariate non-Gaussian distributions, (2) anisotropic space-time covariance function which characterizes the spatial and temporal variations of multisite ten-day periods (TDP) streamflows, and (3) Monte Carlo spatiotemporal simulation of streamflows. The model was applied to the Chia-Nan Irrigation District in southern Taiwan for a multisite spatiotemporal ten-day-period streamflow simulation. Through a rigorous evaluation, the proposed spatiotemporal model is found capable of preserving not only the marginal distributions but also the spatiotemporal correlation structure of the multisite streamflows. An example application which demonstrates utilization of the proposed model for irrigation water shortage risk assessment is also presented.


INTRODUCTION
Streamflow data are widely used in water resource management, hazard risk assessment, modeling climate change scenarios and other applications.Many applications require streamflow estimates at different gauging stations for regional-based decision making.Unfortunately, many gauging stations do not have flow observations with long record lengths to support meaningful analysis and modeling required by the target applications.In order to circumvent such problem, stochastic simulation (or Monte Carlo simulation) of multisite streamflows is an essential need in water resources management.Another factor exacerbating the simulation difficulties is that streamflows at different sites exhibit not only spatial but also temporal correlations, necessitating a spatiotemporal stochastic simulation model for multisite streamflows simulation.
An essential element of stochastic simulation is that statistical properties of the phenomenon under investigation must be preserved in the simulation outputs.While it may be difficult to preserve all statistical properties of complex phenomenon, the minimum goal is to preserve up to the secondmoment properties, i.e., the expectation and variance/covariance.Many environmental and natural variables, such as precipitation amount, storm duration, and streamflow exhibit significant degrees of non-Gaussian randomness (Parada and Liang 2010;Cheng et al. 2011).Univariate simulation of such non-Gaussian random variables can be achieved by using the probability integral transformation (Cheng et al. 2007), and many statistical software packages provide builtin commands for such univariate simulations.Validation of the results of univariate simulation only considers the properties (for examples, the unbiasedness and mean squared error) of parameter estimators (Cheng et al. 2007).Stochastic simulation of two correlated random variables, commonly known as the bivariate simulation, requires preservation of the parameters of marginal densities and the correlation of the two random variables as well.Lee et al. (2010) conducted a bivariate storm frequency analysis by considering a bivariate Gumbel distribution for storm duration and depth.Cheng et al. (2011) developed a frequency-factor based approach for stochastic simulation of the bivariate gamma distribution.
For cases that involve simulation of more than two interrelated random variables, covariance between random variables also need to be preserved, in addition to parameters of individual random variables.A covariance matrix, which characterizes correlations of all possible pairs of two random variables, plays an important role in multivariate simulation.Wu et al. (2012) assessed climate change impacts on basin-average annual typhoon rainfalls by simulating multi-site typhoon rainfalls with consideration of their multisite correlations.Wang and Ding (2007) developed a nonparametric multivariate kernel density model for multivariate daily flow generation.
Stochastic simulations involving spatial and temporal variations have gained more attention in recent years.Westra and Sharma (2009) developed a multivariate estimation model, based on independent component analysis, for multiple reservoir inflow estimations.Lee (2012) proposed a stochastic approach for multivariate streamflow simulation by decomposing multivariate time series into independent components.In this study we developed a stochastic spatiotemporal model for multisite streamflow simulation.The model aims to generate streamflows at different sites which can then be used for drought risk assessment of mitigation measures.Key elements of the model include stochastic simulation of bivariate non-Gaussian distributions and an anisotropic space-time covariance function which characterizes the spatial and temporal variations of multisite streamflows.
The paper is organized as follows.Section 2 describes the study area and streamflow data used in this study.Section 3 gives details of the methods for characterizing the spatial and temporal covariances of streamflows at different sites.Section 4 introduces an algorithm for stochastic simulation of bivariate non-Gaussian distributions using a frequency-factor based approach.Procedures for extending from a bivariate simulation to a multivariate simulation are also described in section 4. Section 5 shows the results of stochastic spatiotemporal simulations using the proposed model.Validation of the simulation results is also presented in section 5. Section 6 demonstrates an application in irrigation management using the simulation results.A few concluding remarks are given in section 7.

STUDY AREA AND DATA
The Chia-Nan Irrigation District (CNID) and its neighboring areas (see Fig. 1) in southern Taiwan were selected for our study area.The study area contains seven sub-basins with twelve streamflow gauging stations.In Taiwan, flows of nominal ten-day periods (TDP) are used for irrigation management, with a total of 36 TDP in each year.Twenty six years (from 1975 to 2000) of TDP streamflow data available at the twelve stations were used for this study.Figure 2 shows the long-term average of TDP flows at individual flow stations.A distinct pattern can be observed with high flow season (May -October) accounting for a majority of the annual total flows.

CHARACTERIZING MULTISITE STREAMFLOWS
TDP streamflow data are available at different gauging stations.The drainage areas for individual gauging stations vary from 83 km 2 for station 4 to 812 km 2 for station 12.As a result, the magnitudes and statistical properties (for examples, the mean and standard deviation) of TDP streamflows differ from one station to another across the study area.The heterogeneity in the expected value and variance of site-specific TDP flows needs to be considered in multisite streamflow simulation.
Let Q i, j, k represent streamflow of the j-th TDP of the i-th year at the k-th station.Site-specific TDP flows are standardized with respect to their expected values and standard deviations through the following equation: where m j, k and s j, k are respectively the mean and standard deviation of streamflows of the j-th TDP at the k-th station, and Q , , * i j k is the standardized TDP streamflows at the k-th station.Standardized TDP flows at different stations have common zero expectation and unit standard deviation.We consider the spatial variation of the standardized TDP flows as an isotropic stationary Pearson Type III (PT3) random field with the following marginal probability density function: where X represents the random variable (i.e., the standardized TDP flows in this study), a , b, and f are respectively The standardized TDP flows at different stations are spatial and temporally correlated and hence their covariances in the spatial and time domains need to be considered simultaneously for a full characterization of the multisite TDP streamflows.Semi-variograms have been used in many studies for modeling of spatial variations of natural and hydrological phenomena (Cheng et al. 2000;Cheng et al. 2003;Franco et al. 2006).A semi-variogram is a function that describes the spatial correlation structure of a random field.Let X(s) be a random variable at location s and { ( ), } X s s !X be an isotropic stationary random field over a spatial domain X.In geostatistics, the characteristics of spatial variation of a random field is often expressed in terms of the semi-variogram ( ) s s cl defined as: where E(X) represents the expected value of a random variable X.For a stationary random field, the semi-variogram is independent of the locations (s and sl ) and can be expressed by ( ) h c with h being the distance between s and sl .Details of the properties of the semi-variogram and its modeling can be found in Journel and Huijbregts (1978).
In this study we first established the spatial semi-variogram and temporal semi-variogram of the standardized TDP flows using only the same-TDP and same-site standardized TDP flows, respectively.Model fitting of both the spatial and temporal semi-variograms of the standardized TDP flows using an exponential model yields the following results: where ( ) s S c D and ( ) t T c D are respectively the spatial and temporal semi-variograms, and Δs and Δt are measured in units of kilometers and ten-day periods, respectively.The covariance functions of the standardized TDP flows in space and time (i.e., C S and C T ) are respectively related to the temporal and spatial semi-variograms and can be represented by the following equations: In modeling the temporal semi-variogram, we consider the time domain (measured in units of TDPs) as a one dimensional space, and the temporal semi-variogram was established in a manner similar to spatial semi-variogram modeling.
The above spatial and temporal semi-variograms (or covariance functions) characterize the random co-variations of the standardized TDP flows in the spatial domain and time domain, respectively.However, co-variation can also exist between standardized TDP flows of different TDPs and at different sites.For example, the standardized TDP flows of the 15 th TDP at station 1 and the standardized TDP flows of the 16 th TDP at station 2 (which is near station 1) are likely to be correlated.However, such co-variation cannot be characterized by either the spatial or temporal semi-variogram (or covariance function) alone.In order to account for such spatiotemporal co-variation (or the spacetime interaction) we consider the spatiotemporal random field being anisotropic in time and in space, and adopted the following spatiotemporal anisotropic covariance function of the standardized TDP flows: where ( , ) C s t , S T D D is the spatiotemporal covariance function and k is a scale factor, also known as the anisotropic ratio, with a value of k = 2.38/54.6.The above spatiotemporal covariance function characterizes isotropic covariance functions in space and time, but anisotropic covariance in space-time (Serban 2013).

SPATIOTEMPORAL SIMULATION OF TDP STREAMFLOWS
In the previous section we demonstrated that the multisite standardized TDP flows can be modeled as a PT3 spatiotemporal random field with an anisotropic covariance function in space-time domain.In this section we describe details of the algorithm and procedures for spatiotemporal simulation of multisite TDP streamflows.
The task of spatiotemporal simulation in this study is to generate a large set of TDP streamflow samples.Each TDP streamflow sample represents one year (36 TDPs) of TDP streamflows at n different flow stations, as demonstrated in Fig. 3a.A spatiotemporal sample is composed of 36 × n standardized TDP flows and can be expressed as X(s, t), s = {s 1 , …, s n }; t = {t 1 , …, t m }.Occurrence of such a sample is characterized by the following joint density: [ ( , ), ( , ), , ( , ), ( , ), ( , ), , ( , ), , ( , ), Based on the conditional probability, the above joint density can be further expressed as: Therefore, simulation of a spatiotemporal random sample can be achieved by a sequential conditional simulation approach.Figure 3b illustrates the column-preference (i.e., space preference) sequential conditional simulation algorithm.Conditional simulation of X(s p , t q ) given other X(s i , t j ) forms the foundation of the sequential simulation.Without loss of generality, simulation of X(s 3 , t 3 ) given X(s 1 , t 1 ) and X(s 2 , t 2 ) is described as follows.
The joint distribution {X(s 1 , t 1 ), X(s 2 , t 2 ), X(s 3 , t 3 )} is non-Gaussian since standardized TDP flows are PT3 distributed.A PT3 distribution can be transformed to a corresponding standard normal distribution using its frequency factor (Cheng et al. 2007): where X is a PT3 random variable with expected value X n , standard deviation X v and coefficient of skewness c , K is the frequency factor of the PT3 distribution, and U represents a standard normal random variable.Similarly, a PT3 spatialtemporal random field with spatiotemporal anisotropic covariance function of Eq. ( 8) can be transformed to a corresponding standard normal random field.Cheng et al. (2011) showed that such transformation requires the following conversion between the correlation coefficient of a bivariate PT3 distribution and the correlation coefficient of a bivariate standard normal distribution: t respectively being the correlation coefficients of the bivariate PT3 and bivariate standard normal distributions, and A 1 6 = `j .For convenience, {X(s 1 , t 1 ), X(s 2 , t 2 ), X(s 3 , t 3 )} is replaced by X = {X 1 , X 2 , X 3 } from here.The correlation matrix of X(Σ X ) can be converted to the correlation matrix of U = {U 1 , U 2 , U 3 }(Σ U ) using Eq. ( 13).The correlation matrix Σ X is calculated using Eq. ( 8) with consideration of the space-time anisotropy.From Eq. (10) we have Techniques for multivariate Gaussian simulation have been well documented and we used a sequential Gaussian simulation (Liou et al. 2011) for stochastic simulation of U = {U 1 , U 2 , U 3 }.Random samples of a standard normal random field can then be transformed to random samples of the spatiotemporal PT3 random field using Eqs.( 11) and ( 12).An illustrative flowchart of the proposed spatiotemporal TDP streamflow simulation approach is shown in Fig. 4.

SIMULATION RESULTS AND VALIDATION
In this study TDP streamflows were generated at eight stations (stations 1, 3, 4, 7, 8, 9, 10, and 11) which have water intake structures.A spatiotemporal sample consists of 36 TDP flows for each of the eight flow stations, making a total of 288 streamflows in each simulated sample.In order to compare sample statistics of the simulation results against that of the 26-year historical data, we generated 1000 blocks of TDP samples, with each block consisting of 26 spatiotemporal samples.Each block of simulated TDP samples can be viewed as a mimic (in a statistical sense) of the 26-year historical data.Using these 1000 blocks of simulated streamflows, we calculated the sample mean and standard deviation of the 26-year simulated TDP streamflows, and evaluated their uncertainties.For example, from each block of simulated samples, an estimate of the mean (or standard deviation) of the i-th TDP streamflows can be calculated for any station.A total of 1000 blocks yields 1000 estimates of the mean streamflow of the i-th TDP, and the mean and standard deviation of these 1000 estimates can be further calculated.flows calculated from 1000 blocks of simulated data was lower (less than 5%) than the standard deviation calculated from historical TDP flows during the high flow season, their differences are minor and will not result in any problems for practical applications.
In addition to the mean and standard deviation of TDP streamflows, it is also important to assess whether the spatiotemporal TDP streamflow simulation can preserve the correlation in time, space, and time-space.A spatiotemporal correlation matrix P is defined as: where P ij itself is a sub-correlation matrix of dimension 8 × 8 which represents the correlation between the i-th TDP flows and the j-th TDP flows, i.e., ( , ) Using 26 years of historical streamflows, we calculated the empirical spatiotemporal correlation matrix of TDP streamflows (see Fig. 6).A sample correlation matrix can also be calculated using each block of simulated samples, and thus the sample mean and standard deviation of the correla-tion matrix P can be derived from a total of 1000 blocks of simulated samples.Figure 7 demonstrates that the mean spatiotemporal correlation matrix based on 1000 blocks of simulated streamflows is almost identical to the theoretical spatiotemporal correlation matrix which can be derived by substituting Eq. ( 8) into Eqs.( 15) and ( 16).Standard deviation map of the spatiotemporal correlation matrix is also shown in Fig. 8.The diagonal elements (from the lower left corner to the upper right corner) have zero standard deviation since they are associated with correlation coefficients of the same site and same TDP streamflows.It can also be observed that moving from the diagonal element, the standard deviation of the spatiotemporal correlation generally increases with the space-time lag from the diagonal element, until the lag reaches approximately 48.Beyond a space-time lag of 48, standard deviations of the spatiotemporal correlation of TDP streamflows seem to stabilize at a value of approximately 0.20.This result can be explained as follows.
The space-time lag of 48 is equivalent to a time lag of 6 TDPs since P ij in Eq. ( 16) represents the correlation matrix of TDP streamflows at eight stations.Notwithstanding the lag in space, the spatiotemporal correlation coefficient [Eq.( 8)] drops to lower than 0.1 when the time lag between streamflows at any two stations becomes greater than 6 TDPs (see Fig. 9).With a very low spatiotemporal correlation coefficient (near zero), the asymptotic standard deviation of the sample correlation coefficient can be approximated by the following equations (Hooper 1958;Priestley 1982): where n represents the sample size (in our study, n = 26) for calculation of the sample correlation coefficient, t is the population correlation coefficient, and r v is the asymptotic standard deviation of the sample correlation coefficient.In our study the sample size n equals 26 and the asymptotic standard deviation of the sample correlation coefficient calculated by either of the two equations yields a value of 0.196.These results also suffice to demonstrate the capability of the spatiotemporal simulation model in preserving the spatiotemporal correlation structure of TDP streamflows.
Using the mean and standard deviation of the spatiotemporal correlation coefficient, we can construct the (mean ± standard deviation) range of the sample correlation   coefficients.Figure 10 demonstrates the sample spatiotemporal correlation coefficients in the first row of the empirical spatiotemporal correlation matrix (see Fig. 6), together with the corresponding mean and (mean ± standard deviation) range of the sample correlation coefficients.It can be seen that the mean profile based on 1000 blocks of simulated data is almost identical to the profile of theoretical spatiotemporal correlation coefficient.In addition, most of the sample spatiotemporal correlation coefficients fall within the (mean ± standard deviation) range.A few points that fall outside of the (mean ± standard deviation) range are found to be within the [mean ± 2 × (standard deviation)] range.
Although not shown in the paper, similar results were found for other profiles of the spatiotemporal correlation matrix.This suggests that the spatiotemporal correlation structure of the historical TDP streamflows can be reflected in the spatiotemporal TDP flow simulation results.

EXAMPE APPLICATION FOR IRRIGATION RISK ASSESSMENT
In order to demonstrate a potential application of the proposed spatiotemporal streamflow simulation technique, an example application for irrigation water shortage risk assessment was carried out.Irrigation water for an irrigation subgroup (see Fig. 11 and the red-marked polygon in Fig. 1) within the CNID is provided by the Bah-Jon River and the Tzi-Shuei River.The irrigation subgroup is composed of 22 subdivisions with irrigation demand varying from 0.02 to 0.3 million cubic meters per TDP.Among the 22 subdivisions, subdivisions A and B were chosen to dem-onstrate irrigation risk assessment based on the results of spatiotemporal TDP streamflow simulation.The irrigation subgroup often experiences irrigation water shortage during the dry season and a mitigation measure using groundwater from a total of 18 groundwater wells (see Fig. 11) has been investigated.The mitigation measure plans to withdraw groundwater at a rate of 0.01 million cubic meters per TDP from each of the 18 groundwater wells.Calculations of the irrigation water amount that can be provided to individual subdivisions within the irrigation subgroup are conducted by an irrigation management model.Within an irrigation subgroup, individual irrigation subdivisions are prioritized by considering their locations in the irrigation water supply network.A general criterion of the irrigation management model is to supply 70% of the irrigation water demand (i.e., the initial supply) to as many subdivisions as possible first, and then the additional amount of water, if existing, is supplied to the subdivisions with higher priorities.The rationale of allocating the initial supply is that with this amount of irrigation water, paddy field rice can endure water stress without significant adverse effects.A detailed description of the irrigation management model is out of the scope of this study and readers are referred to Wen et al. (2007) and for details of the irrigation management model.
One thousand runs of spatiotemporal streamflow simulation were conducted for this example application study.Each simulation run yields a spatiotemporal sample consisting of 36 TDP streamflows at each flow station.Using the results of each simulation run, the irrigation management model calculated the amounts of irrigation water that can be provided to individual subdivisions under two different (1, ) (1, )]; , , } j j j j 1 36 = .
scenarios -(1) without groundwater withdraw and (2) with groundwater withdraw.Under scenario 2, available surface water from a reservoir and rivers are considered as the primary source of irrigation water, i.e., irrigation water drawn from groundwater wells will be provided only when surface water is insufficient for irrigation water supply.
Consider the irrigation water supply (including surface and groundwater) and irrigation demand of a subdivision (for example, subdivision A or B).For the t-th TDP, let X d (t) and X s (t) respectively represent the amounts of irrigation demand and irrigation water supply.From the results of one spatiotemporal TDP streamflow simulation run, the ratio of irrigation shortage (hereinafter referred to as the shortage ratio) is calculated as: Not does the shortage ratio vary with TDPs in one year but also from one year (simulation run) to another.Thus, from the results of 1000 spatiotemporal simulation runs, the mean and empirical cumulative distribution function (ECDF) of each TDP-specific shortage ratio can be obtained.
Figure 12 demonstrates the TDP-specific average shortage ratios and the ECDF of shortage ratios of the 7 th TDP for subdivisions A and B under different scenarios.Average shortage ratios are higher within the 4 th -14 th TDP period (the first-crop period) during which streamflows are low (see Fig. 2) and large amount of irrigation water is needed for paddy field leveling and rice planting.
Under scenario 1, (i.e., without implementation of the groundwater mitigation measure), the mean of average shortage ratios for subdivisions A and B during the 4 th -14 th TDP period are 0.225 and 0.18, respectively (see Fig. 12a).The probabilities for subdivisions A and B to be fully supplied with their irrigation water demands (i.e., zero shortage ratio) in the 7-th TDP are 0.22 and 0.38 (point 1 in Figs.12b and c), respectively.This is equivalently to say that the risk of irrigation water shortage in the 7-th TDP for subdivisions A and B are 0.78 and 0.62, respectively.Similarly, the probability for the 7-th TDP shortage ratio exceeding 0.3 for subdivisions A and B are 0.05 and 0 (point 3 in Figs.12b and c), respectively.
Under scenario 2, the mean of average shortage ratios for subdivisions A and B during the 4 th -14 th TDP period reduce to 0.125 and 0.075, respectively.The probabilities for subdivisions A and B to be fully supplied with their irrigation water demands in the 7-th TDP become 0.46 and 0.72 (point 2 in Figs.12b and c), respectively.The probability for the 7-th TDP shortage ratio exceeding 0.3 becomes 0 for both subdivisions.
The above results demonstrate that the groundwater mitigation measure can help reduce the average irrigation shortage of the first-crop period by approximately 10% of the irrigation water demand (0.225 -0.125 and 0.18 -0.075) for both subdivisions.The mitigation measure can also reduce the risk of irrigation water shortage in the 7-th TDP by 0.24 (0.46 -0.22) for subdivision A and 0.34 (0.72 -0.38) for subdivision B.

SUMMARY AND CONCLUSIONS
In this study we propose a spatiotemporal stochastic simulation model for multisite streamflow simulation.Through a rigorous evaluation, the model is capable of preserving not only the marginal distributions but also the spatiotemporal correlation structure of the multisite streamflows.A few concluding remarks are given here: (1) The spatiotemporal correlation structure plays an essential role in spatiotemporal modeling.An anisotropic spatiotemporal correlation function which characterizes isotropic covariance functions in space and time, but anisotropic covariance in space-time provides an effective tool for modeling the spatiotemporal correlation structure.However, the assumption of isotropic spatial correlation in space implies that standardized TDP flows (PT3 distributed) at different sites have approximately the same values of the coefficient of skewness.For applications which exhibit significant spatial variation of the coefficient of skewness, preprocess of the data may be required.
(2) Stochastic simulation of a non-Gaussian random field

Fig. 1 .
Fig.1.Study area and flow stations.The region enclosed by the red-marked polygon represents the irrigation subgroup of an exemplar application for irrigation water shortage risk assessment (see details in section 6 and Fig.11).
the scale, shape and location parameters of the PT3 distribution.Using the method of L-moments (Hosking and Wallis 1995; Hosking and Wallis 1997) the three parameters were estimated as a c = 1.1059, b c = 0.7453, and f c = -0.8261.Combination of these parameter estimates ensures zero expectation and unit variance of the standardized TDP flows.

Fig. 2 .
Fig. 2. Long-term average TDP flows of 22 individual flow stations.Marked numbers represent flow station numbers in Fig. 1.
Fig. 3. Illustration of the spatiotemporal domain for standardized TDP flow simulation.(a) The spatiotemporal random field is spatially isotropic, but anisotropic in time-space domain.(b) A two-dimensional grid network is used for implementation of the column-preference spatiotemporal simulation.

Fig. 8 .
Fig. 8. (a) Map of standard deviations of the spatiotemporal correlation matrix of simulated TDP streamflows.(b) An enlarged area showing standard deviation maps of the correlation matrix of multisite streamflows lagged by zero to two TDPs.

Fig. 9 .
Fig. 9. Model-based spatiotemporal correlation coefficient of TDP streamflows with respect to lags in time and space.

Fig. 11 .
Fig. 11.Irrigation subgroup of an application for irrigation water shortage risk assessment (also see the red-marked polygon in Fig.1).

Fig. 12 .
Fig. 12.(a) TDP-specific average shortage ratios with and without implementation of the groundwater-withdraw mitigation measure.(b) And (c) ECDF of the 7-th TDP shortage ratios of subdivisions A and B, respectively.