Developing Intensity – Duration – Frequency ( IDF ) Curves under Climate Change Uncertainty : The Case of Bangkok , Thailand

The magnitude and frequency of hydrological events are expected to increase in coming years due to climate change in megacities of Asia. Intensity–Duration–Frequency (IDF) curves represent essential means to study effects on the performance of drainage systems. Therefore, the need for updating IDF curves comes from the necessity to gain better understanding of climate change effects. The present paper explores an approach based on spatial downscaling-temporal disaggregation method (DDM) to develop future IDFs using stochastic weather generator, Long Ashton Research Station Weather Generator (LARS-WG) and the rainfall disaggregation tool, Hyetos. The work was carried out for the case of Bangkok, Thailand. The application of LARS-WG to project extreme rainfalls showed promising results and nine global climate models (GCMs) were used to estimate changes in IDF characteristics for future time periods of 2011–2030 and 2046–2065 under climate change scenarios. The IDFs derived from this approach were corrected using higher order equation to mitigate biases. IDFs from all GCMs showed increasing intensities in the future for all return periods. The work presented demonstrates the potential of this approach in projecting future climate scenarios for urban catchment where long term hourly rainfall data are not readily available.


Introduction
Almost on daily basis we can observe that the world is being overwhelmed by extreme hydro-meteorological disasters such as hurricanes (or typhoons), widespread flooding and droughts, all of which are associated with devastating losses and suffering (e.g., [1][2][3][4]).Therefore, our search for optimal configurations of water infrastructure systems represents a great challenge for researchers and practitioners (e.g., [5][6][7][8][9][10]).The uncertainty brought by climate change (or climate extremes) make the above challenge even more significant.Addressing climate change effects in the development of disaster management plans has been on the agenda of many governments and institutions.Hydrodynamic flood modelling is an important aspect in the development of flood mitigation and climate adaptation strategies, which requires high resolution time series data of continuous precipitation or design storms.Hence, a more effective disaster management planning requires both current and projected climate change scenarios for which a range of optimal measures needs to be sought.The use of physical-based models is invaluable for this purpose [11,12].For example, with numerical models, it is possible to explore the generation and propagation of floods and explore any potential impacts that may be incurred due to floods.However, there is a significant body of literature that suggests that a great care and caution should be applied to data collection and processing activities (e.g., rainfall and hydrological data, system data, terrain data, etc.), selection of a model to reflect the problem at hand and the use of geographic information system (GIS) mapping techniques [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30].
Gaining understanding of the potential change of frequency, intensity and volume of extreme rainfall due to climate change has been put forward by many researchers and practitioners (e.g., [31][32][33][34][35]).The need for such understanding comes from the fact that the existing drainage systems are designed to cope with past records of rainfall events and as such they may be insufficient to accommodate future rainfall characteristics.Hence, development of appropriate present and future Intensity-Duration-Frequency (IDF) curves assist in hydrological study for the urban drainage system performance analysis, design and operation.The IDF relationship is a mathematical relationship between rainfall intensity i, duration d, and the return period T [36].The IDF curve is a convenient form of rainfall information and, as such, it has been an important data set for the design of urban drainage systems [37].However, it has been a well-accepted fact that those IDFs developed from past climatic conditions cannot be valid for future climatic conditions unless the IDFs are updated to the future climate trends (see also [38,39]).A study in Mediterranean area, showed that in spite of the important decreases or lesser increment in corresponding total precipitation, extreme rainfall of shorter duration will increase in future [40].Several methods exist for studying changes in IDF under climate change; for examples, equidistance quantile matching method using spatial and temporal downscaling of global climate models (GCMs) and sub-daily observed rainfall in Canada [41]; spatial downscaling and temporal downscaling based on scaling properties of rainfall in urban area of Spain [40]; and spatial-temporal downscaling approach using scaling generalized extreme value (GEV) distribution in Canada [42].Mailhot and Duchesne [38] strongly argue that there is a need to evaluate the climate induced change of pattern, intensity and frequency of extreme rainfall before embarking on any system performance analysis.In this regard, the process of downscaling climate data has been an important task for many researchers and practitioners; and the use of GCMs has been the most appealing method for studying the implications of climate change.However, there are numerous issues associated with the process of downscaling climate data.Despite being important tools to project the expected future scenarios of climatic parameters, GCMs contain biases when compared to observed data due to their parameterization systems and large grid size, which at the regional scale is insignificant, but are significant at basin scale [43].Bias correction methods are applied to correct one or more statistical aspect, which are referred as statistical bias correction, quantile-based bias correction, and histogram equalization methods [43][44][45][46].The present paper provides contribution in this direction and it presents an approach based on spatial downscaling-temporal disaggregation method (DDM), which appears to give promising results in generating rainfall data for a single site station in an urban catchment.

