Development and evaluation of a stochastic daily rainfall model with long-term variability

The primary objective of this study is to develop a stochastic rainfall generation model that can match not only the short resolution (daily) variability but also the longer resolution (monthly to multiyear) variability of observed rainfall. This study has developed a Markov chain (MC) model, which uses a two-state MC process with two parameters (wet-to-wet and dry-to-dry transition probabilities) to simulate rainfall occurrence and a gamma distribution with two parameters (mean and standard deviation of wet day rainfall) to simulate wet day rainfall depths. Starting with the traditional MC-gamma model with deterministic parameters, this study has developed and assessed four other variants of the MC-gamma model with different parameterisations. The key finding is that if the parameters of the gamma distribution are randomly sampled each year from fitted distributions rather than fixed parameters with time, the variability of rainfall depths at both short and longer temporal resolutions can be preserved, while the variability of wet periods (i.e. number of wet days and mean length of wet spell) can be preserved by decadally varied MC parameters. This is a straightforward enhancement to the traditional simplest MC model and is both objective and parsimonious.


Introduction
Observed rainfall data generally provide a single realisation of a short record, often not more than a few decades.The direct application of these data in hydrological and agricultural systems may not provide the necessary robustness in identification and implication of extreme climate conditions (e.g.droughts, floods).In particular, for urban water security analysis of reservoirs, long-term hydrologic records are required to sample extreme droughts that drive the security of the urban system (Mortazavi et al., 2013).However, the observed data may still be suitable to calibrate stochastic rainfall models that can, in turn, be used to generate long stochastic streamflow sequences for use in reservoir reliability modelling.In addition to historical and current scenarios, the stochastic models are useful to evaluate the climate and hydrological characteristics of future climate change scenarios (Glenis et al., 2015).
There is a major issue in the use of stochastic daily rainfall models.The daily models generally preserve the shortterm daily rainfall variability (since they are calibrated to the daily resolution data) but tend to underestimate the longerterm rainfall variability of monthly and multiyear resolutions (Wang and Nathan, 2007).Such underestimation is critically important for the application of these models in hydrological planning and design.Preserving the long-term variability is important for drought security analysis of reservoirs because the reservoir water levels usually vary at monthly to multiyear resolutions.The underestimation of longer-term A. F. M. K. Chowdhury et al.: Development and evaluation of a stochastic daily rainfall model variability of rainfall may cause an overestimation of reservoir reliability in urban water planning (Frost et al., 2007).Therefore, preserving key statistics of wet and dry spells, and rainfall depths in daily to multiyear resolutions, is important in stochastic rainfall simulation.
Markov chain (MC) models are very common for stochastic rainfall generation.A typical MC rainfall model is composed of two parts: a rainfall occurrence model that uses a transition probability between wet and dry days, and a rainfall magnitude model that uses a probability distribution of wet day rainfall depths (commonly a gamma distribution) fitted to the observed data.The two-part MC-gamma model is one of the most popular parametric models for daily rainfall simulation, primarily proposed by Richardson (1981) and known as WGEN (weather generator).In addition to rainfall, the WGEN also simulates temperature and solar radiation.While other models such as point process models (Cowpertwait et al., 1996) are also used for stochastic rainfall generation, this study has focused on MC-type models.
The first component of the MC model defines wet and dry days.This is determined by the state and order of the Markov process.Most MC models (Richardson, 1981;Dubrovský et al., 2004) use a simple two-state, first-order approach, that is, a day can be either "wet" or "dry" (two-state) and the state of the current day is only dependent on the state of the preceding day (first-order).Other models use higher states and orders; examples include the four-state model (Jothityangkoon et al., 2000), the alternating renewal process model with negative binomial distribution of wet and dry spell lengths (Wilby et al., 1998), the bivariate mixed distribution model (Li et al., 2013), and the multi-order model (Lennartsson et al., 2008).These models are more complex as the number of parameters required in the model increases with the number of states and orders of the Markov process.However, the two-state, first-order MC model can often reproduce the statistics of wet and dry periods just as well as these higher state/order models (Chen and Brissette, 2014).Dubrovský et al. (2004) recommended that, rather than trying an increased order MC, one should consider other approaches for better reproduction of wet and dry days.Mehrotra and Sharma (2007) proposed a modified MC process using memory of past wet periods, which has been found to reproduce the wet and dry spell statistics reasonably well.They also tested a first-order and a second-order process in their modified MC model and found that the second-order process provided only marginal improvements over the firstorder process.Another important finding of Dubrovský et al. (2004) was that the order of MCs generally had no effect on the variability of monthly rainfall depths.
The second component of the MC model is the probability distribution for the wet day rainfall.As the distribution of wet day rainfall is generally right-skewed (Hundecha et al., 2009), it is common practice to use rightskewed exponential-type distributions.Common distributions include the gamma distribution (Wang and Nathan, 2007;Chen et al., 2010), Weibull distribution (Sharda and Das, 2005), truncated normal distribution (Hundecha et al., 2009), and kernel-density estimation techniques (Harrold et al., 2003).A number of other studies fitted a mixture of two or more distributions; for example, the mixed exponential distribution (Wilks, 1999a;Liu et al., 2011), gamma and generalised Pareto distribution (Furrer and Katz, 2008), and transformed normal and generalised Pareto distribution (Lennartsson et al., 2008).However, the gamma distribution is the most commonly used distribution, because it has only two parameters, that can be calculated from the mean and standard deviation (SD) of wet day rainfall, and adequately represent the rainfall probability distribution functions.The parameterisation and application of the distribution in the model is straightforward.Although the gamma distribution has been found to be appropriate for simulating most of the variability in rainfall depth (Bellone et al., 2000), the major drawback of using a gamma distribution is that its tail is too light to capture heavy rainfall intensities (Vrac and Naveau, 2007).Therefore, direct use of a gamma distribution usually causes an underestimation of SD and extreme rainfall depths at monthly to multiyear resolutions.
A number of methods have been developed in an attempt to resolve the underestimation of long-term variability.The major approaches for resolving this issue include (i) models with mixed distributions, (ii) nesting-type models, (iii) models with rainfall-climate index correlation, and (iv) models with modified Markov chains.
The models with mixed distributions concentrate on the upper tail behaviour of the probability distribution of wet day rainfall depths.Since a single component distribution cannot incorporate the extreme rainfall depths well, a mixture of distributions is introduced.In these models, rainfall above a threshold depth is defined as "extreme" and two separate distributions are used to simulate the "extreme" and "small" rainfall amounts.Wilks (1999b) proposed a mixture of two exponential distributions with one shape parameter, but two scale parameters which are used to incorporate the extreme and small rainfall depths.In other models, the "extreme" rainfall depths are modelled by a generalised Pareto distribution (Vrac and Naveau, 2007) or a stretched exponential distribution (Wilson and Toumi, 2005), while small rainfall depths are modelled by a gamma distribution.Nonetheless, these models have difficulty in objectively defining the threshold corresponding to the "extreme value".Wilson and Toumi (2005) defined extreme rainfall as daily totals with an exceedance probability less than 5 %.Although Vrac and Naveau (2007) used a dynamic mixture to avoid choosing a threshold for "extreme", Furrer and Katz (2007) described the method as over-parameterised.Recently, Naveau et al. (2016) proposed a new model with a smooth transition between the "small rainfall" and "extreme rainfall" simulation process to generate low, moderate, and heavy rainfall depths without selecting a threshold.
Nesting models adjust the daily rainfall series at different temporal resolutions to obtain statistics that are optimal for all resolutions.These models initially generate a daily rainfall series, which is then modified to adjust the monthly and yearly statistics.Several models (Dubrovský et al., 2004;Wang and Nathan, 2007;Srikanthan and Pegram, 2009;Chen et al., 2010) use the nesting method.They generally generate a daily rainfall series, then the generated daily rainfall data are aggregated to monthly rainfall values, and these monthly values are modified by a lag-1 autoregressive monthly rainfall model.The modified monthly rainfall values are aggregated to annual rainfall values and these values are then modified by another lag-1 autoregressive annual model (Srikanthan and Pegram, 2009).The nesting-type models generally perform well to reproduce the rainfall variability at all resolutions.Dubrovský et al. (2004) also showed satisfactory performance of their nesting-type model to reproduce the variability of monthly streamflow characteristics and the frequency of extreme streamflow.Although the nesting-type models preserve the daily, monthly and yearly statistics, they are generally based on subjective statistical adjustments and thus have a weak physical basis.
Some parametric models introduced the influence of largescale climate mechanisms such as the El Niño/Southern Oscillation (ENSO) in parameterisation (Hansen and Mavromatis, 2001;Furrer and Katz, 2007).Bardossy and Plate (1992) used the correlation between atmospheric circulation patterns and rainfall in a transformed conditional multivariate autoregressive AR(1) model for daily rainfall simulation.Katz and Parlange (1993) developed a model with parameters conditioned on the ENSO indices.Yunus et al. (2016) developed a generalised linear model for daily rainfall by using ENSO indices as predictors.Although the climate indices were often not strongly correlated to the rainfall, Katz and Zheng (1999) described it as a diagnostic element to detect a "hidden" (i.e.unobserved) index, which could be used to obtain long-term variability.Thyer and Kuczera (2000) developed a hidden state MC model for annual data, while Ramesh and Onof (2014) developed a hidden state MC model for daily data.The major drawback of this model approach is that the hidden index is unobserved and its origin is unknown.
Modified MC models concentrate not only on the order of MC but also introduce modifications to the parameterisation of the MC process to better reproduce the rainfall variability.The transition probabilities are generally modified by considering their long-term variability (i.e.memory of past wet and dry periods), and the wet day rainfall depth is modelled using a nonparametric kernel-density simulator conditional on previous day rainfall (Lall et al., 1996;Harrold et al., 2003).The nonparametric kernel-density techniques usually used resampling of observed data (Rajagopalan and Lall, 1999).While these models perform reasonably well, they usually cannot generate extreme values higher than the observed extremes, because only the original observations are resam-pled in the model (Sharif and Burn, 2006).Mehrotra and Sharma (2007) proposed a semi-parametric Markov model, which was further evaluated by Mehrotra et al. (2015).To incorporate the long-term variability, they modified the transition probabilities of the MC process by taking the memory of past wet periods (i.e.beyond lag-1) into account, while the wet day rainfall depths were simulated by a nonparametric kernel-density process.For rain gauge data around Sydney, the semi-parametric model preserved the rainfall variability at daily to multiyear resolutions (Mehrotra et al., 2015).
The MC models that focus specifically on resolving the underestimation of long-term variability involve subjective assumptions and limitations.In the models with mixed distributions, defining a certain rainfall depth as an extreme value is subjective.The nesting-type models used empirical adjustment factors, generally without physical foundation.The hidden indices of hidden state MC models are unobserved.The models with modified MC parameters, modified the transition probabilities of wet and dry periods to obtain long-term variability, but used the kernel-density technique to resample wet day rainfall depths from observed records.Therefore, they usually cannot generate more extreme values than the observed extremes.
The overarching objective of the research, that this paper forms part of, is to develop a stochastic rainfall generator that can be calibrated to daily rainfall data derived from dynamically downscaled global climate simulations, and which also preserves long-term variability (Evans et al., 2014).A common problem with these simulations is that typical computational CPU limits mean that the length of the simulation is rarely more than a few decades, not long enough to facilitate stochastic assessment of the reliability of water supply reservoirs (e.g.Lockart et al., 2016).Accordingly, we need a rainfall simulator that can be calibrated and run at the daily timescale (to be used as input into a hydrology model at the daily resolution), but which has the right statistical properties (specifically variability about the mean) when averaged over periods up to a decade.In this paper, we develop and test five models using observed rainfall at two sites in Australia with contrasting climates.
Accordingly, this study details the development of a MC model for stochastic generation of daily rainfall.This MC model uses a two-state MC process with two parameters (wet-to-wet and dry-to-dry transition probabilities) for simulating rainfall occurrence and a two-parameter gamma distribution (mean and SD of wet day rainfall) for simulating wet day rainfall depths.Five variants of the MC model, with gradually increasing complexity of parameterisation, are developed and assessed.Starting with a very simple model against which the performances of the other models will be compared.Each of the successive models provide better performance in reproducing the variability and dependence of observed rainfall over the range of resolutions from day to decade, and we assess the incremental improvements in per- formance against the incremental increases in model complexity.

