Development of Monsoonal Rainfall Intensity-Duration-Frequency ( IDF ) Relationship and Empirical Model for Data-Scarce Situations : The Case of the Central-Western Hills ( Panchase Region ) of Nepal

Intense monsoonal rain is one of the major triggering factors of floods and mass movements in Nepal that needs to be better understood in order to reduce human and economic losses and improve infrastructure planning and design. This phenomena is better understood through intensity-duration-frequency (IDF) relationships, which is a statistical method derived from historical rainfall data. In Nepal, the use of IDF for disaster management and project design is very limited. This study explored the rainfall variability and possibility to establish IDF relationships in data-scarce situations, such as in the Central-Western hills of Nepal, one of the highest rainfall zones of the country (~4500 mm annually), which was chosen for this study. Homogeneous daily rainfall series of 8 stations, available from the government’s meteorological department, were analyzed by grouping them into hydrological years. The monsoonal daily rainfall was disaggregated to hourly synthetic series in a stochastic environment. Utilizing the historical statistical characteristics of rainfall, a disaggregation model was parameterized and implemented in HyetosMinute, software that disaggregates daily rainfall to finer time resolution. With the help of recorded daily and disaggregated hourly rainfall, reference IDF scenarios were developed adopting the Gumbel frequency factor. A mathematical model [i = a(T)/b(d)] was parameterized to model the station-specific IDF utilizing the best-fitted probability distribution function (PDF) and evaluated utilizing the reference IDF. The test statistics revealed optimal adjustment of empirical IDF parameters, required for a better statistical fit of the data. The model was calibrated, adjusting the parameters by minimizing standard error of prediction; accordingly a station-specific empirical IDF model was developed. To regionalize the IDF for ungauged locations, regional frequency analysis (RFA) based on L-moments was implemented. The heterogeneous region was divided into two homogeneous sub-regions; accordingly, regional L-moment ratios and growth curves were evaluated. Utilizing the reasonably acceptable distribution function, the regional growth curve was developed. Together with the hourly mean (extreme) precipitation and other dynamic parameters, regional empirical IDF models were developed. The adopted approach to derive station-specific and regional empirical IDF models was statistically significant and useful for obtaining extreme rainfall intensities at the given station and ungauged Hydrology 2018, 5, 27; doi:10.3390/hydrology5020027 www.mdpi.com/journal/hydrology Hydrology 2018, 5, 27 2 of 27 locations. The analysis revealed that the region contains two distinct meteorological sub-regions highly variable in rain volume and intensity.


Background
Nepal has experienced numerous natural hazard events over the past resulting in enormous economic losses and thousands of causalities, compounding the country's already high level of poverty [1][2][3][4][5].The frequent natural and human-induced disasters are due to its fragile geomorphology, active tectonics, and unplanned human activity on steep slopes combined with climate extremes [3,4,[6][7][8].Flooding and landslides occur frequently, mostly triggered by extreme precipitation or by earthquakes, as evidenced by the 2015 Gorkha earthquake.The importance of human interventions, such as road construction, on mass movements was at best underestimated and largely neglected by researchers and authorities in Nepal [9].Landslides are caused by a number of underlying natural and human factors including slope aspect, gradient, soil type, land cover (changes), proximity to rivers and land use (i.e., road construction, quarrying) and most often triggered by rainfall or earthquakes [8,10].Rainfall, while critical to maintaining ecosystem services and supporting livelihoods, is thus a triggering factor that needs to be understood and managed in order to reduce impacts due to flooding and mass movements and erosion.Reducing and managing climate-induced disaster-risk (e.g., damage due to floods and landslides) relies on knowledge of the frequency and intensity of rainfall events [11].Extreme precipitation-induced pluvial flood and shallow landslide disaster management requires the establishment of intensity-duration-frequency (IDF) relationships of extreme events in order to formulate better design guidelines for the development of infrastructure so as to reduce impacts due to disaster-risk, and save lives, property and ecosystems [12].
However, inadequate rainfall data continue to hamper the establishment of reasonable IDF relationships [13][14][15].Researchers have realized that the design of water resources projects, the management of storm water runoff and disaster mitigation planning with scarce and insufficient data are always a challenge [13,16,17].The situation is even more challenging in the least-developed countries such as Nepal.In Nepal, the Department of Hydrology and Meteorology (DHM) is the government organization responsible for meteorological instrumentation and maintaining the climate variable database including rainfall.DHM has reported that there is no instrumentation for measuring fine-scale rainfall (e.g., hourly, sub-hourly precipitation) nor did it develop any systematic IDF relationships for the country.
This research explored the possibility of establishing an IDF relationship in the data-scarce environment of Nepal and developed an empirical IDF model to better estimate extreme rainfall events at some fixed duration of time and recurrence interval for the Panchase region in the Central-Western Hills of Nepal, a region with one of the highest annual rainfall amounts (between 4000-5000 mm, mean = 2984.40mm, standard deviation (SD) = 1497.60).The objectives of the study were to develop monsoon season IDF relationships for the low-resolution rainfall data recording situation followed by the development of location-specific empirical IDF models, and to demonstrate the application of IDF for ungauged locations.For this purpose, the recorded daily and disaggregated hourly rainfall series were the main data sources.Using the IDF, we can better understand the short-and long-term rainfall intensity during the monsoon season that triggers mass movements and other hazards induced by intense rainfall in this region.Furthermore, the IDF relationship can be used in water resources management and infrastructure planning and design.