Use of Global Climate Models (GCMs) in Urban Scale Applications
The outputs from GCMs are typically defined at 150-300 km coarse grids, while regional climate models (RCMs) resolutions are about 12-50 km [47].The different downscaling techniques in practice are dynamical downscaling, statistical downscaling, regression based downscaling, weather typing procedure and the stochastic weather generator (see [47][48][49]).The most commonly applied techniques are dynamic and statistical downscaling methods.Long Ashton Research Station Weather Generator (LARS-WG) is a single site numerical model for simulating time series of daily weather.The Weather Generator is a model which, after calibrating site parameters with observed weather data at that site, is capable of simulating synthetic time series of daily weather that are statistically similar to observed weather [50].LARS-WG incorporates projections from 15 GCMs, as used in the Intergovernmental Panel on Climate Change (IPCC) Assessment Report IV, and climate projections are available for the Special Report on Emissions Scenarios (SRES)-SRB1, SRA1B and SRA2 for most of the 15 GCMs [51].The IPCC's latest Assessment report V categorized future projections defined by Representative Concentration Pathways (RCPs) which is based on different 21st century pathways of greenhouse gases (GHG) emissions and atmospheric concentrations, air pollutant emissions and land use.It includes a stringent mitigation scenario (RCP2.6)which keeps global warming likely below 2 • C above pre-industrial temperatures, two intermediate scenarios (RCP4.5 and RCP6.0) and one scenario with very high GHG emissions (RCP8.5)[52].This study is based on analysis and modelling scenarios based on earlier IPCC Assessment report IV-Special Report on Emission Scenarios (SRES).The Assessment report IV categorized the scenarios of future climate based on carbon emission levels and climate policy.The future scenarios can be summarized into four storylines: A1, A2, B1 and B2, where the A2 family is the most extreme scenario owing to its higher amount of carbon emission [53].In terms of overall forcing, RCP8.5 is broadly comparable to the SRES A2/A1FI scenario, RCP6.0 to B2 and RCP4.5 to B1 [52].
Although it is common practice to downscale the results from GCMs or RCMs to urban scales, and despite considerable improvements in computational power in recent years, climate models still remain relatively coarse in space and temporal resolution, and are unable to resolve significant features at the finer scales of urban drainage systems (see [43,54]).In most places, there are numerous rain gauges that have been in operation for several decades, but most of them measure daily rainfall while recordings at hourly or sub-hourly data are not so common.Any hydrological study requires rainfall data at different time scales, depending on the purpose.The studies of water supply and irrigation require data on a monthly time scale, while studies on rainwater harvesting and river flow estimation require data on a daily time scale.However, for urban drainage and pluvial flood studies, it is necessary to have the rainfall data on sub-daily time scales as hydrological processes in urban drainage systems and small catchments occur at moderately smaller spatial and temporal scales [55].Hydrological study in a small urban catchment may require rainfall data at spatial and temporal resolution of 1 km 2 and 5 min respectively.One of the reasons for the lack of availability of high temporal scale rainfall data is the cost associated with collection of such data.Due to this explicit requirement of small scale and short duration rainfall, studies dealing with climate change impacts on urban drainage are particularly challenging [56].The disaggregation of available coarse rainfall data to a small duration can be one way to address this limitation.Koutsoyiannis [57] developed a method to utilize the available daily rainfall data and disaggregate it to hourly scale.Furthermore, numerous disaggregation models have been developed to disaggregate monthly data to daily scale for a single site stations [57].
Stochastic disaggregation models for rainfall disaggregation can be divided into two categories: first is the point process theory based on Bartlett-Lewis and Neyman-Scott rectangular pulses; the second is the scale invariance theory of cascade, fractal and multifractal approaches [58].Most of the disaggregation techniques developed before 2001 are essentially single site techniques.However, while multiple site rainfall disaggregation (as a means of simultaneous spatial and temporal disaggregation) is of a significant interest to many researchers, this technique comes with issues such as increased mathematical complexity when compared to a single site disaggregation [57].Koutsoyiannis and Onof [59] addressed the problem of single site rainfall disaggregation to a fine scale (of up to 1 h) by developing a stochastic model based on the Bartlett-Lewis rainfall model referred to as Hyetos.Kossieris et al. [60] developed the improved version of Hyetos (Version 2.0) for its application to sub-hourly time scales.In the present paper, we discuss the results obtained from the Hyetos application (Version 0.3, developed by Koutsoyiannis and Onof, in Imperial College London) for continuous series hourly rainfall generation from daily observed and downscaled future data.
An example of applying downscaled data for IDF generation is given in the study by Liew et al. [61] who developed present and future IDF curves for the cities of Singapore, Kuala Lumpur and Jakarta.Since the data required were scarce, dynamic downscaling using RCM, the Weather Research and Forecasting (WRF) model and 1957 to 2002 dataset from ERA-40 (European Centre for Medium range Weather Forecasting) were used.The results from their study showed significant underestimation of rainfall intensities in WRF-ERA-40 generated data.Bias correction was applied for downscaled future data from 2071 to 2100.The study concluded that there is a general increment in rainfall intensities throughout the future period.The study applied linear adjustment of the IDF curve to correct biases from the model's results.The application of linear adjustments is relevant when the duration considered is at coarse scales: for instance, from 6 h to 24 h, the curve will be of a linear nature.However, the mathematical representation of IDFs for shorter duration rainfall requires higher order equations.One of the challenges in urban drainage modelling is that IDFs with very short duration are required owing to the corresponding time of concentration within the catchment.As collecting continuous short duration climate variables requires significant resources, projections of future climate variables using climate downscaling techniques is insufficient as output is generally on daily scale.To generate future IDFs for urban hydrological studies, such coarse rainfall data requires disaggregation to a finer scale.
An approach developed based on spatial downscaling-temporal disaggregation (DDM) is presented in the following sections of this paper and its usefulness is demonstrated on the case study of Bangkok, Thailand.

Case Study Area
Bangkok, the capital and commercial hub of Thailand, is one of the highly urbanized cities in Southeast Asia.The selected Bangkok metropolis rainfall station is located in Sukhumvit area, which is in the central part of Bangkok metropolis (Figure 1).Sukhumvit area is surrounded by major canals on the eastern and northern sides which connects to Chao Phraya River on the southern side.The geographical location of the study catchment is 13 • 44 18.01" N latitudes and 100 • 33 41.31" E longitude.The area is flat, with an average altitude of 0.4 m to 1 m above mean sea level.The climate is tropical wet and dry with long hours of sunshine, fairly high temperatures year-round and high humidity.The rainy season is observed in mid-May which continues until October.The average minimum temperature is approximately 20 • C which occurs during November-February and maximum temperature could reach up to 35 • C during March-June.Bangkok receives an average annual rainfall of 1500 mm and is influenced by seasonal monsoons.The drainage outlets from the Sukhumvit catchment area are connected to canals which discharge further into the Chao Phraya River.The time of concentration in the urban catchment of Sukhumvit area was observed to be 2.5 h (see [62]).Bangkok suffered from severe floods in 1942, 1983, 1995 and 2011.The flood event in 2011 was the worst flood event in Thailand over the past 50 years.Climate risks from flooding and potential sea level rise are causing Bangkok to be one of the most vulnerable coastal cities in the world [63].The temperature and precipitation in dense urban parts of Bangkok are expected to increase in the future [64], with increase in the extent of flood prone areas and flood volume [65].
Water 2017, 9, 145 4 of 21 from the model's results.The application of linear adjustments is relevant when the duration considered is at coarse scales: for instance, from 6 h to 24 h, the curve will be of a linear nature.However, the mathematical representation of IDFs for shorter duration rainfall requires higher order equations.One of the challenges in urban drainage modelling is that IDFs with very short duration are required owing to the corresponding time of concentration within the catchment.As collecting continuous short duration climate variables requires significant resources, projections of future climate variables using climate downscaling techniques is insufficient as output is generally on daily scale.To generate future IDFs for urban hydrological studies, such coarse rainfall data requires disaggregation to a finer scale.
An approach developed based on spatial downscaling-temporal disaggregation (DDM) is presented in the following sections of this paper and its usefulness is demonstrated on the case study of Bangkok, Thailand.