Data and study sites
This study has used daily rain gauge data from Sydney Observatory Hill and Adelaide Airport stations (station number 66062 and 023034, respectively) obtained from the Bureau of Meteorology (BoM), Australia (Fig. 1) for 1979-2008(BoM, 2013)).These two stations have been used because they provide a contrast between a highly seasonal Mediterranean climate with low interdecadal variability in Adelaide and a relatively non-seasonal climate with high interdecadal variability in Sydney (see Fig. 2).Risbey et al. (2009) also showed that the major climate drivers of rainfall (e.g.ENSO) in Sydney and Adelaide are different for all seasons.This paper also used the Oceanic Nino Index (ONI) and the Interdecadal Pacific Oscillation (IPO) index at monthly resolution for the 1979-2008 period (Folland, 2008;NOAA, 2014).These climate indices are used to develop two variants of the MC models discussed in Sect.4.2.2.
3 Model assessment procedures

Statistics for assessment of model performance
Each model developed in this study has been assessed to understand its ability to reproduce the distribution and autocorrelation of observed rainfall.Assessment of the distribution and autocorrelation are generally used to inform the suitability of the model for urban drought security assessment.The assessment criteria of each model consider its ability to re-produce (i) mean, SD, and 95th percentile of rainfall depths at daily to multiyear resolutions; (ii) mean and SD of the number of wet days and mean length of wet spells at monthly to multiyear resolutions; and (iii) month-to-month autocorrelations of monthly rainfall depths and monthly number of wet days.The performances of the MC models for dry period statistics are found to be similar to the wet period statistics (the term "wet period statistics" will hereafter refer to the number of wet days and mean length of wet spells), and hence, only representative results for annual mean length of dry spells are shown.
At daily and monthly resolutions, the distribution statistics are assessed for each month starting from January; while at multiyear resolutions, the distribution statistics are assessed for 1 to 10 overlapping years.Mean length (in days) of wet spells are calculated at monthly, and annual resolution by extracting wet spells of one or more consecutive wet days (two successive wet spells are separated by at-least one dry day) and using Eq. ( 1): mean length of wet spell = (length of wet spells) (occurrences of wet spells) . (1) Similar to wet spells, the mean length of dry spells are also calculated at monthly and annual resolution by extracting dry spells of one or more consecutive dry days.