Rationale
The importance of rainfall and IDF relationships in addressing extreme precipitation-induced hazards such as mass-movement has been emphasized by many researchers [18][19][20].IDF curves are also the graphical representation that summarizes the important statistical properties of extreme precipitation events [21].First established in 1932 by Bernard [17], IDF is a statistical relationship of the intensity, duration and frequency of rainfall derived from historical rainfall data [13,16,22].Since then, many different sets of relationships have been constructed for different parts of the world [13,16,17,[22][23][24][25].
In the literature there are several functions to establish IDF relationships (e.g., [16,22,26]).Chen [26] derived a generalized IDF relationship for any location in the United States using three basic rainfall depths (e.g., 1 h 10-year, 24 h 10-year, and 1 h 100-year rainfall depth).Baghirathan and Shaw [27] and Gert et al. [28] proposed regional IDF formulae for ungauged areas while Kothyari and Garde [29] used daily rain and two-year return period to establish the IDF relation in India.Furthermore, Koutsoyiannis et al. [22] proposed a mathematical approach to formulate IDF relationships followed by the construction of IDF curves utilizing point information of long-term historical rainfall time series in the context of geographical variability and the regionalization of the IDF model.
The concept of regional IDF relationships was executed by Yu and Chen [30] and Madsen et al. [31] who examined regression techniques, while Willems [32], Yu et al. [33] and Langousis and Veneziano [19] applied a scaling method and developed regional IDF curves.Dalrymple [34] proposed a regional frequency analysis (RFA) method for pooling various data samples, also known as the index-flood procedure in hydrology [35].Hosking et al. [36] studied the properties of a probability-weighted moments (PWMs) method based on L-moment, and Hosking and Wallis [11] showed the application L-moments for the RFA.L-moments are an efficient tool used to detect the homogeneous regions, to select suitable regional frequency distribution, and to predict extreme precipitation quantiles at a region of interest.The IDF relationship if regionalized can minimize computational time and effort to obtain the IDF curves for areas where rainfall gauging stations are not installed.
Nepal as a whole is dominated by S-E monsoon [46], where topography has a considerable effect on the rainfall patterns.There are four distinct hydrological seasons: pre-monsoon (April and May: AP), monsoon (June to September: JJAS), post-monsoon (October and November: ON) and winter or dry period (December to March: DJFM) in the country [46].The monsoonal rainfall is highly variable over the country (~1000 mm-~4500 mm) and is intense in nature where more than 80% of the annual rainfall occurs in the 4 months of the monsoon [18,47], resulting in landslides, debris flows, flooding and sedimentation, threatening livelihoods and properties in numerous ways [48].

The Panchase Region
Panchase is a mountainous region in the middle-hills of Central-Western Nepal between latitudes 28 • 12 N to 28 • 18 N and longitudes 83 • 45 E to 83 • 57 E. The region is located in three districts-Kaski, Syangja and Parbat (Figure 1).The elevation varies from 742 m above sea level (masl) (outlet of Phewa Lake, near Pokhara city) to 2523 masl (Panchase Peak) and is characterized by hot, humid summers and cool-temperate winter seasons [6].The Panchase hill range extends from the south-west to the north-west directions, dividing the region into two distinct eastern and western regions.In order to evaluate the rainfall variability and to establish the IDF relationship of the region, the historical rainfall data of 11 weather stations were collected from the DHM (Table 1).The region is one of the most studied areas in the country [48][49][50][51][52][53] due to the importance of Phewa Lake that promotes economic activities and biodiversity.However, the study of rainfall patterns, their intensity and impact have not been properly addressed yet.The Ecosystem Protecting Infrastructures and Communities (EPIC) project established and monitored three community-based rural roadside bioengineering demonstration sites in the region.In addition to capacity building and policy advocacy, the project combined formal and citizen-science by mobilizing local people and exploring local knowledge regarding climate extremes, plant species and the importance of rural access roads for building climate-resilient communities.In order to link local knowledge with formal science, the project installed three continuous recording meteorological stations in the period November-February 2014 next to each demonstration site.
The region is one of the most studied areas in the country [48][49][50][51][52][53] due to the importance of Phewa Lake that promotes economic activities and biodiversity.However, the study of rainfall patterns, their intensity and impact have not been properly addressed yet.The Ecosystem Protecting Infrastructures and Communities (EPIC) project established and monitored three community-based rural roadside bioengineering demonstration sites in the region.In addition to capacity building and policy advocacy, the project combined formal and citizen-science by mobilizing local people and exploring local knowledge regarding climate extremes, plant species and the importance of rural access roads for building climate-resilient communities.In order to link local knowledge with formal science, the project installed three continuous recording meteorological stations in the period November-February 2014 next to each demonstration site.

Data Quality and Rainfall Variability
Historical daily rainfall data of 11 weather stations (Table 1 and Figure 1) in and around Panchase region with a recording period of above 30 years were obtained from DHM.Only the station in Walling had a shorter recording period of 22 years.Frequency analysis of rainfall time series requires that the data are homogeneous and independent [54].In order to evaluate the data quality, we performed homogeneity test and evaluated for missing values.Among several methods for the homogeneity test (see [55]), we implemented simple cumulative deviations of the data from the mean according to Buishand [54] to check whether the sample data were from the same population (Equations ( 1)-( 3), Table 1).The cumulative deviations from the mean can be expressed as: where, X i are the records from the meteorological series 1 X , 2 X , …… n X and X is the mean.
0 k S  and k n S  , respectively, are the initial and final value of the series.
Buishand [54] demonstrated cumulative deviations of the mean of rainfall records are often rescaled dividing k S by the sample standard deviation ( ).By evaluating the maximum (Q) or the range (R) of the rescaled cumulative deviations from the mean, the homogeneity of the data series can be tested, where:  1 for summary).

Data Quality and Rainfall Variability
Historical daily rainfall data of 11 weather stations (Table 1 and Figure 1) in and around Panchase region with a recording period of above 30 years were obtained from DHM.Only the station in Walling had a shorter recording period of 22 years.Frequency analysis of rainfall time series requires that the data are homogeneous and independent [54].In order to evaluate the data quality, we performed homogeneity test and evaluated for missing values.Among several methods for the homogeneity test (see [55]), we implemented simple cumulative deviations of the data from the mean according to Buishand [54] to check whether the sample data were from the same population (Equations ( 1)-(3), Table 1).The cumulative deviations from the mean can be expressed as: where, X i are the records from the meteorological series X 1 , X 2 , . . . . . .X n and X is the mean.S k=0 and S k=n , respectively, are the initial and final value of the series.Buishand [54] demonstrated cumulative deviations of the mean of rainfall records are often rescaled dividing S k by the sample standard deviation (σ).By evaluating the maximum (Q) or the range (R) of the rescaled cumulative deviations from the mean, the homogeneity of the data series can be tested, where: According to Raes et al. [56] high values of Q and R indicate that data of the series are not from the same population and the fluctuations are not purely random.
The inhomogeneous data series were excluded from the analysis whereas homogeneous data series containing less than 5% missing values were further evaluated.To complete the data series, missing values were imputed using the RClimTool [57] and nearest neighborhood method as explained in XLSTAT [58].
To understand the rainfall variability in the Panchase region, we adopted a generic approach of analyzing the completed daily rainfall series in terms of the variation of annual monsoonal rain [JJAS], variation in monsoonal dry days, and extreme events of over 100 mm of rain depth in 24 h.To better estimate the variation, we grouped the rainfall series according to the hydrological year that starts from April and focused on monsoonal rainfall depth.The evaluation of historical monsoonal dry/wet days can be an indicator to understand rainfall variability, as shown by other researchers (e.g., [59]).

Disaggregation of Daily Rainfall
A general framework to generate synthetic rainfall time series at finer time resolution consistent with the given coarser resolution, preserving the statistical characteristics (e.g., mean, standard deviation, lag 1 auto-correlation and percentage of dry days) of both scales, was initially proposed by Koutsoyiannis [42] and Koutsoyiannis and Manetas [60].Koutsoyiannis and Onof [44] extended the methodology to disaggregate daily rainfall into hourly data by coupling with the Bartlett-Lewis rectangular pulse (BLRP) model [41,44] with adjusting procedure [60].The original BLRP model proposed by Rodriguez-Iturbe et al. [61] consists of five parameters (λ, β, γ, η, and µ x ) and modelled the rainfall using rectangular pulses and characterized it by the parameters (Figure 2).The general assumption of the BLRP model as proposed by Rodriguez-Iturbe et al. [61] are:

•
The occurrence of random storm events (t i ) is assumed to be modelled as a Poisson process with rate λ and each event i is associated with a random number of cells.