Description of the Case Study Area
Bangkok, the capital and commercial hub of Thailand, is one of the highly urbanized cities in Southeast Asia.The selected Bangkok metropolis rainfall station is located in Sukhumvit area, which is in the central part of Bangkok metropolis (Figure 1).Sukhumvit area is surrounded by major canals on the eastern and northern sides which connects to Chao Phraya River on the southern side.The geographical location of the study catchment is 13°44′18.01″N latitudes and 100°33′41.31″E longitude.The area is flat, with an average altitude of 0.4 m to 1 m above mean sea level.The climate is tropical wet and dry with long hours of sunshine, fairly high temperatures year-round and high humidity.The rainy season is observed in mid-May which continues until October.The average minimum temperature is approximately 20 °C which occurs during November-February and maximum temperature could reach up to 35 °C during March-June.Bangkok receives an average annual rainfall of 1500 mm and is influenced by seasonal monsoons.The drainage outlets from the Sukhumvit catchment area are connected to canals which discharge further into the Chao Phraya River.The time of concentration in the urban catchment of Sukhumvit area was observed to be 2.5 h (see [62]).Bangkok suffered from severe floods in 1942, 1983, 1995 and 2011.The flood event in 2011 was the worst flood event in Thailand over the past 50 years.Climate risks from flooding and potential sea level rise are causing Bangkok to be one of the most vulnerable coastal cities in the world [63].The temperature and precipitation in dense urban parts of Bangkok are expected to increase in the future [64], with increase in the extent of flood prone areas and flood volume [65].

Methodological Framework for the Spatial Downscaling-Temporal Disaggregation Method (DDM) Approach
To produce future IDF curves for Bangkok a methodological framework was developed (Figure 2).This framework consists of three parts: spatial downscaling of rainfall obtained from GCMs, temporal disaggregation of daily rainfall series into hourly series and the application of generalized extreme value (GEV) distribution with annual maximum series to generate IDF curves for the present and the future.

Methodological Framework for the Spatial Downscaling-Temporal Disaggregation Method (DDM) Approach
To produce future IDF curves for Bangkok a methodological framework was developed (Figure 2).This framework consists of three parts: spatial downscaling of rainfall obtained from GCMs, temporal disaggregation of daily rainfall series into hourly series and the application of generalized extreme value (GEV) distribution with annual maximum series to generate IDF curves for the present and the future.The rainfall data of 3 h and 24 h, from 1981 to 2010 were acquired from the Bangkok metropolis rainfall station of Thai Meteorological Department.These data were utilized for the present IDFs generation, estimation of Bartlett-Lewis Rectangular Pulse (BLRP) parameters, and inputs for LARS-WG and Hyetos.Each of the methodological details is described in the following Sections 4.2 to 4.4.

Spatial Downscaling Using Long Ashton Research Station Weather Generator (LARS-WG)
LARS-WG uses a semi-empirical distribution (SED), the cumulative probability distribution function (PDF) with 23 intervals, to approximate probability distributions of dry and wet series of daily precipitation in each month [51].The semi-empirical distribution Emp = , ; ℎ , = 1, … , 10 is a histogram with 10 intervals, , where < , and ℎ denotes the number of events from the observed data in the i-th interval.Each interval contains a number of events based on the observed data.For each climatic variable v, a value of a climatic variable corresponding to the probability is calculated as [51]: where P() denotes probability based on observed data { }, = 0 and = 1 and these are fixed parameters with corresponding values of = min{ } and = max{ }.The simulation of precipitation occurrence is modeled as alternate wet (precipitation > 0.0 mm) and dry series.The length of each series is chosen randomly from the wet or dry semi-empirical distribution for the month in which the series start.For a wet day, the precipitation value is generated from the semiempirical precipitation distribution for the particular month, independent of the length of the wet series or the amount of precipitation on previous days.The methodology for spatial downscaling using LARS-WG to generate future climate variables includes three major steps: site analysis, calibration i.e. comparison of the statistical properties of the synthetic data with observed data; and generation of climate data.First, an analysis of the observed weather parameters such as temperature, rainfall and solar radiation was done and second, the synthetic daily weather data are generated utilizing these weather parameters.Site analysis involved the LARS-WG model calibration process, which is done to determine the statistical characteristics of the observed data.The synthetic weather of similar time period as observed data is generated using the parameter file derived during the model calibration step, to validate the model.The statistical characteristics determined from the observed weather data are used to generate synthetic weather data, which have the same statistical The rainfall data of 3 h and 24 h, from 1981 to 2010 were acquired from the Bangkok metropolis rainfall station of Thai Meteorological Department.These data were utilized for the present IDFs generation, estimation of Bartlett-Lewis Rectangular Pulse (BLRP) parameters, and inputs for LARS-WG and Hyetos.Each of the methodological details is described in the following Sections 4.2-4.4.