Calculation of Z scores
For the distribution statistics (i.e.mean and SD) of rainfall depths and wet periods (number of wet days and mean length of wet spells), this study has calculated the expected value and error limit (SD) to calculate the Z score of a model simulation.The calculation of the Z score is as follows: 1. Run the model 1000 times using the probability distribution of the parameters calibrated to the observed data, with each run being the same length as the observed data.
2. Calculate the desired statistics (e.g.mean and SD of the daily rainfall depths) in each run, which gives 1000 realisations of each statistic.
3. For each statistic, calculate the mean (expected value) and SD (error limit) of the 1000 realisations.
4. Calculate the Z score of each statistic by comparing the expected value with the respective observed value (calculated from the observed data), as follows: A Z score between −2 and +2 for a statistic indicates that the observed value falls within the 95 % confidence limits of the simulated rainfall, assuming a normal distribution approximates the sampling distribution of Z.A Z score less than −2 or greater than +2 suggests that the statistic is overor under-estimated, respectively, in the model simulation.

Markov chain (MC) models
This study has developed and assessed the following five variants of a Markov chain (MC) model: -

Model 1: Average Parameter Markov Chain (APMC) model
The first MC model -the APMC -is a traditional two-part MC-gamma distribution model.This is similar to the rainfall generator proposed by Richardson (1981), widely known as the WGEN model.The exception is that the parameters in WGEN were smoothed with Fourier harmonics, which has not been done in the case of APMC parameters.Although APMC is not the final model of this study, it is the baseline modelling approach against which the more sophisticated models developed in this study are compared.
The APMC simulates the daily rainfall in two steps: daily rainfall occurrence (i.e.wet and dry day) simulation by firstorder Markov chain, and wet day rainfall depth simulation by gamma distribution.To incorporate the seasonal variability in the model, the APMC uses a separate set of parameters for each month, where the first month of the simulation is January.

Rainfall occurrence simulation
The APMC uses 24 (2 parameters × 12 months) MC parameters, and transition probabilities of dry-to-dry day (P 00 ) and wet-to-wet day (P 11 ) for dry and wet day occurrence simulation.In addition, the unconditional probability of a dry day (π 0 ) in January is used to simulate rainfall occurrence for the first day of the series.In the model calibration, these deterministic MC parameters are calculated from the observed daily rainfall data.To calculate these parameters, a day with rainfall depth of 0.3 mm and above has been considered a wet day, otherwise it was considered a dry day (similar to Mehrotra et al., 2015).In simulation, the MC parameters are used in a Monte-Carlo process to simulate the occurrences of wet and dry days.