•
Each storm event t ij , is assumed as a precipitation rectangular pulse with random duration t d and the origin of storm events t ij of each cell j occurs following a second Poisson process with rate β.The inter-arrival time of two subsequent storm events (i.e., successive cells) is independent, identically distributed and follows an exponential distribution.

•
The cell-generation process terminates after time span of ν i following the exponential distribution rate γ.Also, the number of cells per storm contains a geometric distribution of mean µ c = 1 + β/γ.

•
The random precipitation rectangular pulse duration t d is modelled as w ij and also follows exponential distribution with rate η.• Finally, the cell intensity x ij is assumed to be exponentially distributed with mean µ x .
According to Rodriguez-Iturbe et al. and Onof and Wheater [62,63] the BLRP model reproduced the basic statistics but they observed difficulties in reproducing the temporal characteristics.According to Kossieris et al. [41], in order to improve the model's flexibility in generating a greater diversity of rainfalls, Rodriguez-Iturbe et al. [63] modified the original model so that the parameter η is randomly varied from storm to storm according to gamma distribution with a shape parameter α and rate parameter ν in such way that the ratios of cell origin rate (β) and storm duration rate (γ) to the parameter η (i.e., κ = β/η and φ = γ/η) are kept constant.Also, the parameters β and γ are random variables that follow a gamma distribution with common shape parameter α and rate parameters ν/κ and ν/φ [41].Implementing the modified approach, the model consists of six parameters (λ, α, ν, κ, φ and µ x ), what we called modified BLRP (MBLRP).and rate parameters   and   [41].Implementing the modified approach, the model consists of six parameters (  ,  ,  ,  ,  and x  ), what we called modified BLRP (MBLRP).
Figure 2. Graphical representation of the Bartlett-Lewis rectangular pulse (BLRP) model according to [61].Filled and open circles denote, respectively, the storm origins and cell arrivals [41].
The successful application of the MBLRP model in different climatic regions was reported by several researchers (e.g., [38,39,41,44,62,[64][65][66][67][68][69][70][71].The model was successfully verified and implemented through HyetosMinute, and R-based computer software by Kossieris et al. [41].According to Kossieris et al. [41], HyetosMinute is an extended and improved version of the Hyetos model [43] to disaggregate daily rainfall depth to sub-hourly resolution [72].HyetosMinute implements the MBLRP scheme of the Poisson-clusters model and disaggregates daily rainfall, dividing the rainfall depth into many clusters as sequences of wet days separated by at least one dry day.Also, each sequence is considered to be an independent event [39,41]. We used the HyetosMinute [72] to generate the required level of rainfall resolution according to Kossieris et al. [41].As far as we are aware, this is the first attempt to apply a rainfall disaggregation approach in order to establish IDF relationships in Nepal and, thus, we evaluated the model parameter for each station first for the annual and second for the seasonal rainfall series.We focused on available daily monsoonal rainfall series to parameterize the model adopting a global optimization algorithm.The algorithm was implemented in a computer package called Evolutionary Annealing-Simplex (EAS) in R software [41,72] utilizing the historical statistical characteristics in terms of mean, variance, covariance, and percentage of dry days.While performing the EAS, suitable boundary conditions of the model parameters were estimated executing iterations, [41] in such a way that the model preserves the statistical characteristics of the historical rainfall.
The global optimization algorithm generated statistical parameters (  ,  ,  ,  ,  and x  ) derived from the daily available monsoonal rainfall data, and rainfall depth as an input disaggregation was executed in HyetosMinute [41,72].

Evaluation of Disaggregation Model
In the case of Nepal, we are facing an inadequate data situation; hence we are limited in our ability to tune the model.We performed the model evaluation using limited fine-resolution (hourly) rainfall data recorded within the study region.We implemented the model performance evaluation scheme according to Kossieris et al. [41].The EPIC project has established three tipping-bucket type weather stations within the test area which were calibrated to measure rainfall volume of 0.2 mm per tip with temporal resolution of an hour.Hourly rainfall data (JJAS 2016) of Gharelu, one of the test sites of EPIC, were obtained and used to evaluate the disaggregation model performance.For this purpose, we first aggregated the hourly rainfall to daily series then the series was disaggregated adopting the procedure discussed in Section 2.2.The statistical characteristics of the aggregated daily series was calculated and fitted into the EAS model to estimate the parameters (Table 2).[61].Filled and open circles denote, respectively, the storm origins and cell arrivals [41].
The successful application of the MBLRP model in different climatic regions was reported by several researchers (e.g., [38,39,41,44,62,[64][65][66][67][68][69][70][71].The model was successfully verified and implemented through HyetosMinute, and R-based computer software by Kossieris et al. [41].According to Kossieris et al. [41], HyetosMinute is an extended and improved version of the Hyetos model [43] to disaggregate daily rainfall depth to sub-hourly resolution [72].HyetosMinute implements the MBLRP scheme of the Poisson-clusters model and disaggregates daily rainfall, dividing the rainfall depth into many clusters as sequences of wet days separated by at least one dry day.Also, each sequence is considered to be an independent event [39,41]. We used the HyetosMinute [72] to generate the required level of rainfall resolution according to Kossieris et al. [41].As far as we are aware, this is the first attempt to apply a rainfall disaggregation approach in order to establish IDF relationships in Nepal and, thus, we evaluated the model parameter for each station first for the annual and second for the seasonal rainfall series.We focused on available daily monsoonal rainfall series to parameterize the model adopting a global optimization algorithm.The algorithm was implemented in a computer package called Evolutionary Annealing-Simplex (EAS) in R software [41,72] utilizing the historical statistical characteristics in terms of mean, variance, covariance, and percentage of dry days.While performing the EAS, suitable boundary conditions of the model parameters were estimated executing iterations, [41] in such a way that the model preserves the statistical characteristics of the historical rainfall.
The global optimization algorithm generated statistical parameters (λ, α, ν, κ, φ and µ x ) derived from the daily available monsoonal rainfall data, and rainfall depth as an input disaggregation was executed in HyetosMinute [41,72].

Evaluation of Disaggregation Model
In the case of Nepal, we are facing an inadequate data situation; hence we are limited in our ability to tune the model.We performed the model evaluation using limited fine-resolution (hourly) rainfall data recorded within the study region.We implemented the model performance evaluation scheme according to Kossieris et al. [41].The EPIC project has established three tipping-bucket type weather stations within the test area which were calibrated to measure rainfall volume of 0.2 mm per tip with temporal resolution of an hour.Hourly rainfall data (JJAS 2016) of Gharelu, one of the test sites of EPIC, were obtained and used to evaluate the disaggregation model performance.For this purpose, we first aggregated the hourly rainfall to daily series then the series was disaggregated adopting the procedure discussed in Section 2.2.The statistical characteristics of the aggregated daily series was calculated and fitted into the EAS model to estimate the parameters (Table 2).We calculated the mean (E n ), variance (Var n ) and skewness (Skew n ) of the disaggregated series and compared the characteristic of aggregated series for each of the monsoonal months.Furthermore, we evaluated the recorded and disaggregated hourly intensity for the test data series for which a goodness-of-fit (GoF) test was performed.

Fitting of Probability Distribution Function and Construction of Reference Intensity-Duration-Frequency (IDF)
The probability distribution of time series with daily rainfall data is important in the field of hydrology and meteorology [73].It is also important for the construction of IDF curves that requires the fitting of a probability distribution function (PDF) according to Koutsoyiannis et al. [22].Several distribution models (e.g., Gumbel Extreme Value Type I (EV I), Generalized Extreme Value (GEV), Log Pearson Type III, Beta, Gamma, log-normal, Normal, etc.) are particularly useful for hydrological and meteorological time-series data analysis [22,[74][75][76][77].In Nepal there is currently limited use of IDF curves/model, and thus no preferred distribution model to be fitted for the rainfall data series.However, in some flood-frequency studies, the Gumbel Extreme Value Type I and Log-Pearson Type III were used [52].In this research, the PDF was evaluated and the best fit was chosen for which EasyFit statistical software developed by MathWave-Technology [78] was used as employed in Gamage et al. [79] and Misic [80] and as discussed in Hosking et al. [81], Stedinger et al. [82], Koutsoyiannis [83,84] and Millington et al. [75].The homogeneous daily rainfall data of eight weather stations discussed earlier (see Table 1) were examined by fitting the PDF, in particular GEV and EV I (Equations ( 4) and ( 5)).
where, κ > 0, λ > 0 and ψ are shape, scale and location parameters, respectively.For κ = 0, the GEV distribution turns into the Gumbel distribution [82]: While analyzing the daily time series we recognized that the data contain some degree of seasonality effects, leading to separation of the annual series into seasonal to better understand the rainfall characteristics.The PDF for the monsoonal daily and annual extreme time series rainfall was assessed using the test statistics adopting the Kolmogorov-Smirnov (K-S), Anderson-Darling (A-D) and chi-squared test (χ 2 ) as illustrated in the EasyFit software [78].Accordingly, the best-fitted PDF was chosen and model parameters (κ, λ and ψ) were estimated.
We proposed the construction of reference IDF utilizing the recorded monsoonal daily and disaggregated hourly extreme rainfall and computed the frequency of precipitation depth P T , (in mm) for the given rainfall duration t d (in hour) with specified return period T r (in Years) according to Wilson [85]; where, P T is frequency of precipitation depth, P ave is mean rainfall of monsoonal annual series (recorded daily and disaggregated hourly) extremes, S d is standard deviation of annual series, and K T is Gumbel frequency factor for the given duration (Equation ( 8)): where, T r is return period and 0.5772 is Euler's constant.
From Equations ( 6)-( 8), we determined the rainfall intensities for a given duration and given return period, what we called IDF scenarios, representing given locations (i.e., stations).These IDF scenarios were later used as a reference to evaluate station-specific fine resolution IDF curves developed by implementing the mathematical model proposed by Koutsoyiannis et al. [22] expressed as below: According to Koutsoyiannis et al. [22], the above expression (Equation ( 9)) can be expressed in terms of a(T) and b(d) as below (Equations ( 10) and ( 11)): where, λ and ψ are scale and location parameters of the distribution, respectively, and T r is return period for the best-fitted PDF (i.e., EV I, see Section 3.4).
In order to evaluate the PDF, three sets of parameters were derived from three different levels of analysis, respectively, from annual daily extremes, hourly extremes and monsoonal time-series rainfall depth.The main objective was to establish monsoonal IDF relationship; we focused on the monsoonal rainfall series and thus the PDF parameters.For this reason, the best-fitted PDF was adopted and we obtained λ and ψ as discussed in Section 2.4.The empirical parameter θ and η were estimated adopting Koutsoyiannis et al. [22].Accordingly, separate IDF curves for fine time resolution were developed representing each station.In order to estimate the performance of the mathematical model, the reference IDF curves (t d = 1 and 24 h) were fitted to the mathematical model and results were evaluated using the standard error of prediction (SEP) (Equation ( 12)) and the GoF test.

SEP =
where, SEP is the standard error of prediction, I s and I p are, respectively, the intensities (at given duration t d and return period T r ) observed in the reference IDF and mathematical model predicted intensities, and n is the number prediction (i.e., return periods).
In order to improve the prediction and better fitting of reference IDF, the mathematical model was calibrated, adjusting the parameter by minimizing the SEP implemented in SOLVER add-ins available in MS-EXCEL.This process of analysis and adjustment of the parameters led to the development of individual station-specific empirical models.

Estimation of IDF for Ungauged Location
Researchers have demonstrated several ways to estimate extreme storms in terms of IDF for ungauged locations.Among them logistic and multiple linear regression, spatial interpolation techniques and regional frequency analysis (RFA) are particularly in use (e.g., [11,16,22,86]).The RFA technique was initially introduced for extreme flood estimation, known as the index-flood method based on L-moments, and is also in use to estimate extreme rainfall events in the homogeneous region [11,35,82,87].In this study, the L-moment based RFA method was used to establish regional IDF relationship and develop the regional empirical IDF model.
Details about the L-moments method can be found in Hosking and Wallis [11].In brief, the L-moments are defined as the linear function of the probability-weighted moments (PWMs) explored by Greenwood et al. [88], with the benefit of offering a description of the shape of a probability distribution by L-skewness and L-kurtosis.The PWM can be defined as (Equation ( 13)): The rth L-moment λ r is related to the rth PWM, according Hosking [89]: Therefore, the first four L-moments are: The L-moments are independent of units of measurement, called L-moment ratios, and are defined to the quantities, according Hosking [89]: where, τ is L-coefficient of variations (L-C v ), τ 3 is L-coefficient of skewness (L-C s ) and τ 4 is the L-coefficient of kurtosis (L-C k ).According to Hosking and Wallis [11], the L-moments can be defined uniquely and no two distributions can have the same L-moments.
L-moments have superior abilities to conventional moments in discriminating between different distributions, because the L-moment ratio estimators of location, scale and shape are nearly unbiased, regardless of the probability distribution from which the observations arise and efficient estimators of the characteristic of climatic data and of the parameters of the distribution.The L-moment based RFA consist of the following steps:

•
Screening of data through discordancy measure; • Identification of homogeneous regions; • Selection of regional distribution and goodness of fit measure; • Estimation of regional growth curve using index-flood procedure.

Screening of Data and Discordancy Measure
Screening of data comprises evaluation of the data in terms of gross error, homogeneity, consistency and stationarity.In addition, according to the assumptions of the index-flood method, the serial and spatial independency of sites is required.For the detection of data quality, Hosking and Wallis [11] proposed evaluation of discordancy measures D i as an indicator for a site i.Mathematically, D i can be expressed as (Equation ( 15)): where, µ i = t (i) , t T are vectors where t, t 3 and t 4 are L-moment ratios µ is mean of µ i .The measure D i indicates how far µ i is from the center of the region relative to the size of the region.Hosking and Wallis [11] recommended that the site is discordance if D i > 3.

Identification of Homogeneous Region
This is the core part of RFA and based on the grouping of sites to regions in which the samples are from the same population and consist of the same frequency distributions.We implemented the L-moment based RFA by dividing the eight rain gauge stations distributed into the heterogeneous region to two homogeneous sub-regions (eastern region: Lumle, Bhadaure, Pokhara-Airport and Khairenitar; western region: Kusma, Karki-Neta, Syangja and Walling).The division was made considering the meteorological and physical characteristics such as mean annual precipitation, altitude, latitude and distance to the gauging stations to Lake Phewa.According to Hosking and Wallis [11], heterogeneity measure estimates the degree of heterogeneity in a group of sites and assesses whether they might reasonably be treated as a homogeneous region.The heterogeneity measure (Equation ( 16)) compares the observed and simulated dispersion of L-moments for N sites under consideration.
The regions are regarded as "acceptably homogeneous" when H j < 1, "possibly heterogeneous" when 1 < H j < 2, and "definitely heterogeneous" when H j > 2. The details of the calculation of H j are given in Hosking and Wallis [11].

Selection of Regional Distribution and Goodness of Fit
L-moment ratio diagrams were constructed using the unbiased estimators of L-moments according to Hosking [89].The curves show the theoretical relationships between L-Skewness (L-C s ) and L-Kurtosis (L-C k ) of various candidate distributions.In addition, Z-statistics (Z DIST ) defined by Hosking and Wallis [90] compare simulated L-C s and L-C k of fitted distribution with the regional average L-C s and L-C k values obtained from observed data.The following relation (Equation ( 17)) defines the Z DIST : where, τ DIST   4   is the L-C k of fitted distribution, and β 4 and σ 4 simulated regional bias and simulated regional standard deviation of t R 4 .The simulation was made with the fitted Kappa distribution to regional L-moments.The fit is regarded as adequate if Z DIST is close to zero and acceptable if Z DIST ≤ 1.64 at a confidence level of 90%.
The satisfactory distributions were obtained as depicted in the Z DIST and L-moment ratio diagram.Among the satisfactory distributions, simpler distributions were chosen, and accordingly regional growth curves were developed.

Estimation of Regional Growth Curves using Index-Flood Procedure
The important assumption of the index-flood or the growth curves procedure was that the frequency distributions of N stations in the homogeneous region are identical apart from a site-specific scaling factor of the ith site which can be written as (Equation ( 18)): where, Q i (F) is the precipitation at site i at a given return period, µ i is the index-flood [11], N is the number of sites, and q(F) is the regional growth curve, a dimensionless quantile function common to the homogeneous site.The index-flood µ i also known as the scaling factor is the sample means of the data at site i.Since, we were concerned with the short duration, extreme precipitation and thus hourly (extreme) mean precipitation was considered as the index-flood (µ i = µ i(hr) ).The following equations (Equations ( 19) and ( 20)) presents the regional growth curve defined by the quantile function of chosen distribution (Gen.Logistic (GLO) and GEV) for the western and eastern regions under consideration: where, λ, ψ and κ are, respectively, the scale, location and shape parameter and T r is the return period.
Utilizing the above two regional quantile functions together with the µ i , we developed a nominator [a(T)] of Equation (9).For the denominator [b(d) = (t d + θ) η ] of Equation ( 9), we made a use of the reference IDF and solved the following (Equation ( 21)) relation to obtain two parameters (θ and η) of b(d) at a given duration (t d = 1 and 24 h) by implementing SEP (Equation ( 12)) as discussed in Section 2.4:

Data Quality and Rainfall Variability
The evaluation of the historical daily rainfall series from the available 11 stations led to the exclusion of three stations (Chapakot, Nr. 810; Lamachaur, Nr. 818 and Ghandruk, Nr. 821) from further analysis as these data passed the threshold of inhomogeneity or contained more than five percent of missing values.Therefore, only eight stations were considered for further analysis.Table 1 presents the characteristics of the rainfall records, with missing values and homogeneity test results.
The daily rainfall series of the eight stations showed that 81.40% of rainfall occurred during the monsoon season followed by pre-monsoon (11.10%), winter (4.0%) and post-monsoon (3.50%).Out of the eight stations, three are in the eastern part (Bhadaure-Deurali, Pokhara-Airport and Khairenitar), four in the western part (Syangja, Walling, Karki-Neta and Kusma) and Lumle on the N-W part of the Panchase hill range.The mean monsoonal rain of the eastern and western part was 2657 mm and 2027 mm, respectively.The analysis also revealed that over a 30-year period in the eight stations the total number of storms exceeding 100 mm in 24 h was 998 of which 57 exceeded 200 mm (Figure 3b).The highest recorded annual total precipitation in the region was 5631 mm in 1984 in Lumle with the mean monsoonal sum of 4681 mm (Figure 3a) and the recorded maximum rain in 24 h was 357 mm at Pokhara-Airport.This analysis clearly indicated that the rainfall in Panchase region is highly variable and the hill range distinctly divided the area into two meteorological regions.

Disaggregation of Daily Rainfall Depth
The rationale behind the rainfall disaggregation was to split the daily rainfall into hourly to subhourly [13,39,41,42,44,71,91] when the resolution of the available data is not fine enough.We derived MBLRP model parameters (Table 3) utilizing the statistical characteristics of historical monsoonal daily rainfall series for all eight stations and implemented the HyetosMinute software in R. We observed that the higher skewness value of the daily rainfall series was possibly due to strong variability in the daily rainfall time series, whereas the annual rainfall time series is rather smoothened out leading to low skewness, which was expected.The HyetosMinute effectively utilized the MBLRP model and disaggregated the monsoonal daily rainfall to synthetic hourly series even for long clusters of wet days separated by at least 1 dry day.

Evaluation of Disaggregation Model
The model discussed in Section 2.2 was evaluated utilizing the short duration hourly rainfall depth (June-September, 2016) available from the EPIC project.For this purpose, we implemented the global optimization algorithm and obtained the parameters discussed in Sections 2.3 and 2.4 according to Kossieris et al. [41] by comparing the original and disaggregated rainfall series in terms of statistical characteristics (e.g., mean, variance and skewness).We summarized the results for each month demonstrating that the disaggregated series preserved the statistical characteristics of the recorded rainfall series (Table 4) for the continuous weather station of the EPIC project at Gharelu in Kaski.However, we noticed the hourly synthetic rainfall series contained small differences in variance whereas the mean and skewness were similar to that of the original daily series.

Disaggregation of Daily Rainfall Depth
The rationale behind the rainfall disaggregation was to split the daily rainfall into hourly to sub-hourly [13,39,41,42,44,71,91] when the resolution of the available data is not fine enough.We derived MBLRP model parameters (Table 3) utilizing the statistical characteristics of historical monsoonal daily rainfall series for all eight stations and implemented the HyetosMinute software in R. We observed that the higher skewness value of the daily rainfall series was possibly due to strong variability in the daily rainfall time series, whereas the annual rainfall time series is rather smoothened out leading to low skewness, which was expected.The HyetosMinute effectively utilized the MBLRP model and disaggregated the monsoonal daily rainfall to synthetic hourly series even for long clusters of wet days separated by at least 1 dry day.

Evaluation of Disaggregation Model
The model discussed in Section 2.2 was evaluated utilizing the short duration hourly rainfall depth (June-September, 2016) available from the EPIC project.For this purpose, we implemented the global optimization algorithm and obtained the parameters discussed in Sections 2.3 and 2.4 according to Kossieris et al. [41] by comparing the original and disaggregated rainfall series in terms of statistical characteristics (e.g., mean, variance and skewness).We summarized the results for each month demonstrating that the disaggregated series preserved the statistical characteristics of the recorded rainfall series (Table 4) for the continuous weather station of the EPIC project at Gharelu in Kaski.However, we noticed the hourly synthetic rainfall series contained small differences in variance whereas the mean and skewness were similar to that of the original daily series.In order to evaluate how well the model disaggregated the extreme intensity of the synthetic hourly rainfall depth, we compared the recorded and disaggregated extreme intensity for each of the monsoonal months.The GoF test showed that there was a good agreement among the recorded and synthetic hourly intensities of the rainfall series (χ 2 = 0.932, alpha = 0.05 and degree of freedom = 3).

Selection of Probability Distribution Function (PDF) and Parameter Estimation
The chosen PDF (EV I and GEV) were evaluated for which K-S, A-D and the chi-squared test were performed (Table 5).The K-S and A-D statistics ranked the three parameter (κ, ψ and λ) GEV distribution model better over the EV I, whereas the chi-squared ranked the two parameter (ψ and λ) EV I better than the GEV.Stedinger et al. [82] expressed that the EV I distribution is obtained when κ = 0, and the general shape of GEV turns to EV I where κ < 0.3.In our case, for the monsoonal annual extremes the κ value was always less than 0.3 (Table 5) whereas this was not fully satisfied for the daily series.Also, according to Koutsoyiannis et al. [22] the performance of GEV is not very satisfactory with a small number of samples as in our case.Based on the analysis at various levels (e.g., annual and monsoonal daily series, annual extremes and disaggregated hourly extremes series), we concluded that in this case the EV I distribution was better than GEV.

Construction of Reference and Empirical IDF Relationship
The reference IDF curves were developed for the coarser time resolution of 1 and 24 h duration utilizing the extreme rainfall from the recorded daily and disaggregated hourly monsoonal rainfall series.The frequency of extreme rainfall depth (P T ) and dimensionless Gumbel frequency factor (K T ) also known as Gumbel growth curve [92,93] were derived for the eight stations.The results showed that the frequency of precipitation is highest in the Bhadaure-Deurali station and the lowest in Khairenitar.Utilizing the method of Wilson [85], we developed the IDF curves for 1 and 24 h durations, which we called the reference IDF curves.These curves were the means to verify the finer time resolution IDF relationship.
In order to establish the IDF relationship for finer time resolution, the mathematical model proposed by Koutsoyiannis et al. [22] was parameterized.Accordingly, the best-fitted PDF parameters scale (λ) and location (ψ) and other empirical constants θ and η were estimated to be fitted into the mathematical model discussed in Section 2.4.The estimated parameters are shown in Table 6 which were used to develop the station-specific mathematical IDF relationship.

Evaluation of the IDF Relationship and Development of Empirical Model
The mathematically computed IDF relationship (Section 3.5) was evaluated by fitting reference IDF scenarios, and SEP and chi-square tests were performed.The test statistics indicated that the mathematical model overestimated the IDF values for the Lumle (Nr.814) station for both scenarios (i.e., 1 and 24 h).Similarly the Walling (Nr.826) and Bhadaure-Deurali (Nr.813) stations were poorly estimated for 1 h duration.The possible reason for this could be either due to the PDF that we chose or the limitation of the mathematical model where intense and prolonged rainfall occurs.This leads to the calibration of the model for which we again performed minimization of SEP by adjusting all four parameters (λ, ψ, θ and η) at a time.This reconstruction of the IDF significantly reduces the SEP and achieved significant test statistics except in Khairenitar station, for which both SEP and GoF were increased but yet significant.Table 7 compared the test statistics before and after the calibration of the mathematical model.Utilizing the calibrated empirical constant (λ , ψ , θ , η ) and fitting them into Equations ( 9)-( 11), we developed the station-specific empirical IDF model.
Interpretation of the IDF relationship indicated that the eastern part generally received higher rain than the western part.Variability was also observed in the rainfall intensity from east to west. Figure 4 shows computed IDF relationship of four stations (Lumle, Pokhara-Airport, Syangja and Kusma) from the equations shown in Table 8.The basis of this IDF relationship was the theoretical probabilistic approach of the rainfall data distribution to compute the parameters of a(T) and the optimization procedure for the parameters of b(d) that represent the dynamics of the rainfall pattern discussed in Section 2.4 and demonstrated in Koutsoyiannis et al. [22].However, there was some adjustment on the empirical constants that demonstrated better statistical significance and let the IDF curves pass through or much closer to the reference IDF point at given td and Tr.The reason for such an adjustment was also explained in Koutsoyiannis et al. [22] and Van de Vyver and Demaree [17] who stated that empirical considerations are not always consistent with the theoretical probabilistic approach of the IDF relationship.

Reganalization of IDF for Ungauged Locations
We executed RFA based on L-moments for the heterogeneous regions divided into two homogeneous regions (eastern region: Lumle, Pokhara-Airport, Bhadaure-Deurali and Khairenitar and western region: Kusma, Karki-Neta, Syangja and Walling) for which the heterogeneity measure was evaluated.The heterogeneity measure of the L-moment demonstrated that the sub-regions were found to be statistically homogeneous.The estimated heterogeneity measures (H) for the eastern and  The basis of this IDF relationship was the theoretical probabilistic approach of the rainfall data distribution to compute the parameters of a(T) and the optimization procedure for the parameters of b(d) that represent the dynamics of the rainfall pattern discussed in Section 2.4 and demonstrated in Koutsoyiannis et al. [22].However, there was some adjustment on the empirical constants that demonstrated better statistical significance and let the IDF curves pass through or much closer to the reference IDF point at given t d and T r .The reason for such an adjustment was also explained in Koutsoyiannis et al. [22] and Van de Vyver and Demaree [17] who stated that empirical considerations are not always consistent with the theoretical probabilistic approach of the IDF relationship.

Reganalization of IDF for Ungauged Locations
We executed RFA based on L-moments for the heterogeneous regions divided into two homogeneous regions (eastern region: Lumle, Pokhara-Airport, Bhadaure-Deurali and Khairenitar and western region: Kusma, Karki-Neta, Syangja and Walling) for which the heterogeneity measure was evaluated.The heterogeneity measure of the L-moment demonstrated that the sub-regions were found to be statistically homogeneous.The estimated heterogeneity measures (H) for the eastern and western regions were −0.97 and 0.34 respectively.Similarly, the L-moment ratio diagram and Z DIST demonstrated reasonably acceptable regional distribution functions (Table 9 and Figure 5).For the development of the regional growth curve, we chose one set up based on simple but statistically acceptable distributions, (1) Gen. Extreme Value (GEV) and ( 2) Gen. Logistic (GLO) respectively for the eastern and western region.Utilizing the empirical parameters (Table 10), the quantile function (Equations ( 19) and ( 20)) of the chosen distribution, and hourly (extreme) mean precipitation (see Table 3) as an index-flood, we established the IDF relation for an hour (t d = 1 h) and a given return period (T r = 2, 5, 10, 25, 50, 100, 200 years).
Hydrology 2018, 5, x 18 of 26 western regions were −0.97 and 0.34 respectively.Similarly, the L-moment ratio diagram and Z DIST demonstrated reasonably acceptable regional distribution functions (Table 9 and Figure 5).For the development of the regional growth curve, we chose one set up based on simple but statistically acceptable distributions, (1) Gen. Extreme Value (GEV) and ( 2) Gen. Logistic (GLO) respectively for the eastern and western region.Utilizing the empirical parameters (Table 10), the quantile function (Equations ( 19) and ( 20)) of the chosen distribution, and hourly (extreme) mean precipitation (see Table 3 21) and from the regional distribution we obtained the regional parameters shown in Table 10.We then evaluated the regional IDF with the station-specific reference IDF for the given duration (td = 1 and 24 h) and return period (Tr = 2, 5, 10, 25, 50, 100 and 200 years) for which SEP and GOF were performed.The test depicted that the stations Khairenitar (Nr.815), Lumle (Nr.814) and Walling (Nr.826) demonstrated insignificant IDF values (i.e., the regional model either overestimated or underestimated the IDF at the given duration (td = 1 h and 24 h)), indicating the existence of some degree of heterogeneity, leading to the adjustment/calibration in the empirical parameters in Table 10.In Table 11, we present the test statistics before and after the adjustment.By solving Equation (21) and from the regional distribution we obtained the regional parameters shown in Table 10.We then evaluated the regional IDF with the station-specific reference IDF for the given duration (t d = 1 and 24 h) and return period (T r = 2, 5, 10, 25, 50, 100 and 200 years) for which SEP and GOF were performed.The test depicted that the stations Khairenitar (Nr.815), Lumle (Nr.814) and Walling (Nr.826) demonstrated insignificant IDF values (i.e., the regional model either overestimated or underestimated the IDF at the given duration (t d = 1 h and 24 h)), indicating the existence of some degree of heterogeneity, leading to the adjustment/calibration in the empirical parameters in Table 10.In Table 11, we present the test statistics before and after the adjustment.The test statistics showed that the two homogeneous sub-regions (eastern and western) can be divided into further sub-regions such as three sub-regions in the eastern region [eastern sub-region 1: Khairenitar, eastern sub-region 2: Pokhara-Airport and Bhadaure-Deurali and eastern sub-region 3: Lumle] and two in the western region [western sub-region 1: Kusma, Karki-Neta and Syangja and western sub-region and 2: Walling].Utilizing the best-fitted regional distribution and parameters, we developed regional IDF formula for the region, as shown in Table 12.The adjustment of the empirical parameters helps to improve the test statistics except for Walling, indicating that the distribution we chose was not feasible for the western sub-region2.However, we could still say that the regional empirical model is useful in evaluating the IDF of rainfall for other sub-regions.

Rainfall Intensity in the Region
The computed rainfall intensity of eastern sub-region-3 was found to be the highest in terms of shorter duration rainfall (i.e., t d = 0.5, 1 and 2 h), whereas for the longer duration (i.e., t d = 24 h) the rainfall was relatively intense in western sub-region 2 followed by eastern sub-region 2 and western sub-region 2. Overall, rainfall in eastern sub-region 2 was less intense followed by western sub-region 2 for the shorter duration and return periods.Table 13 presents the computed rainfall intensity of some important return periods (T r ) and duration (t d ) for example.

Discussion
Geographically, the Panchase region can be divided into two sub-regions, eastern and western parts, where the Panchase hill range is sitting in the midway extending N-W to S-E direction.This study demonstrated that the region can also divide into two parts in terms of rainfall variability.Monsoonal rainfall over the region is highly variable where the eastern part (i.e., Kaski District) received higher rain (>3000 mm, except Kahirenitar) than the western part (i.e., Syangja and Parbat Districts), where mean monsoonal rainfall was below 2300 mm.The geographical setting of the eastern part is relatively wider in comparison to the western part, which is a valley containing several lakes (e.g., Phewa Lake, Begnas Lake and Rupa Lake, etc.).Rainfall is a complex process and the complexity is compounded due to the presence of lakes and high mountain topography.This complex process requires more detailed analysis considering wind direction; variation in daily temperature, solar radiation, etc. in order to better understand the rainfall variability.
Use of the historical characteristics of the recorded daily rainfall to parametrize the disaggregation model in order to generate finer resolution synthetic rainfall series is an important step forward for data-scarce regions [15,39,41,45].The method is useful especially for those locations where fine time-resolution rainfall data is not available, leading to better estimation of the finer resolution rainfall up to a minute [41].However, we performed the disaggregation procedure only for time periods of one hour or more since the possible error accumulation for finer resolution disaggregation is yet to be known, and has to be investigated (Kossieris, 2017:personal communication).
Although many studies (e.g., [15,37,[39][40][41][42]44,45,64,94,95] have reported on the generation of synthetic rainfall series and their effectiveness, we noticed that the attention is less on the intensities.We observed that the synthetic rainfall intensities were unable to represent the recorded data of a particular day since the model implements a Poisson process and assumed that the intensity varies exponentially.In our case, however, the statistics suggest that the synthetic intensities were statistically significant for the monsoonal months and useful for developing reference IDF scenarios for the study region.In order to better understand the synthetic rainfall intensities, a rigorous analysis is required where long-term fine-resolution time-series rainfall data are available.
The application of the Gumbel frequency factor technique is popular for evaluating frequencies of rain and flood storms [92,93,96,97].The technique is also frequently used to establish IDF relationships [22,93] using fine-resolution rainfall data [22].In our case, due to inadequate data resolution we demonstrated the use of synthetic rainfall intensities and developed reference IDF curves utilizing the technique of Gumbel frequency factor as scenarios to be fitted into the mathematically-derived empirical IDF model for better performance.In the literature, the use of scenario IDF computed from the disaggregated and recorded rainfall extremes to establish IDF relationship is limited.This is also the case for the approach of parameter adjustment in developing IDF relationships.
The L-moment based RFA method implemented in this study was able to demonstrate that the study area is heterogeneous in terms of geography and meteorology, which we observed while constructing the station-specific IDF through Gumbel.This has led to the division of the study area into sub-regions for the RFA.Application of L-moment ratios and other statistics in identifying reasonably acceptable distributions of the data set is robust, as was depicted in the results.Therefore, the technique is useful to better understand the IDF at the ungauged locations and can save computational time and resources.In this study, we observed that while fitting the regional IDF to the reference IDF, the chosen distribution slightly underestimated the intensity for the shorter duration and return periods.However, in most cases the result was statistically significant and useful.
The present study indicated that the Panchase region is highly variable in terms of rainfall amount and intensity.Short duration and intense monsoonal rain may have different effects on floods and landslides than less intense but extended duration monsoonal rainfall, respectively noticed in the western and eastern sub-regions of the area.
This research was carried out with limited data sources in terms of resolution, length of rainfall records and sparse observational network, which is common in the country.The low-resolution rainfall data leads to an exploration of alternative approaches of disaggregation to generate synthetic hourly rainfall.For this reason, the empirical model we proposed for the region shall be used with caution.Continuous monitoring and recording of finer resolution (e.g., hourly/sub-hourly) rainfall is important in helping to improve the IDF relationship by improving storm water management.Also, in order to better understand the regional rainfall trend, non-stationary modelling will be executed and, accordingly, the regional IDF curves will be developed.

Conclusions
This work demonstrated monsoonal rainfall variability over a geographic region in the Central-Western hills of Nepal and attempted to establish monsoon season IDF relationships where fine-resolution rainfall data was scarce, common to many developing countries.The reported methodology and available tools are useful and applicable for ungauged locations within the study region.
The available daily data were evaluated by grouping them in subsets according to the hydrological year.This approach was helpful to better understand the rainfall variability and better fitting of PDFs, as there was a strong seasonal effect on the rainfall in the region where more than 80% of rain occurred in the four monsoonal months.The Panchase hill geographically divides the region into an eastern (i.e., Kaski District) and western (i.e., Parbat and Syangja Districts) sub-regions.Analysis of the above 30 years of homogeneous daily rainfall of eight stations in the region showed that the monsoonal rainfall amount was higher in the eastern part than was recorded in the western part.Also, the numbers of monsoonal dry days were found to be increased over the period of 30 years in the eastern part, with the constant amount of monsoonal rain indicating that the rainfall intensity is increased; whereas in the western part no clear trend was noticed in the annual monsoonal rainfall amount except in Syangja, where dry days are decreased with an increased number of storms.However, because of geographical complexity, in order to better understand the complex rainfall process in the region more detailed analysis is needed.
The proposed methodology to develop the IDF relationship and empirical model in a data-scare situation was found to be statistically significant, for which the available daily data were disaggregated to hourly synthetic series.The recorded daily and disaggregated hourly rainfall data were useful for computing the reference IDF for all the stations.The reference IDF was effective in evaluating the mathematically computed station-specific IDF relationships.The mathematical model parameters (λ, ψ, θ and η) were estimated according to the established methods demonstrated by other researchers (e.g., [22]).Moreover, the adjustment of the empirical constants performed better statistical significance of the IDF and demonstrated the usefulness of the model for the ungauged locations within the region.To regionalize the IDF relationships, the adopted L-moment based RFA method was implemented, dividing the region into two sub-regions first; later it was noticed that the region contains some degree of heterogeneity, leading to the production of five sets of empirical models.The models were used to estimate the regional IDF for the given duration and return periods.Although the regional models underestimated the rainfall intensity for the shorter duration and return periods while fitting them into the reference IDF, they were significant in most cases.This procedure can be extended for larger areas to establish the IDF relationship in the data-scarce situation of Nepal, leading to better management of storm water, roadside drainage design, management and the mitigation of various hazards-risks induced by mass movement.
The station-specific (Gumbel) empirical IDF model demonstrated that the short duration (i.e., 5, 10, 30, 60 min) rainfall intensity in the western part was higher than in the eastern part (highest in Karki-Neta) whereas this was not the case in the regional IDF model.The regional model revealed that the eastern sub-region 2 (Lumle area) received the intense rain.The reason behind the different intensities could be the chosen distribution model.In addition, the regional empirical IDF model generalized the distribution model to be fitted for the region and used the regional mean (extreme) precipitation leading to variation in the intensities that come from the station-specific IDF model.The interpretation of the IDF clearly indicated that the Panchase hill range distinctly divides the region into two meteorological regions.The variability of the rainfall in terms of rain volume and intensity may have different effects on mass movement, soil erosion and floods in the region, which have to be investigated for better and effective hazard-risk management.

Figure 1 .
Figure 1.Study Region indicating the location (name and station number) of weather stations of the DHM and those installed from the Ecosystem Protecting Infrastructures and Communities (EPIC) project (Table 1 for summary).

Figure 1 .
Figure 1.Study Region indicating the location (name and station number) of weather stations of the DHM and those installed from the Ecosystem Protecting Infrastructures and Communities (EPIC) project (Table 1 for summary).

Figure 2 .
Figure 2. Graphical representation of the Bartlett-Lewis rectangular pulse (BLRP) model according to[61].Filled and open circles denote, respectively, the storm origins and cell arrivals[41].

Figure 3 .
Figure 3. (a) Seasonal rainfall pattern and (b) monsoonal season extreme events and number of dry days (1982-2015).The numbers in the x-axis represent the stations.

Figure 3 .
Figure 3. (a) Seasonal rainfall pattern and (b) monsoonal season extreme events and number of dry days (1982-2015).The numbers in the x-axis represent the stations.

Figure 5 .
Figure 5. L-moment ratio diagram for six distributions (note the black plus and square block representing the western region whereas the blue plus and red block is for eastern region).

Figure 5 .
Figure 5. L-moment ratio diagram for six distributions (note the black plus and square block representing the western region whereas the blue plus and red block is for eastern region).

Table 1 .
Summery statistics of historical daily rainfall (1982-2015), missing values, and homogeneity test of all 11 stations.

Table 2 .
Estimated modified BLRP (MBLRP) parameters to be used in HyetosMinute derived from the test rainfall data series.

Table 3 .
Adopted MBLRP model parameters implemented in HyetosMinute derived from the historical monsoonal (daily) rainfall series of the homogeneous data set of eight stations.

Table 3 .
Adopted MBLRP model parameters implemented in HyetosMinute derived from the historical monsoonal (daily) rainfall series of the homogeneous data set of eight stations.

Table 4 .
Comparison of statistical characteristics of recorded and disaggregated rainfall series; respectively, the mean, variance and skewness (E n , Var n , Skew n and E n , Var n , Skew n ).

Table 5 .
Probability distribution function (PDF), estimated parameters and test statistics for the monsoonal daily time series and annual extremes.

Table 7 .
Standard error of prediction (SEP) and chi-square test-statistics (alpha = 0.05, test statistics = 12.59 for degree of freedom = 7) before and after the calibration of the mathematical model.

Table 9 .
List of evaluated distribution indicating the goodness-of-fit (GoF) test statistics (Z DIST < 1.64 to accept the distribution).

Table 9 .
List of evaluated distribution indicating the goodness-of-fit (GoF) test statistics (Z DIST < 1.64 to accept the distribution).

Table 10 .
Chosen distribution for their fitted and estimated empirical parameters.

Table 10 .
Chosen distribution for their fitted and estimated empirical parameters.

Table 11 .
Standard error of prediction (SEP) and chi-square test-statistics (alpha = 0.05, test statistics = 12.59 for degree of freedom = 7) before and after the calibration of the regional model.

Table 12 .
L-moment based regional empirical IDF model for the Panchase region.
t d = duration in hour, T r = return period in year and µ hr = hourly (extreme) mean precipitation.

Table 13 .
Example of computed rainfall intensity (mm/h) for some durations and return periods of regions/sub-regions.