Spatial Downscaling Using Long Ashton Research Station Weather Generator (LARS-WG)
LARS-WG uses a semi-empirical distribution (SED), the cumulative probability distribution function (PDF) with 23 intervals, to approximate probability distributions of dry and wet series of daily precipitation in each month [51].The semi-empirical distribution Emp = {a 0 , a i ; h i , i = 1, . . ., 10} is a histogram with 10 intervals, [a i−1 , a i ] where a i−1 < a i , and h i denotes the number of events from the observed data in the i-th interval.Each interval contains a number of events based on the observed data.For each climatic variable v, a value of a climatic variable v i corresponding to the probability p i is calculated as [51]: where P() denotes probability based on observed data {v obs }, p 0 = 0 and p n = 1 and these are fixed parameters with corresponding values of v 0 = min{v obs } and v n = max{v obs }.The simulation of precipitation occurrence is modeled as alternate wet (precipitation > 0.0 mm) and dry series.The length of each series is chosen randomly from the wet or dry semi-empirical distribution for the month in which the series start.For a wet day, the precipitation value is generated from the semi-empirical precipitation distribution for the particular month, independent of the length of the wet series or the amount of precipitation on previous days.The methodology for spatial downscaling using LARS-WG to generate future climate variables includes three major steps: site analysis, calibration i.e. comparison of the statistical properties of the synthetic data with observed data; and generation of climate data.First, an analysis of the observed weather parameters such as temperature, rainfall and solar radiation was done and second, the synthetic daily weather data are generated utilizing these weather parameters.Site analysis involved the LARS-WG model calibration process, which is done to determine the statistical characteristics of the observed data.The synthetic weather of similar time period as observed data is generated using the parameter file derived during the model calibration step, to validate the model.The statistical characteristics determined from the observed weather data are used to generate synthetic weather data, which have the same statistical properties as observed weather.Climate data are generated based on site parameters and climate scenarios from GCMs of interest.A scenario file is required for this process which can be manually edited and prepared or can be derived automatically from GCMs' data in grid.The latest version of the stochastic weather generator was developed as LARS-WG 5.0, which is able to generate high resolution climate change scenarios over a region using direct outputs from GCMs.LARS-WG incorporates predictions from 15 GCMs, namely BCM2, CGMR, CNCM3, CSMK3, FGOALS, GFCM21, GIAOM, HADCM3, HADGEM, INCM3, IPCM4, MIHR, MPEH5, NCCCSM and NCPCM; and climate predictions are available for IPCC Assessment Report IV's SRES emissions scenarios SRB1, SRA1B and SRA2 for most of the GCMs (see [51]).In this study, two scenarios, SRA1B and SRA2, were considered and future projections of these scenarios are available from 9 GCMs, namely, CNCM3, GFCM21, HADCM3, HADGEM, INCM3, IPCM4, MPEH5, NCCCSM and NCPCM (Table 1).

Disaggregation to Hourly Time Scale Using Hyetos
Hyetos is based on the Bartlett-Lewis Rectangular Pulse (BLRP) process theory which has been developed into a computer program and it is a proven disaggregation technique which adjusts finer scale data in order to obtain the coarser scale data without affecting stochastic structure by the model (see [57]).The general assumptions of the model are: storm origins occur following a Poisson process (rate lambda-λ), origins of cells of each storm arrive following a Poisson process (rate Beta-β), arrivals of each storm terminate after a time exponentially distributed (parameter gamma-γ), each cell has a duration exponentially distributed (parameter eta-η), and each cell has a uniform intensity with a specified (exponential or gamma) distribution.The parameter η is randomly varied from storm to storm with a gamma distribution with a shape parameter alpha-α and a scale parameter v. Subsequently, parameters β and γ also vary such that the ratios k = β/η and φ = γ/η stay constant.The distribution of uniform intensity X ij is typically assumed to be exponential with the parameter 1/µX.Alternatively, two-parameter gamma can be assumed with mean µX and standard deviation σX.The relationship for mean, variance, lag 1 autocovariance and proportion of dry periods, depending upon these variables of the Bartlett-Lewis parameters at the particular time had been derived into equations as presented in [57].Hyetos adopts the Bartlett-Lewis model for the generation of synthetic rainfall along with an adjusting procedure to preserve consistency with observed daily data.For a sequence of wet days, model runs several times and generates a sequence that best matches with the observed daily data.The sequence of synthetic hourly rainfall is modified according to the proportional adjusting procedure, which aggregates to daily data, such that certain statistics are preserved among daily and hourly data [60].In addition to adjusting procedure, the model uses repetition to ensure similarity in observed and modelled data.For the longer series, the model divides the wet series into independent cluster separated by at least one day, which reduces the computational time.The model runs separately for each wet cluster series, until the departures of the sequence of daily sums from the given sequence of daily rainfall becomes lower than an acceptable limit.For the sequence of L wet days generated, the daily rainfall depths are calculated and compared with the observed ones by logarithmic distance, which is mathematically given as [59]: where Z i and Z i are observed and generated daily rainfall depth of day i of wet day sequence, respectively, and c is constant (=0.1 mm).The model verifies if departure d is smaller than acceptable limit d a , which continues for the allowed number of repetitions and if the condition is not satisfied, the cell is discarded and new one is generated.