Rainfall depth simulation
After simulation of the rainfall occurrence using MC parameters, the next step is to generate rainfall depths for the wet days.The rainfall depth for dry days is zero.The APMC rainfall depth simulation process assumes that (i) daily rainfall depth for wet days follows a gamma distribution, and (ii) the rainfall depth for a wet day is independent of the rainfall depth of the preceding day.
The gamma distribution has two parameters, α (shape parameter) and β (scale parameter), with mean µ = αβ and variance σ 2 = αβ 2 .Since both α i and β i are directly proportional to and can be derived from µ i and σ i of wet day rainfall of the month i, then during calibration of the model it is convenient to calculate µ i and σ i values from the daily rainfall observed data.The appropriate ratios of µ i and σ i can then be used in the rainfall depth generation process using the gamma distribution.Therefore, µ i and σ i will be referred to as the gamma distribution parameters in further discussions of this paper.
In the calibration of APMC, deterministic average values of µ i and σ i are calculated from the entire period of the data record for each month.This gives 12 values of µ and σ each.In simulations, the rainfall depth for each wet day of a month i is generated using the µ i and σ i values of the respective month using the gamma distribution.In generating the rainfall depth for a wet day, if a random sample from the gamma distribution gives a rainfall depth less than 0.3 mm then the rainfall for that day is set to 0.3 mm (i.e. the threshold rainfall depth), while the rainfall depths for dry days are set to 0.0 mm.Chowdhury (2017) showed that setting rainfall below 0.3 mm to 0.3 mm for the lowest rainfall depth does not significantly affect the overall distribution of modelled rainfall depths.

Independence of rainfall depths in successive wet days
The APMC assumes that the rainfall depth for a particular day is independent of the rainfall depth of the preceding day.
To validate this assumption, this study examined the autocorrelation of wet day rainfall depths, and only found very weak lag-1 autocorrelations (r 2 < 0.1) for both Sydney and Adelaide.This finding is consistent irrespective of seasonal variations.The conclusion is that the underlying assumption of daily independence of the APMC is consistent with the respective characteristic of the observed data.

Model 2: Decadal Parameter Markov Chain (DPMC) model
Section 6 will show that the APMC significantly underestimates the rainfall variability at monthly to multiyear res- olutions.The DPMC assumes that the interannual rainfall variability can be captured by the decade-to-decade variability of the parameters that APMC failed to capture.The idea is to divide the observed rainfall sample into subsamples of 10 years duration (similar models with climate-based subsamples are discussed in Sect.4.2.2).For example, a 30-year rainfall sample is divided into three subsamples of 10 years in duration.Then, 4 × 12 parameters of P 00 , P 11 , µ, and σ (one set of four parameters for each of the 12 months) are calculated from each of the subsamples.The simulation proceeds in a way similar to the APMC, except that the deterministic, decadal average, parameters of DPMC are varied from decade to decade.

Decadal variability of DPMC parameters
Figure 2 shows the DPMC values of P 11 and µ for each decade along with APMC values (i.e. the 30-year averages) for Sydney and Adelaide.For Sydney, DPMC values of P 11 and µ show clear variability between the three decadal samples and deviations from the APMC values.However, DPMC values of P 11 and µ for Adelaide show less variability between the decadal samples.The use of decadally varied parameters in DPMC is subject to the question of how significant the decadal variability of these parameters is -is the decadal variability statistically significant or just sampling variability?Therefore, the statistical significance of the decadal variability of DPMC parameters were examined by Monte-Carlo simulations as per Sect.3.2.This examination suggested that the sampling variability of DPMC parameters in decadal samples is mostly within the sampling variability of their corresponding APMC values (not shown).This suggests that the decadal variability of DPMC parameters is not statistically significant.

Potential impact of climate modes
This study has also investigated other subsampling approaches of the MC-gamma parameters similar to the DPMC.In these models, this study has calibrated the MCgamma parameters to subsamples of rainfall time series divided according to the phases of IPO (e.g.positive and negative) and ENSO (La Niña, Neutral, and El Niño).Since previous studies (Verdon-Kidd et al., 2004) found that the interannual variabilities of east-Australian rainfall are influenced by these large-scale climate drivers, the idea behind these models was to introduce more interannual variability to the model by simulating rainfall for different phases of climate drivers with parameters calibrated to respective phases.These climate-based models are very similar to DPMC, except that the subsamples are different.The following two types of climate-based models have been tested: -The IPO based model: the observed data for every month was divided into two subsamples according to the positive and negative phases of the monthly IPO index (e.g. for January, data of the years with positive IPO index and data of the years with negative IPO index are separated).Then, for each month, the MC-gamma parameters (P 00 , P 11 , µ, and σ ) are calibrated to each subsample.In simulation, the rainfall for the months of each IPO phase were modelled by using parameters of the respective phase.
-The ENSO based model: the observed data for every month was divided into three subsamples according to monthly ONI index: La Niña (ONI ≤ −0.5), neutral (−0.5 < ONI < 0.5), and El Niño (ONI ≥ 0.5).Then, the MC-gamma parameters are calibrated to each subsample and the rainfall for the months of each ENSO phase were modelled by using parameters of the respective phase.

Model 3: Compound Distribution Markov Chain (CDMC) model
The results in Sect.6 will show that neither APMC nor DPMC can satisfactorily reproduce the SD of rainfall depths for monthly to multiyear resolutions.Therefore, in the third MC model -the CDMC -this study has incorporated the long-term variability of rainfall depths by introducing random variability in µ and σ .However, for wet and dry period simulation, the CDMC still uses the deterministic parameters of P 00 and P 11 , as in the APMC.Thus, this model stochastically varies the rainfall depth model, but not the rainfall occurrence model.In the CDMC, µ i and σ i are randomly sampled for each month of each year.The random sampling was done independently of the sampling for the preceding month(s).To estimate the distribution of µ i and σ i , this study has calculated µ i and σ i for every month of every year from the observed data.For example, from the 30-year observed data, for January (i = 1), this study has calculated 30 samples of µ 1 and σ 1 values each.
By testing the probability distributions of µ i and σ i values for each month, this study has found that both µ i and σ i values for each month are lognormally distributed.Figure 3 shows the lognormal probability plots of µ i and σ i values for July (i = 7), which is representative of the other months.The r 2 for log µ i and log σ i are generally above 0.90, indicating a very good fit of the lognormal distributions.Additionally, the hypothesis that log µ i and log σ i are normally distributed is supported by the Kolmogorov-Smirnov test at 5 % significance level.In addition to the lognormally distributed µ i and σ i values, this study has also found that the log µ i and log σ i values for each month are strongly correlated with each other with correlation coefficient r c,i around 0.90 (Fig. 4).Therefore, for each month i, this study has fitted a bivariate-normal distribution to the log µ i and log σ i values with parameters (λ µ,i , ζ µ,i ), (λ σ,i , ζ σ,i ), and r c,i .The λ and ζ parameters denote the mean and SD of the log variate, while r c is the correlation coefficient between log µ and log σ .
At the start of each month of each year of the simulation, the log µ i is sampled from its fitted normal distribution log µ i ∼ N (λ µ i ζ 2 µ i ) for month i.Then, the log σ i is sampled from the fitted conditional normal distribution: (3) These stochastically sampled µ i and σ i values are then used to generate rainfall in the wet days for the month in question, while the sequence of wet and dry days is determined using the deterministic APMC values of P 00,i and P 11,i .However, the sampled µ i and σ i values of a month (i) are not correlated to the µ i−1 and σ i−1 of the preceding month (i − 1) as this study has found that the month-to-month autocorrelations of µ and σ values are not strong (Fig. 5).

Model 4: Hierarchical Markov Chain (HMC) model
The results in Sect.6 will show that the CDMC cannot satisfactorily reproduce the SD of wet periods for monthly to multiyear resolutions.Therefore, in the fourth MC modelthe HMC -we introduce stochastic variation in both MC and the gamma distribution models to incorporate long-term variability of rainfall depths as well as wet and dry periods.In the calibration, for month i, the P 00,i and P 11,i are calculated for each month of each year from the observed data.For month i, these P 00,i and P 11,i values (e.g. 30 P 11,7 values for July from the 30-year observed data) are found to be normally distributed with values between 0 and 1 (Fig. 6).Therefore, this study has fitted a truncated normal distribution, bounded by 0 and 1, to the calculated P 00 and P 11 values.In simulation, for each year, the P 00,i and P 11,i are sampled from their truncated normal distributions.This procedure is similar to what was done for µ i and σ i to sample from bivariate-lognormal distribution.However, it does not include a bivariate distribution because the correlation between P 00,i and P 11,i was weak.

Impact of autocorrelations on stochasticity of MC parameters
In the HMC, the sampled MC parameters of each month are independent of the parameters of the preceding month.How-ever, this study has found strong month-to-month autocorrelations of the P 00 and P 11 for Adelaide (Fig. 5a), although the autocorrelations are weak for Sydney (Fig. 5b).Therefore, this study has tested an alternative to the HMC (referred to as HMC2), which uses a lag-1 autocorrelation equation (a similar equation was used by Wang and Nathan (2007) in their rainfall depth model) in the stochastic sampling of P 00,i and P 11,i from the truncated normal distribution.The following lag-1 autocorrelation equation has been used to modify the randomly sampled P 00,i (same method used for P 11,i ) for month i by correlating with the P 00,i−1 of month i − 1 (preceding month): P 00,i − mean(P 00,i ) sd(P 00,i ) = r × P 00,i−1 − mean(P 00,i−1 ) sd(P 00,i−1 ) + (1 − r 2 ) 1/2 P 00,i − mean(P 00,i ) sd(P 00,i ) , where for a month i (e.g.January), mean(P 00,i ) is the mean of the yearly varied parameter values calculated from observed data for month i (e.g.mean of 30 P 00,7 values for July), sd(P 00,i ) is the SD of the yearly varied parameter values calculated from observed data for month i, mean(P 00,i−1 ) is the mean of the parameter values calculated from observed data for month i − 1, sd(P 00,i−1 ) is the SD of the parameter values calculated from observed data for month i − 1, r is the lag-1 autocorrelation coefficient for observed month-to-month autocorrelation of P 00 (constant for all month), -P 00,i is the stochastic parameter value sampled from a truncated normal distribution fitted to the yearly varied observed parameter values for month i, -P 00,i−1 is auto-correlated parameter value for month i − 1 (used to simulate the dry days of the preceding month), -P 00,i is the final auto-correlated parameter value which was used in simulation of dry days for month i.
P 11,i for month i is sampled using a similar process.

Impact of cross-correlations on stochasticity of MC parameters
We observed a strong positive correlation between P 11,i , and log µ i and log σ i , although the correlations between P 00,i , and log µ i and log σ i are weak.Therefore, another alternative to HMC (referred to as HMC3) was tested by using a multivariate sampling system for the P 11,i , µ i , and σ i , while P 00,i remains independent.

Model 5: Decadal and Hierarchical Markov Chain (DHMC) model
Section 6 will show that the CDMC, with APMC values of MC parameters, significantly underestimates the wet period variability at multiyear resolutions, while the HMC (including the two alternatives HMC2 and HMC3) with stochastic MC parameters, significantly overestimates the wet period variability at monthly resolution.However, we found that the DPMC can satisfactorily preserve the wet period variability at both monthly and multiyear resolutions, although it underestimates the rainfall depth variability.Therefore, in the DHMC model, this study has used the DPMC values of MC parameters (the parameter values vary for each decade of simulation) for simulation of wet and dry days, while the stochastic parameters of the gamma distribution (same as CDMC) are used for simulation of wet day rainfall depths.

Methodological comparison of five MC models
The following points discuss the key common features in the five MC models of this study, while other key methodological comparisons are shown in Table 1: -All models use first-order MC parameters to simulate the rainfall occurrences and gamma distribution to simulate rainfall depths in wet days.
-Simulation of rainfall depth for each wet day is independent of the rainfall depth of the preceding day.
-Separate sets of parameters are used for each month (e.g. 12 sets of MC and gamma parameters) to incorporate seasonal variability.

Model comparison for distribution statistics
This section compares the performances of the five MC models for the mean and SD of rainfall depths and wet period statistics (i.e.number of wet days and mean length of wet spells).