Intensity-Duration-Frequency (IDF) Generation and Correction
There are different distribution functions for IDF analysis: Extreme Value Type I, i.e., Gumbel (EVI) distribution, Generalized extreme value (GEV) distribution, Gamma distribution, Log Pearson III distribution, Lognormal distribution, Exponential distribution, and Pareto distribution.Gumbel distribution function (see for example, [36]) and GEV distribution (see [42]) are the most commonly used function for IDF analysis.In this present work, the annual maximum series from observed rainfall data (1981-2010) were fitted into GEV, Log Pearson III, Lognormal, Gamma and Gumbel distributions using method of moments and L-moments.The mathematical expression for GEV probability distribution is given by, where k (shape parameter), σ (scale parameter, σ > 0) and µ (location parameter) are to be determined.For k = 0, GEV reduces to Gumbel (EVI) distribution and for k < 0, it reduces to Extreme Value Type II distribution, which implies an upper bound of the variable, which is not the case in maximum rainfall intensity [36].
The method of moments, first developed by Karl Pearson in 1902, considers that the good estimates of the parameters of a probability distribution are those for which moments of the probability density function about the origin are equal to the corresponding moments of the sample data [66].L-moments are the improvement over the conventional moments, which has major advantage, being linear functions of the data, suffer less from the effects of sampling variability and are more robust than conventional moments to outliers in the data.It is also less subjected to bias in estimation and enable more secure inferences to be made from small samples about an underlying probability distribution [67].Let X 1:n ≤ X 2:n ≤ . . .≤ X n:n be an order statistics of the random sample of size n from distribution X, and X is the real valued random variable with cumulative distribution function F(x), the L-moments can be estimated by, r = 1, 2, . . ., where the expectation of an order statistic can be expressed as, For the details of L-moments and parameters for different probability distributions, readers are referred to Hosking and Wallis [68].
For the selection of distribution functions, the goodness of fit of the probability distributions are carried out by comparative test of theoretical and sample values of the relative frequency or the cumulative frequency function, using Kolmogorov-Smirnov (K-S) test and Chi-Squared test to determine if the annual maximum rainfall data follow the specified distribution.The K-S test is based on the largest vertical difference between the theoretical and the empirical cumulative distribution function.For the given sample size of x 1 , ... , x n , of the distribution with cumulative distribution function (CDF) F(x), the K-S test statistic is expressed as, The Chi-squared test statistics is defined as, where Observed i is observed number of occurrence in interval i, is the corresponding expected number of occurrences in interval i, F is the cumulative distribution function being tested and k = 1 + log 2 N(sample size).The p-value from the goodness of fit tests are calculated based on test statistics of each fitted distribution, which compares the suitability of each distribution.Finally, the graphical comparison of IDFs from observed and disaggregated data and the corrections carried out using suitable equation in this study are discussed in Section 5.

Spatial Downscaling
In this study, the observed precipitation data were analysed to determine their statistical characteristics and to compute site parameters.Weather generation was done for 1981-2010.The total monthly mean rainfall and maximum daily rainfall for each month were selected as statistics for comparing observed and generated data (Figure 3).
The Chi-squared test statistics is defined as, where Observed is observed number of occurrence in interval i, Expected = F( ) − F( ) is the corresponding expected number of occurrences in interval i, F is the cumulative distribution function being tested and = 1 + log ( sample size) .The p-value from the goodness of fit tests are calculated based on test statistics of each fitted distribution, which compares the suitability of each distribution.Finally, the graphical comparison of IDFs from observed and disaggregated data and the corrections carried out using suitable equation in this study are discussed in Section 5.

Spatial Downscaling
In this study, the observed precipitation data were analysed to determine their statistical characteristics and to compute site parameters.Weather generation was done for 1981-2010.The total monthly mean rainfall and maximum daily rainfall for each month were selected as statistics for comparing observed and generated data (Figure 3).The Root Mean Square Errors (RMSE) and Efficiency Index (EI) have been used as standard metrics to show that the model is capable of reproducing past rainfall.RMSE and EI are mathematically expressed as: = 1 (10) where is observed data at time i, is mean of observed data, is synthetic data at time i and n is the total number of data points.The RMSE for both mean and maxima are in acceptable range and The Root Mean Square Errors (RMSE) and Efficiency Index (EI) have been used as standard metrics to show that the model is capable of reproducing past rainfall.RMSE and EI are mathematically expressed as: where X i is observed data at time i, X is mean of observed data, Y i is synthetic data at time i and n is the total number of data points.The RMSE for both mean and maxima are in acceptable range and further EI are over 90%.Developing a set of IDFs generally depends on extreme rainfall and the LARS-WG results showed its ability to project extreme rainfall as close to the observed rainfall maxima, hence, LARS-WG was used for future downscaling of climate variable (i.e., rainfall series).
The above-mentioned nine GCMs (in Section 4.4) (for SRA1B and SRA2) were used in the analysis to estimate range of likely changes in future IDFs.

Temporal Disaggregation
The method follows the assumption of the BLRP process that can be derived into expressions defined by parameters-λ, k, φ, α, v, µX and σX-which mathematically elucidate the rainfall event (see also, [57,59]).The method aims to preserve statistical characteristics such as mean, variance, lag 1 autocovariance and proportion of dry periods of observed time series with the synthetically generated disaggregated time series.The estimation of the BLRP parameters was done by calculating these statistical characteristics from the observed 3 h, 12 h and 24 h data from 1986 to 2000 for each month, separately (Figure 4).A period between 1986 and 2000 was selected for the analysis since maximum daily rainfall occurred in the same period.The parameters λ, k, φ, α, v, µX and σX were optimized to minimize relative error between synthetic and observed values.The synthetic rainfall generated preserves the statistical properties of observed rainfall series.
Water 2017, 9, 145 9 of 21 further EI are over 90%.Developing a set of IDFs generally depends on extreme rainfall and the LARS-WG results showed its ability to project extreme rainfall as close to the observed rainfall maxima, hence, LARS-WG was used for future downscaling of climate variable (i.e., rainfall series).The above-mentioned nine GCMs (in Section 4.4) (for SRA1B and SRA2) were used in the analysis to estimate range of likely changes in future IDFs.

Temporal Disaggregation
The method follows the assumption of the BLRP process that can be derived into expressions defined by parameters-λ, k, φ, α, v, μX and σX-which mathematically elucidate the rainfall event (see also, [57,59]).The method aims to preserve statistical characteristics such as mean, variance, lag 1 autocovariance and proportion of dry periods of observed time series with the synthetically generated disaggregated time series.The estimation of the BLRP parameters was done by calculating these statistical characteristics from the observed 3 h, 12 h and 24 h data from 1986 to 2000 for each month, separately (Figure 4).A period between 1986 and 2000 was selected for the analysis since maximum daily rainfall occurred in the same period.The parameters λ, k, φ, α, v, μX and σX were optimized to minimize relative error between synthetic and observed values.The synthetic rainfall generated preserves the statistical properties of observed rainfall series.Optimization of the variables to minimize weighted error between synthetic and historical observed data is shown in Figure 4.The overall error in simulating disaggregated data is below 2 mm for 3 h, 12 h and 24 h, with satisfactory optimization of mean, variance, lag 1 autocovariance and proportion of dry periods.The optimized parameters for each month (Table 2) were used to disaggregate present observed data and future downscaled daily data.The same parameters have been used for past and future disaggregations.The maximum 3 h precipitation for each month is compared between disaggregated and observed data for the period 1981-2010 (Figure 5).The disaggregated data indicated under prediction, especially in the months of April, May, September and October.Optimization of the variables to minimize weighted error between synthetic and historical observed data is shown in Figure 4.The overall error in simulating disaggregated data is below 2 mm for 3 h, 12 h and 24 h, with satisfactory optimization of mean, variance, lag 1 autocovariance and proportion of dry periods.The optimized parameters for each month (Table 2) were used to disaggregate present observed data and future downscaled daily data.The same parameters have been used for past and future disaggregations.The maximum 3 h precipitation for each month is compared between disaggregated and observed data for the period 1981-2010 (Figure 5).The disaggregated data indicated under prediction, especially in the months of April, May, September and October.Optimization of the variables to minimize weighted error between synthetic and historical observed data is shown in Figure 4.The overall error in simulating disaggregated data is below 2 mm for 3 h, 12 h and 24 h, with satisfactory optimization of mean, variance, lag 1 autocovariance and proportion of dry periods.The optimized parameters for each month (Table 2) were used to disaggregate present observed data and future downscaled daily data.The same parameters have been used for past and future disaggregations.The maximum 3 h precipitation for each month is compared between disaggregated and observed data for the period 1981-2010 (Figure 5).The disaggregated data indicated under prediction, especially in the months of April, May, September and October.

Generation of IDF Curves
The annual maximum rainfalls at different durations from 1981 to 2010 were fitted into Generalized extreme value (GEV) distribution, Gamma distribution, Gumbel distribution, Log Pearson III distribution and Lognormal distribution.The goodness of fit test (Table 3) using Kolgomorov-Smirnov (K-S) test showed that GEV distribution fits better for annual maximum rainfall at 3 h, 6 h and 24 h durations, while for 12 h duration Log normal distribution fits better, although GEV also fits as closely to Log normal distribution.The Chi-squared (χ2) test results showed GEV fits better for annual maximum rainfall at 12 h and 24 h durations, while GEV and Gamma distribution fits better at 3 h duration and Gumbel distribution fits better at 6 h duration.In overall comparison, GEV seemed to fit better than other distribution so we have considered GEV distribution to generate present and future IDFs.The future annual maximum rainfall for 2011-2030 and 2046-2065 at different rainfall durations of 1 h, 3 h, 6 h, 12 h and 24 h also fit well with GEV distribution, the distribution analysis for future are not shown here due to brevity of the paper.Table 3 also shows the comparison of rainfall intensities at 2-, 5-and 20-years return periods when annual maximum rainfalls are fitted into different distribution.The result showed minimal differences in the rainfall intensities.The standard deviation of rainfall intensities however increases accordingly to higher return periods.
For the comparison of disaggregated rainfall (daily to hourly) with the observed 3 h rainfall, the disaggregated 1 h data were aggregated to 3 h data, from 1981 to 2010, and IDFs were generated.The disaggregated data (i.e., modelled data) showed significant underestimation as compared to the observed data for the duration of 3 h; however, the underestimation was observed to be relatively lower for longer durations (Figure 6).Thus, the correction factors were used to correct the bias for modelled rainfall intensities (Table 4).The differences and ratios between observed and simulated intensities were found to be larger for shorter durations and smaller for durations of longer than 3 h.For all durations, the differences increased with return periods.Nevertheless, the ratio remained almost equal for different return periods for a specific duration of rainfall.