Mean and SD of rainfall depths
Figures 7, 8, and 9 compare the five MC models for the mean and SD of rainfall depths at daily, monthly, and multiyear resolutions, respectively.Figure 9 also compares the 95th percentile of multiyear rainfall depths.For mean and SD of rainfall depths, the performances of APMC and DPMC are similar.The performances of CDMC, HMC, and DHMC are also similar, but different from APMC and DPMC.All five models preserve the mean (i.e.satisfactorily reproduce the observed mean) rainfall depths at all resolutions with Z scores between −2 and +2.However, the CDMC, HMC, and DHMC show a tendency to underestimate the mean with mostly positive Z scores (between 0 and +2).The APMC and DPMC preserve the SD of rainfall depths only at daily resolution and significantly underestimate the SD at monthly and multiyear resolutions for Sydney but preserve the SDs at all resolutions for Adelaide (Figs. 7, 8, and 9).The CDMC, HMC, and DHMC preserve the SD of rainfall depths at all resolutions for both stations except a slight tendency to underestimate the SD for February and November at daily resolution in Sydney.We conclude that those models with stochastic parameters for the gamma distribution (i.e.CDMC, HMC, and DHMC) best preserve SDs at all resolutions for both stations.For the 95th percentile of rainfall depths, we found that models which can preserve the SD at a given resolution can also preserve the 95th percentile at that resolution and vice versa.In Fig. 9, the representative results at multiyear resolution (average of the absolute values of Z scores for daily and monthly resolutions are shown in Table 2) show that the CDMC, HMC, and DHMC preserve the 95th percentile for both stations but the APMC and DPMC underestimate the statistic for Sydney.

Mean and SD of number of wet days
Figures 10 and 11 compare the five MC models for the mean and SD of number of wet days at monthly and multiyear resolutions, respectively.All five models preserve the mean of number of wet days for both monthly and multiyear resolutions.For the SD of the monthly number of wet days, all models except HMC can satisfactorily reproduce the SD with Z scores between −2 and +2 for almost all months of both stations, while the HMC tends to overestimate the SD (Fig. 10).For the SD of multiyear number of wet days, the APMC and CDMC significantly underestimate the SD for Sydney but preserve the statistic for Adelaide.The DPMC and DHMC perform similarly and satisfactorily to preserve the SD of multiyear number of wet days for both Sydney and Adelaide, while HMC also preserves the statistic for both stations.We conclude that the models with decadally varied MC parameters (i.e.DPMC and DHMC) perform satisfactorily at reproducing the variability of the number of wet days at monthly and multiyear resolutions for both stations.

Mean and SD of mean length of wet and dry spells
Figure 12 compares the five MC models for the mean and SD of mean length of wet and dry spells at annual resolution.The averages of the absolute values of the Z scores for monthly resolution are shown in Table 2.The comparative performances of the five MC models for the mean and SD of mean length of wet spells at monthly (Table 2) and annual (Fig. 12) resolutions are mostly consistent with their respective performances for mean and SD of number of wet days.All models except HMC preserve the mean and SD of mean length of wet spells, while the HMC tends to overestimate the SD (Fig. 12).
For the mean and SD of mean length of dry spells, we found that models that can preserve the wet spells distributions also preserve the dry spells distributions and vice versa.As a representative result, the Z scores for the mean and SD of annual mean length of dry spells shown in Fig. 12 indicate that all models except HMC preserve both mean and SD, while HMC overestimates the SD. Figure 12   cates that the DPMC and DHMC perform better (Z scores closer to zero) than the APMC and CDMC to reproduce the SD of annual mean length of dry spells.
We conclude that models with decadally varied MC parameters (i.e.DPMC, DHMC) perform relatively better and more satisfactorily at reproducing the variability of the length of wet and dry spells.The HMC introduces too much variability in the length of wet and dry spells, while the APMC and CDMC tend to underestimate the variability.

Potential impact of climate modes
Since the DPMC significantly underestimates the SD of rainfall depths at monthly and multiyear resolutions, the major target of the models with subsamples according to cli-mate modes such as IPO and ENSO indices (discussed in Sect.4.2.2) was to preserve the SD of rainfall depths at monthly and multiyear resolutions.However, we found that these climate-based models also significantly underestimate the SD of rainfall depths at month and multiyear resolutions with performances similar to the DPMC, and are therefore not considered further.