Generation of IDF Curves
The annual maximum rainfalls at different durations from 1981 to 2010 were fitted into Generalized extreme value (GEV) distribution, Gamma distribution, Gumbel distribution, Log Pearson III distribution and Lognormal distribution.The goodness of fit test (Table 3) using Kolgomorov-Smirnov (K-S) test showed that GEV distribution fits better for annual maximum rainfall at 3 h, 6 h and 24 h durations, while for 12 h duration Log normal distribution fits better, although GEV also fits as closely to Log normal distribution.The Chi-squared (χ2) test results showed GEV fits better for annual maximum rainfall at 12 h and 24 h durations, while GEV and Gamma distribution fits better at 3 h duration and Gumbel distribution fits better at 6 h duration.In overall comparison, GEV seemed to fit better than other distribution so we have considered GEV distribution to generate present and future IDFs.The future annual maximum rainfall for 2011-2030 and 2046-2065 at different rainfall durations of 1 h, 3 h, 6 h, 12 h and 24 h also fit well with GEV distribution, the distribution analysis for future are not shown here due to brevity of the paper.Table 3 also shows the comparison of rainfall intensities at 2-, 5-and 20-years return periods when annual maximum rainfalls are fitted into different distribution.The result showed minimal differences in the rainfall intensities.The standard deviation of rainfall intensities however increases accordingly to higher return periods.
For the comparison of disaggregated rainfall (daily to hourly) with the observed 3 h rainfall, the disaggregated 1 h data were aggregated to 3 h data, from 1981 to 2010, and IDFs were generated.The disaggregated data (i.e., modelled data) showed significant underestimation as compared to the observed data for the duration of 3 h; however, the underestimation was observed to be relatively lower for longer durations (Figure 6).Thus, the correction factors were used to correct the bias for modelled rainfall intensities (Table 4).The differences and ratios between observed and simulated intensities were found to be larger for shorter durations and smaller for durations of longer than 3 h.For all durations, the differences increased with return periods.Nevertheless, the ratio remained almost equal for different return periods for a specific duration of rainfall.The average ratios of observed to modelled data at different rainfall duration were fitted into a power and polynomial equations and the resulting equations are the function of rainfall duration, which is independent from the return period (Figure 7).Coefficient of Determination (R 2 ) for the fitted polynomial equation showed better accuracy than other equation.Using the equation, where y is the correction factor for a specific rainfall duration of x in hours, the downscaled-disaggregated IDFs were graphically corrected.IDFs generated from nine GCMs were corrected using the respective correction factors.The corrected IDFs for 1981-2010 are shown in Figure 6.The average ratios of observed to modelled data at different rainfall duration were fitted into a power and polynomial equations and the resulting equations are the function of rainfall duration, which is independent from the return period (Figure 7).Coefficient of Determination (R 2 ) for the fitted polynomial equation showed better accuracy than other equation.Using the equation, = 0.0017 − 0.0606 + 1.5417 (11) where y is the correction factor for a specific rainfall duration of x in hours, the downscaleddisaggregated IDFs were graphically corrected.IDFs generated from nine GCMs were corrected using the respective correction factors.The corrected IDFs for 1981-2010 are shown in Figure 6.(i) ( Figure 8 presents the future IDFs projections shown by box and whiskers plot for three quartile regions, which is compared with baseline (1981-2010) and Table 5 shows the detail changes in future intensities compared to baseline period.Comparison of present and future IDFs for return periods (2-, 5-and 20-years), showed that intensities increases under climate change at most of the rainfall durations.In comparison to the rainfall intensities the uncertainties among GCMs are found to be higher for smaller duration, which then increases with higher return periods.Even the median of the projected nine GCMs data showed the increases in future rainfall intensities when compared to the present IDFs, in all durations except for 1 h duration in 20 years return period.The maximum projected by some GCMs shows future rainfall intensity increments also in 1 h, while even minimum projected by some GCMs is higher than the present scenario, except at smaller durations in 20 years return period.There is no significant difference among the results from SRA1B and SRA2 scenarios in 2-and 5-years return periods, while SRA2 showed maximum increments in 20 years return period.It can be observed that few of the GCMs in 20 years return period at 1 h duration projected decrease in rainfall intensities.Compared to 1981-2010, during 2011-2030 in 20 years return period, the maximum decrease in 1 h rainfall intensities is projected by NCPCM in SRA1B (−7.8%) and GFCM21 in SRA2 (−10.1%).Similarly, for 2046-2065 in 20 years return period, the maximum decrease in 1 h rainfall intensities are projected by IMCM3 and HADGM3 in SRA1B (−10.2%) and HADGEM in SRA2 (−18.3%).The maximum increase in 1 h rainfall intensities during 2011-2030 are projected by NCCCSM in SRA1B (20.4%) and IMCM3 in SRA2 (20.4%) in 2 years return period, NCCCSM in SRA1B (13.9%) and INCM3 in SRA2 (18.7%) in 5 years return period, GFCM21 in SRA1B (3.7%) and NCCCSM in SRA2 (25.9%) in 20 years return period.The maximum increase in 1h rainfall intensities during 2046-2065 are projected by NCCCSM in SRA1B (25.6%) and IMCM3 in SRA2 (15.7%) in 2 years return period, NCCCSM in SRA1B (18.7%) and HADCM3 in SRA2 (14.3%) in 5 years return period, IPCM4 in SRA1B (10.2%) and HADCM3 in SRA2 (14.5%) in 20 years return period.Besides, it is noticeable that all GCMs showed only positive increments in rainfall intensities for higher   5 shows the detail changes in future intensities compared to baseline period.Comparison of present and future IDFs for return periods (2-, 5-and 20-years), showed that intensities increases under climate change at most of the rainfall durations.In comparison to the rainfall intensities the uncertainties among GCMs are found to be higher for smaller duration, which then increases with higher return periods.Even the median of the projected nine GCMs data showed the increases in future rainfall intensities when compared to the present IDFs, in all durations except for 1 h duration in 20 years return period.The maximum projected by some GCMs shows future rainfall intensity increments also in 1 h, while even minimum projected by some GCMs is higher than the present scenario, except at smaller durations in 20 years return period.There is no significant difference among the results from SRA1B and SRA2 scenarios in 2-and 5-years return periods, while SRA2 showed maximum increments in 20 years return period.It can be observed that few of the GCMs in 20 years return period at 1 h duration projected decrease in rainfall intensities.Compared to 1981-2010, during 2011-2030 in 20 years return period, the maximum decrease in 1 h rainfall intensities is projected by NCPCM in SRA1B (−7.8%) and GFCM21 in SRA2 (−10.1%).Similarly, for 2046-2065 in 20 years return period, the maximum decrease in 1 h rainfall intensities are projected by IMCM3 and HADGM3 in SRA1B (−10.2%) and HADGEM in SRA2 (−18.3%).The maximum increase in 1 h rainfall intensities during 2011-2030 are projected by NCCCSM in SRA1B (20.4%) and IMCM3 in SRA2 (20.4%) in 2 years return period, NCCCSM in SRA1B (13.9%) and INCM3 in SRA2 (18.7%) in 5 years return period, GFCM21 in SRA1B (3.7%) and NCCCSM in SRA2 (25.9%) in 20 years return period.The maximum increase in 1 h rainfall intensities during 2046-2065 are projected by NCCCSM in SRA1B (25.6%) and IMCM3 in SRA2 (15.7%) in 2 years return period, NCCCSM in SRA1B (18.7%) and HADCM3 in SRA2 (14.3%) in 5 years return period, IPCM4 in SRA1B (10.2%) and HADCM3 in SRA2 (14.5%) in 20 years return period.Besides, it is noticeable that all GCMs showed only positive increments in rainfall intensities for higher duration in all return periods and scenarios (Figure 8).This is possibly due to the limitation and potential biases in the methods in preserving extreme values of disaggregated data at finer scales and limitation of extrapolating correction factor to 1 h using fitted correction equation.
Figure 8 showed the range of uncertainty, through box and whisker plots, in agreement among GCMs outputs.By considering multiple GCMs and two future scenarios we quantified the range of future uncertainty.Again considering 3 h rainfall, the future rainfall intensities will increase compared to present condition by 2.4% to 27.8% in 2 years return period and 1.3% to 30.4% in 20 years return period during 2011-2030; and 5.4% to 27.0% in 2 years return period and 8.5% to 41.2% in 20 years return period during 2046-2065.These increments in future rainfall intensities will affect performance of urban drainage systems which increase further challenges in managing urban storm water in Bangkok.The GEV distribution fitted better for annual maximum series in all duration for present as well as future.A study using Gumbel distribution for same station [62], the results of the study which is not presented here showed similar behaviour from the selected GCMs, however with some changes in projected rainfall intensities.The selected future scenarios SRA1B and SRA2 represents medium to high greenhouse gases emission scenarios, and RCP8.5 is broadly comparable to SRA2.However, the individual GCMs outputs, in two scenarios, in term of fitted rainfall intensities showed no precise differences.The highest future increments in rainfall intensities were projected under SRA2 scenario in 2011-2030 and under SRA1B scenario in 2046-2065.The temporal disaggregation of daily to hourly rainfall series is important part of the DDM approach.One of the major advantage of using temporal disaggregation is its ability to support urban hydrological studies and its practical applications, using only short term higher resolution data together with fairly coarser resolution long term data that are available in most cities.The application of robust disaggregation techniques and analysis under newer future climate scenario will further improve this approach.The DDM approach is applicable and relevant to other urban catchments in different climatic conditions where availability of higher resolution data is constraint.

Conclusions
The present paper presents an approach based on spatial downscaling-temporal disaggregation method (DDM)-to develop future IDFs using stochastic weather generator, Long Ashton Research Station Weather Generator (LARS-WG) and the rainfall disaggregation tool, Hyetos.The work was carried out for the case of Bangkok, Thailand.The results obtained show that the DDM approach is capable of producing promising results that can be applied to other catchments.Downscaling with LARS-WG gave reasonable results in estimating extreme climate variable (i.e., rainfall series).The downscaled-disaggregated IDFs showed underestimation in intensities, especially for short durations; hence, it was found necessary to correct the bias in IDFs generation.
The results obtained show that IDFs for shorter durations, which exhibit non-linear properties compared to general IDF in higher durations (e.g., 12 h to 24 h) which exhibit linear properties, are corrected using higher order equation.With the availability of only 3 h interval observed rainfall data for the study, the correction factor for 1 h is extrapolated with the same equation.Nine GCMs were analysed which incorporate future projections in SRA1B and SRA2, and all of them projected higher intensities for durations from 3 h to 24 h.The few discrepancies at 1 h rainfall are due to limitations of spatial downscaling and temporal disaggregation method in preserving maxima.The overall findings indicate an increase in future rainfall for Bangkok region.Most of the rainfall stations globally lack long-term higher temporal resolution rainfall data for long periods while most GCMs outputs are in daily time scale.In places where such data availability is an issue, the DDM approach can provide quantified evidence of climate change implication on rainfall intensities.The short term sub-daily scale rainfall at the particular rainfall station or nearby stations is needed for generation of BLRP parameters using Hyetos.However, BLRP parameters will vary in relation to the location and nature of rainfall data.Our future research in this direction will focus on improving temporal disaggregation methods and achieving more robust outputs of future climate downscaling models in higher resolutions.We also intend to incorporate the subsequent IPCC emissions scenarios in our future research work.

Figure 1 .
Figure 1.Location of the case study area: Bangkok, Thailand.

Figure 1 .
Figure 1.Location of the case study area: Bangkok, Thailand.

Figure 7 . 8 Figure 7 .
Figure 7. Correction factor for modelled rainfall intensities fitted into higher order equations.

Figure 8
Figure8presents the future IDFs projections shown by box and whiskers plot for three quartile regions, which is compared with baseline (1981-2010) and Table5shows the detail changes in future intensities compared to baseline period.Comparison of present and future IDFs for return periods (2-, 5-and 20-years), showed that intensities increases under climate change at most of the rainfall durations.In comparison to the rainfall intensities the uncertainties among GCMs are found to be higher for smaller duration, which then increases with higher return periods.Even the median of the projected nine GCMs data showed the increases in future rainfall intensities when compared to the present IDFs, in all durations except for 1 h duration in 20 years return period.The maximum projected by some GCMs shows future rainfall intensity increments also in 1 h, while even minimum projected by some GCMs is higher than the present scenario, except at smaller durations in 20 years return period.There is no significant difference among the results from SRA1B and SRA2 scenarios in 2-and 5-years return periods, while SRA2 showed maximum increments in 20 years return period.It can be observed that few of the GCMs in 20 years return period at 1 h duration projected decrease in rainfall intensities.Compared to 1981-2010, during 2011-2030 in 20 years return period, the maximum decrease in 1 h rainfall intensities is projected by NCPCM in SRA1B (−7.8%) and GFCM21 in SRA2 (−10.1%).Similarly, for 2046-2065 in 20 years return period, the maximum decrease in 1 h rainfall intensities are projected by IMCM3 and HADGM3 in SRA1B (−10.2%) and HADGEM in SRA2 (−18.3%).The maximum increase in 1 h rainfall intensities during 2011-2030 are projected by NCCCSM in SRA1B (20.4%) and IMCM3 in SRA2 (20.4%) in 2 years return period, NCCCSM in SRA1B (13.9%) and INCM3 in SRA2 (18.7%) in 5 years return period, GFCM21 in SRA1B (3.7%) and NCCCSM in SRA2 (25.9%) in 20 years return period.The maximum increase in 1 h rainfall intensities during 2046-2065 are projected by NCCCSM in SRA1B (25.6%) and IMCM3 in SRA2 (15.7%) in 2 years return period, NCCCSM in SRA1B (18.7%) and HADCM3 in SRA2 (14.3%) in 5 years return period,

Table 1 .
Global climate models applied in the study.

Table 3 .
Goodness of Fit between different distributions for annual maximum precipitation, showing parameters, Kolgomorov-Smirnov (K-S) test and Chi-squared (χ2) test statistics and p-values, and intensities at different return periods for each distribution and different durations.

Table 4 .
Ratio and difference between observed and simulated rainfall intensities.

Table 4 .
Ratio and difference between observed and simulated rainfall intensities.