Impact of stochasticity on MC parameters
Since the HMC significantly overestimates the SD of monthly wet periods (i.e.number of wet days and mean length of wet spell), the major target of the HMC2 and HMC3 models (with a lag-1 autocorrelation equation and a multivariate sampling system, respectively; see Sect.  was to better preserve the SD.However, these models also significantly overestimate the SD of monthly wet periods with performances similar to the HMC (negative Z scores less than −2 for all months).We conclude that the models with stochastic, yearly varied, parameters for the MC part of the model (i.e.HMC, HMC2, and HMC3) consistently overestimate the variability of monthly wet periods.

Overall performances
Table 2 shows the average of the absolute values of Z scores (average of 12 values at daily and monthly resolutions and 10 values at multiyear resolution) for the distribution statistics of rainfall depths, and wet and dry periods at daily, monthly, annual, and multiyear resolutions.It shows that models 1-4 (APMC, DPMC, CDMC, and HMC) fail to preserve the following statistics: the APMC fails to preserve the SD and 95th percentile of rainfall depths and SD of number of wet days at multiyear resolution for Sydney, the DPMC fails to preserve the SD and 95th percentile of rainfall depths at multiyear resolution for Sydney, the CDMC fails to preserve the SD of number of wet days at multiyear resolution for Sydney, the HMC fails to preserve SD of mean length of wet spell at annual resolution for both Sydney and Adelaide.
However, model 5, DHMC, has preserved all of the statistics for both stations.We conclude that the DHMC is better than the other four models at reproducing the distributions of rainfall depths, and wet and dry periods for resolutions varying from daily to multiyear.

Reproduction of seasonal autocorrelations
Figure 13 compares how the five MC models reproduce the month-to-month autocorrelations of the monthly number of wet days and monthly rainfall depths.For Adelaide (Fig. 13a), the lag-1 and lag-12 autocorrelations are strong with systematic seasonal variation, which have been reproduced very well in the corresponding APMC, DPMC, CDMC, and DHMC simulations, while the HMC (the model with stochastic MC parameters) tends to underestimate the autocorrelations.For Sydney (Fig. 13b), the month-to-month autocorrelations of monthly number of wet days and monthly rainfall depths in the data are weak and all models perform well.

Discussion
The primary motivation of this study was to develop a stochastic rainfall generation model that can reproduce not only the short resolution (daily) variability, but also the longer resolution (monthly to multiyear) variability of observed rainfall.Preserving long-term variability in rainfall models has been a difficult challenge for which a number of solutions have been proposed in the stochastic rainfall generation literature.The solutions developed and tested by this study are relatively simple MC models: two MC parameters (P 00 and P 11 ) of two-state, first-order processes defining the wet and dry days, and two gamma-distribution parameters (µ and σ ) defining the rainfall depths in wet days.For seasonal variability, the models operate at daily time step with monthly varying parameters for each of the 12 months.Starting with the simplest MC-gamma modelling approach with deterministic parameters (similar to Richardson, 1981), this study has developed and assessed four other variants of the MC-gamma modelling approach with different parameterisations.The key finding is that if the parameters of the gamma distribution are randomly sampled from fitted distributions prior to simulating the rainfall for each year, the variability of rainfall depths at longer resolutions can be preserved, while the variability of wet periods (i.e.number of wet days and mean length of wet spell) can be preserved by decadally varying parameters for the MC model.This is a straightforward enhancement to the traditional simplest MC model, and the enhancement is both objective and parsimonious.
The overall comparative performances of the models to reproduce the distribution and autocorrelation characteristics of observed rainfall are as follows: -For the simulation of the distribution of rainfall depths, the performances of the APMC and DPMC with deterministic gamma parameters are similar, although DPMC with more parameters (e.g. the decadally varying MC parameters) performs slightly better.The performances of CDMC, HMC, and DHMC are similar as they use the same stochastic sampling for the parameters of the gamma distribution.
-For the mean and SD of daily rainfall depths, all five models perform satisfactorily.Good reproduction of daily statistics is expected as the model parameters are calibrated to daily time series.While the APMC and DPMC reproduce the statistics almost exactly, the CDMC, HMC, and DHMC show a slight tendency to underestimate the SD.This indicates that the stochastic parameters of these three models slightly affected their performances at daily resolution compared to the APMC and DPMC with deterministic parameters.
-For the monthly to multiyear resolution, the APMC and DPMC reproduce the mean of rainfall depths well, but significantly underestimate the SD of rainfall depths.The underestimation of rainfall variability at monthly to multiyear resolutions by APMC-like models with deterministic parameters is a well-known limitation of these models (Wang and Nathan, 2007).Although the DPMC uses more parameters than the APMC, the DPMC has not significantly improved performance in reproducing the SD of rainfall depths at monthly to multiyear resolutions.Other models similar to DPMC (e.g.models with parameters varying for phases of IPO or ENSO) show similar performances to the DPMC and still systematically underestimate the SD of rainfall depths at monthly to multiyear resolutions.This suggests that the use of deterministic parameters in DPMC-like models might not be adequate to satisfactorily reproduce the SD of rainfall depths at longer resolutions.
-While the APMC and DPMC, with deterministic parameters for the gamma distribution, significantly un-Hydrol.Earth Syst.Sci., 21, 6541-6558, 2017 www.hydrol-earth-syst-sci.net/21/6541/2017/ derestimate the SD of rainfall depths at monthly to multiyear resolutions, the CDMC, HMC, and DHMC, with stochastic parameters for the gamma distribution, preserve the SD of rainfall depths at monthly to multiyear resolutions.This indicates that the stochastic parameters for the gamma distribution are useful to incorporate the longer-term variability of rainfall depths.However, these three models show a tendency to underestimate the mean of rainfall depths at all resolutions.
-The models that can preserve the SD of rainfall depths can also preserve the 95th percentile of rainfall depths.
-For the simulation of the distribution of wet periods, the performances of the APMC and CDMC are similar as both models use the same deterministic MC parameters.With a similar trend, the DPMC and DHMC perform better than the APMC and CDMC, while DPMC and DHMC use more deterministic MC parameters.The performance of the HMC, with stochastic MC parameters, is different (discussed below) from the other four models (that use deterministic MC parameters).
-For the mean of wet period statistics (e.g.number of wet days and mean length of wet spells) at monthly to multiyear resolutions, all models except HMC perform similarly and satisfactorily, while the HMC tends to overestimate the mean.We conclude that introducing stochasticity from year to year into the MC parameters, as in HMC, degrades the performance.
-For the SD of monthly wet period statistics, all models except HMC perform similarly and satisfactorily, while the HMC significantly overestimates the SD.This indicates that the stochastic MC parameters of the HMC introduce excessive variability in the wet period simulation at monthly resolution.This study has further examined two other variants of the HMC with different stochastic parameterisation of the MC process, but they did not perform better than the HMC.We conclude that introducing stochasticity from year to year into the MC parameters, as in HMC, degrades the ability to reproduce the variability about the mean of all of the wet period statistics.
-For the SD of wet period statistics at annual and multiyear resolutions, the APMC and CDMC tend to underestimate the SDs.This suggests that the APMC values of MC parameters (same monthly parameter values for each year of simulation) limits the reproduction of the wet period variability at multiyear resolutions.However, the APMC and CDMC preserved the multiyear SDs for Adelaide, where the interdecadal variability of MC parameters is less variable.This suggests that for locations with less variability of wet-to-wet and dry-todry day transitions, a single set of deterministic MC parameters is adequate, however for locations with more transition variability, a single set of MC parameters (i.e.not varying with time) is insufficient, as it cannot introduce enough variability.
-The DPMC and DHMC with decadally varied MC parameters show a better ability to reproduce the SD of annual mean length of wet spells and SD of multiyear number of wet days.This suggests that the decadally varied MC parameters can significantly improve the simulation of wet period variability, although the decadally varied gamma parameters cannot improve the simulation of rainfall depth variability.However, the HMC preserves the SD of multiyear number of wet days but overestimates the SD of annual mean length of wet spells.This suggests that the monthly and annually varying stochastic MC parameters can improve the simulation of wet period (i.e.number of wet days and mean length of wet spell) variability at multiyear resolutions, although they significantly overestimate the wet period variability at monthly and annual resolutions (i.e. they introduce too much variability).
-The models that can preserve the wet spells distributions can also preserve the dry spells distributions and vice versa probably because the wet and dry days are modelled using similar transition probabilities of wet-to-wet and dry-to-dry days, respectively.
-The autocorrelation assessments have shown that the APMC, DPMC, CDMC, and DHMC can satisfactorily reproduce the strong lag-1 and lag-12 monthly autocorrelations of monthly number of wet days and monthly rainfall depths.However, the HMC (the only model with monthly and annually varying MC parameter values) tends to underestimate the autocorrelations, which is possibly due to excessive variability in wet period simulation at monthly resolution.

Conclusions
Each model developed in this study has advantages and disadvantages.The APMC and DPMC with deterministic parameters significantly underestimate the variability of rainfall depths at monthly to multiyear resolutions.This systematic underestimation of the rainfall depth variability at monthly to multiyear resolutions is critical for using the models in urban water security assessment as the reservoir water levels usually vary at these longer resolutions.The CDMC, HMC, and DHMC with stochastic parameters of the gamma distribution preserve the rainfall depth variability at all resolutions, but the CDMC and HMC have limitations in reproducing the variability of wet periods.The CDMC with APMC values of MC parameters tends to underestimate the multiyear variability of wet periods, while the HMC with stochastic MC parameters tends to overestimate the monthly variability of wet periods.However, the DHMC with decadally varied MC parameters (same as DPMC) performs better than the CDMC and HMC, and preserves the wet period and dry period variability at monthly to multiyear resolutions.Among the five MC models of this study, the overall performance of the DHMC is the best.The DHMC model has (1) monthly varying MC parameters that vary from decade to decade and (2) stochastic parameters for the gamma rainfall distribution, where the parameters are randomly varied from year to year using a probability distribution function that is derived for each month of the year.While the DHMC has great potential to be used in hydrological and agricultural impact studies (e.g.urban drought security assessment), there are two important limitations of the DHMC: -The DHMC tends to underestimate the mean of multiyear rainfall depths, which is probably linked to the use of stochastic gamma parameters.A more sophisticated stochastic sampling strategy for the gamma parameters might overcome this limitation.
-The performance of the DHMC suggests that the use of decadally varied MC parameters are effective to incorporate the long-term variability of wet periods (although the use of decadally varied gamma parameters in DPMC was not effective to incorporate the long-term variability of rainfall depths).However, other climatebased subsamples (e.g. according to the ENSO phases) instead of decadal samples can be used for parameter calibration.This study tested the subsamples according to the phases of IPO and ENSO climate modes with a focus on incorporating the long-term variability of rainfall depths, but has not incorporated climate-based sub sampling into DHMC because DHMC had not been developed at the time this analysis was performed.A more comprehensive assessment of such ideas might improve the wet period simulation of the DHMC.
In a subsequent paper, the performances of the CDMC, HMC, and DHMC will be compared against the semiparametric model of Mehrotra and Sharma (2007) using rain gauge data from 30 stations around Sydney (those used in Mehrotra et al., 2015) and the 12 stations (Fig. 1) around Australia.

Figure 1 .
Figure 1.Location map of 12 rain gauge stations around Australia.This study has presented the assessment results of the developed models for Sydney and Adelaide stations (red circled) only.The shaded green, yellow, and red colours indicate the coastal, inland, and monsoonal areas, respectively.

Figure 2 .
Figure 2. Comparison of the decadal variability of the DPMC parameters P 11 and µ (mm) with the APMC parameters.

Figure 3 .
Figure 3. Lognormal probability plots of µ and σ for July (typical of other months).

Figure 4 .
Figure 4. Correlation between log µ and log σ for July (typical of other months).

Figure 6 .
Figure 6.Normal probability plots of P 00 and P 11 for July (typical of other months).

Figure 7 .
Figure 7.Comparison of the mean and SD of daily rainfall depths for the five MC models.

Figure 8 .
Figure 8.Comparison of the mean and SD of monthly rainfall depths for the five MC models.
Figure 12 compares the five MC models for the mean and SD of mean length of wet and dry spells at annual resolution.The averages of the absolute values of the Z scores for monthly resolution are shown in Table2.The comparative performances of the five MC models for the mean and SD of mean length of wet spells at monthly (Table2) and annual (Fig.12) resolutions are mostly consistent with their respective performances for mean and SD of number of wet days.All models except HMC preserve the mean and SD of mean length of wet spells, while the HMC tends to overestimate the SD (Fig.12).For the mean and SD of mean length of dry spells, we found that models that can preserve the wet spells distributions also preserve the dry spells distributions and vice versa.As a representative result, the Z scores for the mean and SD of annual mean length of dry spells shown in Fig.12indicate that all models except HMC preserve both mean and SD, while HMC overestimates the SD.Figure12also indi-

Figure 10 .
Figure 10.Comparison of the mean and SD of monthly number of wet days for the five MC models.

Figure 11 .
Figure 11.Comparison of the mean and SD of multiyear number of wet days for the five MC models.

Figure 12 .
Figure 12.Comparison of the mean and SD of annual mean length of wet and dry spells for the five MC models.

Figure 13 .
Figure 13.Comparison of the autocorrelations of monthly number of wet days and monthly rainfall depths for the five MC models for (a) Adelaide and (b) Sydney.The shadings indicate 95 % confidence limits.

Table 1 .
Methodological comparison of the five MC models.

Table 2 .
Average of the absolute values of Z scores (average of the Z scores for all 12 months at daily and monthly resolutions, and average of the Z scores for 1 to 10 years at multiyear resolution) for Sydney (SY) and Adelaide (AD).The averaged Z scores greater than 2 are shown in bold.