# Identification of Design Rainfall Changes Using Regional Frequency Analysis: A Case Study in Ontario, Canada

## Abstract

Historical rainfall records at 32 gauges located in Southern Ontario, Canada were analyzed to assess if temporal changes in design rainfall intensities were evident in this geographical area. Rainfall records were split into two time periods and design rainfalls of 2 y, 5 y and 10 y return periods were estimated using partial duration series instead of the conventionally used annual maximum series. However, due to limited available record lengths, uncertainties in design rainfall estimates were substantial and determinations of whether storm intensities were increasing were challenging. To reduce uncertainties, the gauges were grouped into regions and a regional L-moment algorithm was used for each region.

We found that the regional algorithm produced more accurate rainfall estimates (i.e. lower relative root mean square error) and less uncertainty (i.e. narrower confidence intervals) than at-site models. This approach shows, for example, that for 10 y return storms, the regional model reduced relative root mean square error by 26% for the first period (1960–1983) and 35% for the second period (1984–2007). The findings indicate that this procedure indeed identified increasing trends in rainfall intensities for the design rainfalls for 30% of the sites (11 of 32). Further, the regional models showed 66% of the sites have at least one location of increasing design rainfalls.

## 1 Introduction

Climate change research has been catapulted to the fore in recent years. As a result of the increased energy present in the hydrologic cycle, more intense precipitation events are expected, and are evident (Alexander et al. 2006; Adamowski et al. 2010; Burn and Taleghani 2012). If changes are occurring, they are expected to be evident in the mean values and variations of precipitation intensities through changes in the extreme rainfall recurrence frequency. Consequently, real challenges exist in the design and management of urban stormwater systems. For example, a stormwater sewer system may be designed to convey a 5 y event specified in the applicable stormwater design manual but the actual design rainfall intensities are changing. Further, estimating design rainfall intensities for the future requires a comprehensive understanding of changes that have been observed in relation to extreme rainfall events over the past several decades.

Types of temporal changes can be characterized as step changes or gradual changes. In frequency analysis, especially when the objective is to assess the design rainfall intensity, an independent and identically distributed record is normally assumed, and gradual changes are removed (de-trended) if identified. Therefore, in frequency analysis, using periods of record to estimate extreme values and to identify step changes rather than gradual changes over time is more beneficial to the assessment of design rainfall intensities.

Statistics used for identifying changes over time include expected values and variances of design rainfalls, for both a return period and duration. Cumulative density functions of extreme rainfall records may also be analyzed to discover step changes. Descriptive indices of extremes have also been used to describe changes (WMO 2009), such as R95p, the number of days with rainfall above the 95th percentile of daily accumulations.

Various assessments of changes in precipitation rates have been identified across Canada. Groisman (1999) observed a 50% increase in mean summer precipitation over the past century based on daily precipitation records. Zhang et al. (2001) characterized heavy precipitation events for the period 1900–1998 and found increasing trends of heavy rainfall during spring periods in Eastern Canada. Vincent and Mekis (2006) separately analyzed daily precipitation records for periods 1950–2003 and 1900–2003, and found that in the latter half of the twentieth century the number of days with precipitation increased, while daily intensity and maximum number of consecutive dry days have both decreased. Burn and Taleghani (2012) identified more increasing trends than decreasing trends based on records of 51 stations across Canada and found that the design rainfalls of various return periods have increased in the most recent 20 y period, compared to the entire record.

Rainfall trends in Southern Ontario have also been studied. Stone et al. (2000) grouped the Great Lakes and St. Lawrence area as one homogeneous region when analyzing daily precipitation events across Canada. The paper identified seasonally increasing trends in total precipitation in Southern Ontario and concluded that more extreme precipitation events in autumn and winter were related to the negative Pacific–North American teleconnection (PNA). Adamowski and Bougadis (2003) discovered both increasing and decreasing trends of extreme rainfall events in various rain gauge stations in Ontario.

Changes in design rainfalls are discovered by comparing estimates of design rainfall values from the records of two different time periods. Burn and Taleghani (2012) compared the estimate from the most recent 20 years to the estimate from the entire timeframe. Since the record lengths of rain gauge stations in Ontario are usually short (mostly starting from the 1960s), splitting the record into two segments results in even shorter records (approximately 20 y in each segment). A short record will entail more uncertainties in the estimates of design rainfall values. To solve this problem, Burn and Taleghani (2012) used resampling techniques to estimate the uncertainty (lack of precision) of quantile estimates.

As seen from previous research, frequency analysis is a commonly used technique in identifying changes. The analysis may be performed for individual sites (referred to as at-site analysis) or for a region consisting of multiple sites possessing similar characteristics (referred to as regional frequency analysis, RFA). RFA uses rainfall records from adjacent rainfall stations that are in a statistically homogeneous region. It assumes that sites in a region have identical frequency distributions except for site-specific scale factors. Hosking and Wallis (1997) introduced several comprehensive measures to delineate homogeneous regions and to select frequency distributions, the regional L-moment algorithm, and the methods used to assess the accuracy of estimated quantiles. This analysis has been demonstrated to be an effective way to reduce the uncertainty of quantile estimates for design rainfall intensity.

RFA has been used in various applications. For example, Madsen et al. (1998, 2002) analyzed regional variability in extreme rainfall statistics by using linear regression between site statistics (e.g. index–flood) and site characteristics (e.g. mean annual precipitation) and developed a regional estimation model for precipitation in Denmark. Trefry et al. (2005) updated intensity–duration–frequency (IDF) curve estimates for Michigan using regional frequency analyses based on both annual maximum series (AMS; the data series consisting of the maximum value in each year) and partial duration series (PDS; the data series consisting of all values exceeding a threshold) data. Results show that RFA can provide reliable rainfall IDF estimates. Bonnin et al. (2006) describe the RFA approach in a precipitation–frequency atlas of the United States. Ngongondo et al. (2011) compared the accuracy of design rainfall estimates between at-site and regional analysis for 23 rainfall stations in Malawi and concluded that the regionally based estimates have smaller uncertainties. Sveinsson et al. (2002) analyzed regional extreme precipitation frequencies in northeastern Colorado, concentrating on an extraordinary storm that occurred in 1997. Darwish et al. (2018) used regional frequency analysis in the United Kingdom, and Khan et al. (2017) used L-moments and partial L-moments in assessing the regional frequency analysis of extreme events in Pakistan. Dad and Benabdesselam (2018) used RFA on maximum daily precipitation for northeastern Algeria and showed that an at-site model underestimated the quantities having a high return period for >80% sites in their study. Hao et al. (2019) assessed RFA of precipitation extremes using L-moments in China, and Forestieri et al. (2018) performed RFA for Sicily using AMS from ~130 rain gauges.

More research papers using RFA techniques are based on AMS. Laurenson (1987) states the advantage of using PDS over AMS for at-site analyses: the rainfall model using PDS data evaluates rainfall intensity of the average recurrence interval between storm events instead of between two hydrologic years, in which a given rainfall intensity is exceeded regardless of the number of exceedances. The difference between rainfall intensity models based on PDS and AMS is discussed by Vrban et al. (2018), who used an RFA approach to improve the accuracy of quantile estimates for heavy events (i.e. small return periods) based on PDS data.

## 2 Methodology

### 2.1 Differences Between the Use of PDS and AMS Data

lmomRFA is an RFA package implemented in the programming language R (Hosking 2012; R Core Team 2013). The simulation algorithms in the package need to be updated when using PDS data in RFA, since the record length is not the same as the number of years of the rainfall record, as currently required by lmomRFA. PDS is a data series extracted from the historical record by selecting rainfall events, with a corresponding time of occurrence, that exceed a given threshold *xT*; the record length is likely different from the count of values in PDS. Reference is made to Vrban et al. (2018) for details of PDS construction. An average annual arrival rate *λ* is calculated by dividing the number of years of the record into the number of values in the record. The original regional L-moment algorithm uses AMS, and storms of the same return period have the same non-exceedance probability *F* for all gauges. This consistency in *F* is not preserved when using PDS data due to the differences in *λ* between rain gauges. In this case, the projection from the regional frequency distribution to the at-site storm estimate depends on both the scale factor *ℓ*_{1} (mean rainfall intensity at gauge) and *λ* (average annual arrival rate); that is, *Q**̂** _{l}*(

*T*) =

*ℓ*

_{1}

^{(}

^{i}^{)}

*q*

*̂*(1 − 1/(

*λ*

*T*)) for the design rainfall intensity with a return period

*T*y at station

*i*(where

*q*

*̂*is the quantile function of the fitted regional frequency distribution). The non-exceedance probability is calculated as

*F*= 1 − 1/(

*λ*

*T*).

The PDS data were extracted with thresholds *x** _{τ}* specific to the individual rain gauges. The sensitivity of distribution parameters (using a generalized Pareto distribution) and design rainfall estimates with respect to the thresholds was analyzed following Coles (2001). The threshold was selected from a range when distribution parameters were relatively stable and rainfall intensity estimates were close to the interpolated values of ranked rainfall records using Gringorten’s plotting position formula (Gringorten, 1963). Vrban et al. (2018) discusses the choice of selection of intensity thresholds.

Gross errors introduced by instrument malfunction or mistakes in transcription need to be eliminated before the data can be used for frequency analyses. Hosking and Wallis (1997, p. 48) suggest a discordancy measure to identify gauges that are discordant with other gauges in a group, which we used.

### 2.2 Identifying Homogeneous Regions

The RFA algorithms assume that rainfall records within a homogeneous region have similar frequency distributions, excepting for scale factors. The homogeneity of a region is in relation to the statistical similarity of frequency distributions between records, which indicates gauges within a homogeneous region are not necessarily close in terms of statistical proximity. If the rainfall pattern of an area is highly affected by geographical factors such as oceans, lakes or mountains, then geographical characteristics should also be considered in grouping the gauges. Other site characteristics (e.g. mean annual precipitation, MAP, and time of year at which extreme events mostly occur) are also widely used (e.g. Adamowski et al. 1996; Trefry et al. 2005; Ngongondo et al. 2011).

In this study, the gauges being considered were located in Southern Ontario which has a continental climate markedly modified by the Great Lakes (Hare and Thomas 1979). Therefore, the gauge distances to three lakes, Lake Huron, Lake Erie and Lake Ontario, were used as site characteristics in the grouping of gauges.

Clustering methods are practical when dealing with large datasets. Gauges are partitioned or aggregated into groups based on similarities of their at-site characteristics or statistics. Trefry et al. (2005) used a *k*-means clustering method to group gauges, and Ngongondo et al. (2011) used *k*-means and Ward’s hierarchical method to cluster rain gauges in southern Malawi. This study uses Ward’s hierarchical clustering method to initially group rain gauges, and then uses a *k*-means clustering algorithm to adjust clusters.

Regional homogeneity is assessed using a measurement of the degree of heterogeneity in a group of gauges. This measurement compares two components: between gauge dispersion of the sample L-moment ratios and between gauge dispersion if the group of gauges is heterogeneous. Between gauge dispersion is the average difference between gauge statistics and average group statistics weighted by the length of record (in years) at each gauge. The dispersion of a homogeneous region is calculated using Monte Carlo simulation. In each realization, records are simulated with the same length as their real world counterparts, using a parent frequency distribution (usually a four parameter Kappa distribution).

Between gauge dispersion of sample coefficients of L-variation (*L*-*CVs*) is calculated by:

(1) |

where:

V |
= | dispersion of L-CVs, |

N |
= | number of rain gauges, |

t^{i} |
= | L-CVs at gauge i, and |

t^{R} |
= | group average of L-CVs, weighted by record length n._{i} |

The heterogeneity measure (*H*) is calculated by:

(2) |

where:

μ and _{V}σ_{V} |
= | mean and standard deviation of the dispersions calculated by Monte Carlo simulations. |

The classification of *H* is *acceptable homogeneity* when *H* <1, *possible heterogeneity* when 1 ≤ H < 2, and *definite heterogeneity* when H ≥2.

### 2.3 Selection of Frequency Distribution

The goodness-of-fit measure introduced in Hosking and Wallis (1997) is designed to select between candidate frequency distributions. It assumes that the variations of L-moment ratios in a homogeneous region are due to sampling variability; therefore, the candidate distributions are evaluated by how well the fitted L-skewness and L-kurtosis match the regional average L-skewness and L-kurtosis of the observed data. For three-parameter candidate distributions, the L-skewness is fitted to regional averages; thus only the difference of L-kurtosis between fitted distribution *τ*_{4}* ^{DIST}* and regional average

*t*

_{4}

*is evaluated:*

^{R}(3) |

where:

Z^{DIST} |
= | the goodness of fit measure, |

t_{4}^{R} |
= | regional average L-kurtosis, |

DIST |
= | candidate distribution, and |

σ_{4} |
= | standard deviation of t_{4}.^{R} |

As for the heterogeneity measure, the Monte Carlo simulation is used to quantify the bias *B*_{4} and variability *σ*_{4} of *t*_{4}* ^{R}*:

(4) | |

(5) |

where:

m |
= | index in the N replications of the Monte Carlo simulation, and_{sim} |

t_{4}^{[m]} |
= | regional average L-kurtosis in m-th replication. |

With the bias of *t*_{4}* ^{R}*, the goodness-of-fit measure is modified by:

(6) |

The criterion |*Z ^{DIST}*| ≤ 1.64 is used (Hosking and Wallis1997) to judge if

*Z*is sufficiently close to zero with a 90% confidence level (i.e. the L-kurtosis of the fitted distribution is close to the regional average L-kurtosis of the observed data). A lower value of

^{DIST}*Z*suggests a better distribution fit.

^{DIST }### 2.4 Regional L-moment Algorithm

The regional L-moment algorithm is based on the index–flood method, which averages the statistics of data at the gauge to form regional estimates. Instead of using conventional moments in the index–flood method, the regional L-moment algorithm uses L-moment ratios of the data. The regional L-moment algorithm assumes that no serial correlation for data is observed at the same gauge, and there is no dependence between observations at different gauges. The index flood at gauge *i*, or the scale factor, is estimated by the sample mean of record and denoted as *l*_{1}^{(}^{i}^{)}. Data at each gauge are rescaled using this index flood; therefore, the regional average mean rainfall intensity is unity.

The L-moment ratios at gauge *i* are denoted as *t*^{(}^{i}^{)}, *t*_{3}^{(}^{i}^{)} and *t*_{4}^{(}^{i}^{)} and the regional average L-moment ratios are denoted by *t*^{(}^{R}^{)}, *t*_{3}^{(}^{R}^{)} and *t*_{4}^{(}^{R}^{)}; they are calculated by:

(7) | |

(8) |

The parameters of the selected frequency distribution are estimated based on the regional average L-moment ratios, and the quantile function with non-exceedance probability *F* is calculated by:

(9) |

where:

= | quantile estimate at gauge i at non-exceedance probability F, and | |

= | quantile estimate for the region at non-exceedance probability F. |

**2.5 Assessment of Estimated Quantile Accuracy**

The accuracy of the estimated quantile in regional L-moment algorithm is estimated using Monte Carlo simulation. The simulated samples should have the same statistical characteristics as those of observed data to maintain the heterogeneity, dependence and other statistical characteristics of the observed data as well as the possibility when the frequency distribution is wrongly specified.

Hosking and Wallis (1997) provide a detailed description of the assessment procedure, so we only discuss the modifications applied for using PDS.

A correlation matrix *R* is used to indicate the between site dependencies. It is originally calculated by:

(10) |

where:

= | , | |

n_{i} |
= | record length of gauge i; and |

Q_{ik} |
= | data value for gauge i at time point k. |

The time point *k* extends over all time points for which both gauges *i* and *j* have values. Equation 10 works well for AMS data since there is only one value per year. However, PDS data may have several values recorded in the same year. Thus *Q _{ik}* is redefined as the maximum value for gauge

*i*at time point

*k*.

We note that in Equation 10 data values recorded in a year at one gauge without counterparts at another gauge are excluded from the calculation of between-site dependence. However, in the simulation algorithm, these data values are generated together with the data values included in the dependence calculation using the correlation matrix *R*. When using PDS data, the second largest values in a year are excluded when calculating between-site-dependence. These data values are also calculated using the correlation matrix *R*.

The correlation matrix *R* needs to be positive definite to calculate correlated variables using Cholesky decomposition. Nonpositive–definite correlation matrices are updated by changing negative eigenvalues to small positive values (0.000 000 01) and normalized (Brissette et al. 2007).

Over *M* repeated simulations, the relative root mean square error (rRMSE) is approximated by:

(11) |

where:

R_{i}(F) |
= | rRMSE at gauge i, |

Q̂_{i}^{[}^{m}^{]}(F) |
= | quantile estimate at the non-exceedance probability F of replication m, and |

Q(_{i}F) |
= | estimated quantile based on observed data, at gauge i. |

The at-site rRMSE can be averaged over all gauges within a region to obtain a regional averaged rRMSE by:

(12) |

where:

N |
= | the number of gauges in the region. |

The rRMSE calculated in Equation 12 is useful to quantify the variances in estimates of the regional model, and to compare to its counterpart in the estimates of the at-site model. However, confidence intervals of design rainfall estimates are needed to identify step changes in the design rainfall intensity over time. As discussed in Wang et al. (2015), a step change is statistically significant if the confidence intervals of design rainfall estimates for different time periods are not overlapping and the significance level is greater than or equal to the confidence level.

The calculation of the confidence interval requires assumptions such as independence between rainfalls from different gauges, statistically homogeneous regions, and properly selected regional statistical distributions. However, it is hard to satisfy all these assumptions. The rainfall data usually present some degree of violation of the assumptions. Instead, the empirical quantiles of the distribution of estimates are useful assessments of errors. In Monte Carlo simulation, the ratio of the estimated value to the true value [*Q**̂** _{l}*(

*F*)/

*Q*(

_{i}*F*)] at site

*i*are accumulated over each realization, and the upper and lower 90% percentiles are found and denoted as

*U*

_{.05}(

*F*) and

*L*

_{.05}(

*F*). The true value

*Q*(

_{i}*F*) is expressed as:

(13) |

where:

= | estimated quantity, | |

U_{.05}(F) |
= | 95th percentile of , and |

L_{.05}(F) |
= | 5th percentile of . |

Hosking and Wallis (1997) referred to these bounds as *90% error bounds* for *Q**̂** _{l}*(

*F*) and indicated that they can be confidence intervals only if the distribution of

*Q*

*̂*

*(*

_{l}*F*)/

*Q*(

_{i}*F*) is independent of the at-site means and the regional average L-moment ratios. In practice, “the independence does not hold, and confidence statements are at best approximate” (Hosking and Wallis 1997). The error bounds can be accurate estimates of the magnitude of errors when the number of repetitions

*M*is large (e.g.

*M*= 1 000).

## 3 Application of a Regional Frequency Model

### 3.1 Data Collection and Screening

Records for 270 rainfall gauges located in Ontario, Canada were obtained from Environment Canada. The gauge data have daily maximum 1 h rainfall amount records for various time spans, from as early as 1937 to as late as 2009. Records of 44 gauges were combined with records from other gauges at the same locations or in very close proximity that had consecutive rainfall records. The research period selected was 1960–2007 and the dataset was split at the end of 1983 (as the mid-point) to identify changes in extreme rainfall intensities by comparing estimates from rainfall records pre- and post-1983. The rainfall records for the period April 1–October 31 (seven months) were extracted, and data for any seven month period that had >20% values missing were excluded from the analyses. We analysed data from 32 gauges in Southern Ontario (Figure 1 and Table 1). Due to the high density of gauges in Southern Ontario, all the gauges had records for >10 y for both pre- and post-1983 periods. The data had previously been subjected to quality control by Environment Canada, so it was assumed that any gross errors had been eliminated.

Table 1 Information for 32 gauges in Southern Ontario.

Climate ID |
Gauge Location | Longitude | Latitude | Record Length (y) | Region | ||

1960–1983 | 1984–2007 | ||||||

1 | 6100971 | BROCKVILLE PCC | −75.67 | 44.60 | 16 | 17 | 1 |

2 | 6104027 | KEMPTVILLE CS | −75.63 | 45.00 | 13 | 17 | 1 |

3 | 6104175 | KINGSTON PUMPING STATION | −76.48 | 44.24 | 22 | 23 | 1 |

4 | 6105978 | OTTAWA CDA RCS | −75.72 | 45.38 | 20 | 20 | 1 |

5 | 6106000 | OTTAWA MACDONALD-CARTIER INT’L A | −75.67 | 45.32 | 14 | 20 | 1 |

6 | 6127514 | SARNIA AIRPORT | −82.30 | 42.99 | 16 | 21 | 3 |

7 | 6131415 | CHATHAM WPCP | −82.22 | 42.39 | 17 | 17 | 3 |

8 | 6131983 | DELHI CS | −80.55 | 42.87 | 19 | 20 | 4 |

9 | 6133362 | HARROW CDA AUTO | −82.90 | 42.03 | 16 | 12 | 3 |

10 | 6136606 | PORT COLBORNE | −79.25 | 42.88 | 16 | 14 | 4 |

11 | 6137287 | ST CATHARINES A | −79.17 | 43.20 | 11 | 18 | 4 |

12 | 6137362 | ST THOMAS WPCP | −81.21 | 42.77 | 20 | 21 | 3 |

13 | 6139148 | VINELAND STATION RCS | −79.40 | 43.18 | 13 | 10 | 4 |

14 | 6139525 | WINDSOR A | −82.96 | 42.28 | 23 | 20 | 3 |

15 | 6140954 | BRANTFORD MOE | −80.23 | 43.13 | 22 | 11 | 4 |

16 | 6142400 | FERGUS SHAND DAM | −80.33 | 43.73 | 20 | 21 | 3 |

17 | 6143090 | GUELPH TURFGRASS CS | −80.22 | 43.55 | 20 | 14 | 5 |

18 | 6144478 | LONDON CS | −81.15 | 43.03 | 23 | 23 | 3 |

19 | 6146714 | PRESTON WPCP | −80.35 | 43.38 | 10 | 12 | 3 |

20 | 6148105 | STRATFORD MOE | −81.00 | 43.37 | 16 | 17 | 3 |

21 | 6149387 | WATERLOO WELLINGTON A | −80.38 | 43.45 | 13 | 18 | 3 |

22 | 6150689 | BELLEVILLE | −77.39 | 44.15 | 17 | 24 | 1 |

23 | 6150830 | BOWMANVILLE MOSTERT | −78.67 | 43.92 | 14 | 17 | 2 |

24 | 6151042 | BURKETON MCLAUGHLIN | −78.80 | 44.03 | 13 | 16 | 2 |

25 | 6153194 | HAMILTON A | −79.93 | 43.17 | 13 | 19 | 4 |

26 | 6153301 | HAMILTON RBG CS | −79.91 | 43.29 | 17 | 22 | 4 |

27 | 6155878 | OSHAWA WPCP | −78.83 | 43.87 | 12 | 16 | 2 |

28 | 6158355 | TORONTO CITY | −79.40 | 43.67 | 24 | 24 | 2 |

29 | 6158665 | TORONTO ISLAND A | −79.40 | 43.63 | 13 | 11 | 2 |

30 | 6158733 | TORONTO LESTER B. PEARSON INT’L A | −79.63 | 43.68 | 24 | 22 | 2 |

31 | 6158875 | TRENTON A | −77.53 | 44.12 | 17 | 15 | 1 |

32 | 6166418 | PETERBOROUGH A | −78.37 | 44.23 | 11 | 18 | 2 |

Figure 1 Thirty-two gauges in 5 clusters located in Southern Ontario.

### 3.2 Identification of Regions

The regions were identified by clustering site characteristics, which include distances to Lake Huron, Lake Erie, and Lake Ontario. In Southern Ontario, sources of wet and warm air could be as close as the Great Lakes, or as far as the Atlantic Ocean. The distances of rain gauges to the Atlantic Ocean are not sufficiently distinct to be the clustering Z variable; therefore, the distances to three of the Great Lakes were used.

Using Ward’s hierarchical clustering method and the k-means clustering algorithm, five regions were initially formed. However, some regions are either too small in size, or statistically heterogeneous. The heterogeneity measure tends to indicate false homogeneity for small regions. Thus, regions are adjusted to obtain larger regions and to achieve statistical homogeneity. In particular, the rainfall record (1984–2007) at the Guelph gauge has a negative L-kurtosis ratio and a small L-skewness ratio, which is statistically distinct from any other gauges; therefore, this gauge was split off as a single-site region. Four regions were identified and the number of gauges in each region and corresponding measures of heterogeneity are listed in Table 2. Regions 2 and 4 show some degree of heterogeneity; therefore the regenerated samples should preserve equivalent heterogeneity.

Table 2 Number of gauges in each region and related heterogeneity measures.

Region ID | # of Gauges | Heterogeneity Measure (1960–1983) |
Heterogeneity Measure (1983–2007) |

1 | 7 | 0.58 | −0.54 |

2 | 7 | 1.66 | 1.25 |

3 | 10 | −0.07 | 0.28 |

4 | 7 | 1.23 | 1.68 |

### 3.3 Choice of Frequency Distribution

Candidate frequency distributions are selected on the L-moment ratio diagram. Figures 2 and 3 show separate diagrams for each region for the pre- and post-1983 periods. In each diagram, a frequency distribution is selected as a candidate if it is proximate to the center of the at-site L-moment ratios. The five frequency distributions under consideration include: generalized logistic (GLO), generalized extreme value (GEV), generalized Pareto (GPA), generalized normal (GNO) and Pearson type III (PE3). Figures 2 and 3 show that all five distributions are acceptable for selection. The goodness-of-fit measure *Z ^{DIST}* (Equation 6) is computed for each of the five candidate distributions for each region and compiled in Table 3. The bias of regionally averaged L-kurtosis is considered since all regions have <20 gauges. For both pre- and post-1983 time periods, PE3 is chosen for use in regions 1 and 2, and GPA is chosen for use in regions 3 and 4, since the corresponding measures for goodness-of-fit are smaller in comparison with other types of distributions.

Table 3 Measures of goodness-of-fit for candidate distributions in each region.

Region ID | 1960–1983 | 1984–2007 | ||||||||

GLO | GEV | GNO | PE3 | GPA | GLO | GEV | GNO | PE3 | GPA | |

1 | 4.44 | 3.56 | 2.58 | 0.69 | 0.82 | 3.97 | 3.29 | 2.13 | 0.18 | 0.92 |

2 | 3.28 | 2.5 | 1.7 | 0.17 | 0.44 | 3.23 | 2.43 | 1.58 | 0.09 | −0.12 |

3 | 3.07 | 2.33 | 0.92 | −1.34 | −0.31 | 4.16 | 3.2 | 1.77 | −0.83 | 0.09 |

4 | 3.7 | 2.54 | 1.78 | 0.53 | −0.01 | 2.41 | 1.67 | 0.8 | −0.67 | −0.37 |

### 3.4 Assessment of Quantile Estimate Accuracy

In the Monte Carlo simulation, a series of steps were applied to reproduce regions with the same statistical characteristics as the original region. The heterogeneity was preserved by randomly assigning a range of L-CV values to gauges in the region. The L-skewness ratios are all replaced with the regional averaged L-skewness ratios. Interested readers are referred to Hosking and Wallis (1997, p. 95). For example, period 1960–1983 of region 1 has seven gauges with L-CVs varying in the range 0.134–0.205, and the heterogeneity measure is *H* = 0.58. The simulated gauges have L-CVs varying in the range 0.137–0.202 with a heterogeneity measure *H* = 0.55.

In Figures 4 and 5, between-gauge correlations are relatively low in all regions, with average values ranging from −0.02 to 0.21. Thus the samples generated in the Monte Carlo simulation were not correlated. For each period, 1000 replications of each region were simulated, and rRMSE and error bounds (90%) were calculated.

### 3.5 At-site PDS Model

The design rainfall intensity estimates based on at-site PDS data were also calculated for comparison with estimates from the regional model. Statistical distributions used for at-site PDS models were consistent with the distribution used for the region in which the gauge was located. The rRMSE and 90% error bounds were estimated in the same way as the regional L-moment algorithm (Equation 11 and Equation13) except for the use of at-site estimates instead of regional estimates.

## 4 Results of Accuracy Improvement

Comparisons between the regional model and the at-site model show improvement in the accuracy of design rainfall intensity estimate. Table 4 shows the comparison of regional averaged rRMSE between two models. The regional averaged rRMSE is calculated using Equation 12, and the ratios in Table 4 are [*R ^{R}(F)* for at-site model /

*R*for regional model]. The reduction in rRMSE is shown in all regions, for rainfall intensity estimates of all three evaluated return periods (2 y, 5 y and 10 y).

^{R}(F)Table 4 Ratio of regional averaged root mean square error of the regional models versus the at-site models.

Region ID | 1960–1983 | 1984–2007 | ||||

2 Year | 5 Year | 10 Year | 2 Year | 5 Year | 10 Year | |

1 | 0.82 | 0.75 | 0.69 | 0.78 | 0.66 | 0.61 |

2 | 0.88 | 0.85 | 0.81 | 0.86 | 0.78 | 0.71 |

3 | 0.88 | 0.79 | 0.68 | 0.83 | 0.74 | 0.64 |

4 | 0.87 | 0.88 | 0.81 | 0.91 | 0.90 | 0.81 |

Figure 6 Histograms of ratio of error bound extent between the regional model and the at-site model.

Table 5 Ratio of regional model root mean square error against at-site model root mean square error.

Gauge Location | 1960–1983 | 1984–2007 | ||||

2 Years | 5 Years | 10 Years | 2 Years | 5 Years | 10 Years | |

Brockville PCC | 0.59 | 0.54 | 0.47 | 0.68 | 0.63 | 0.62 |

Kemptville CS | 0.80 | 0.64 | 0.50 | 0.82 | 0.65 | 0.56 |

Kingston Pumping Station | 0.93 | 0.93 | 0.89 | 0.76 | 0.65 | 0.63 |

Ottawa CDA RCS | 0.77 | 0.75 | 0.75 | 0.96 | 0.81 | 0.74 |

Ottawa Macdonald-Cartier Int’l A | 0.89 | 0.89 | 0.87 | 1.01 | 0.84 | 0.80 |

Belleville | 0.93 | 0.90 | 0.89 | 0.76 | 0.79 | 0.79 |

Trenton A | 1.01 | 0.85 | 0.80 | 0.58 | 0.46 | 0.37 |

Bowmanville Mostert | 1.09 | 1.05 | 0.94 | 1.13 | 1.14 | 1.12 |

Burketon McLaughlin | 1.69 | 1.76 | 1.70 | 0.84 | 0.73 | 0.68 |

Oshawa WPCP | 0.78 | 0.67 | 0.60 | 1.00 | 1.02 | 1.02 |

Toronto City | 0.85 | 0.88 | 0.87 | 0.69 | 0.63 | 0.58 |

Toronto Island A | 0.72 | 0.66 | 0.63 | 0.58 | 0.43 | 0.36 |

Toronto Lester B. Pearson Int’l A | 0.79 | 0.84 | 0.84 | 1.25 | 1.32 | 1.26 |

Peterborough A | 0.75 | 0.76 | 0.74 | 0.95 | 0.93 | 0.91 |

Sarnia Airport | 0.92 | 0.94 | 0.87 | 0.66 | 0.57 | 0.47 |

Chatham WPCP | 0.66 | 0.58 | 0.52 | 1.07 | 1.23 | 1.23 |

Harrow CDS Auto | 0.87 | 0.71 | 0.55 | 0.75 | 0.69 | 0.62 |

St Thomas WPCP | 0.83 | 0.82 | 0.71 | 0.83 | 0.80 | 0.75 |

Windsor A | 1.14 | 1.19 | 1.13 | 1.13 | 1.02 | 0.90 |

Fergus Shand Dam | 0.68 | 0.54 | 0.42 | 0.81 | 0.66 | 0.51 |

London CS | 0.87 | 0.81 | 0.74 | 1.00 | 0.98 | 0.89 |

Preston WPCP | 1.09 | 0.96 | 0.85 | 0.87 | 0.83 | 0.75 |

Stratford Moe | 0.87 | 0.83 | 0.79 | 0.73 | 0.59 | 0.49 |

Waterloo Wellington A | 0.98 | 0.81 | 0.69 | 0.70 | 0.54 | 0.45 |

Delhi CS | 1.45 | 1.93 | 2.08 | 0.77 | 0.86 | 0.82 |

Port Colborne | 0.98 | 1.20 | 1.30 | 0.69 | 0.69 | 0.62 |

St Catharines A | 0.61 | 0.54 | 0.44 | 1.09 | 1.38 | 1.36 |

Vineland Station RCS | 0.74 | 0.69 | 0.61 | 0.82 | 0.83 | 0.80 |

Brantford Moe | 1.10 | 0.94 | 0.81 | 1.42 | 1.13 | 0.91 |

Hamilton A | 0.95 | 1.25 | 1.34 | 0.77 | 0.64 | 0.53 |

Hamilton RBG CS | 0.80 | 0.73 | 0.68 | 1.21 | 1.48 | 1.53 |

Guelph Turfgrass CS | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 |

Amongst the 32 gauges analyzed, 22 gauges had smaller rRMSEs in the regional model than in the at-site model for 1960–1983, and 22 gauges (not the same set of gauges as referred to above) for 1984–2007. Figure 6 shows the histograms of ratios of extent of rRMSE between two models for both time periods. Comparisons of the error bounds between the two models show almost identical results.

The gauges that did not show improvement in accuracy in the regional model are mostly in regions with possible heterogeneities. In Table 5, for regions 2 and 4, where both regions have *H* >1, a substantial number of gauges have larger rRMSE for the regional model than the at-site models. These gauges are indicated by symbols consisting of crosses and circles on the L-moment ratio diagrams in Figures 2 and 3 for the two time periods. In Figure 3 (L = Logistic, N = Normal and U = Uniform), the L-moment ratios of three gauges (23, 27 and 30) are all plotted close to the regional distribution curve (PE3); in contrast, gauge 24, which was plotted away from the PE3 curve, was estimated more accurately in the regional model. This indicates that closeness to the regional distribution curve on the L-moment ratio diagram is not an indicator of more accurate estimates.

## 5 Changes Identified in Design Rainfall Intensities

As discussed in Vrban et al. (2018), a step change is identified if the confidence intervals of design rainfall estimates for the two time periods do not overlap, with a significance level no less in magnitude than the confidence level. In addition, the 90% error bounds are good approximations of the confidence interval when the number of simulations is large (*M* = 1000). Therefore, we used the error bounds to identify changes in the design rainfall intensity.

For the 1 h design rainfall, based on the at-site PDS model, changes were identified at five gauges in Southern Ontario. In detail, the gauges located at Ottawa Macdonald-Cartier International Airport and Toronto Pearson International Airport show statistically significant decreases; while the Delhi, Belleville, and Bowmanville gauges increased significantly. The annual change rate is calculated by:

(14) |

where:

i and _{1}i_{2} |
= | design rainfall intensities estimated from periods 1 and 2, separately; and |

t and _{1}t_{2} |
= | the mid-point years of periods 1 and 2. |

The complete results of changes identified based on the at-site PDS model are listed in Table 6, including changes which are not significant. Among the 32 gauges, both the at-site models and the regional models identify increasing trends for the gauge Kemptville CS (i.e. positive values 0.79%/y, 0.92%/y, 0.96%/y, 0.62%/y, 0.80%/y and 0.93%/y) along with 10 other gauges that represent 30% of the sites. At gauge Brockville PCC, the regional model corrects the change direction (i.e. from 0.23%/y decreasing for the at-site model to 0.12%/y increasing for the regional model). The regional models also correct the direction changes for six other gauges. The regional models estimate that 21 sites (66%) will have at least one increasing design rainfall.

Table 6 Changes in design rainfall intensity in Southern Ontario (%/y, using at-site PDS).

Gauge Location | At-Site Model | Regional Model | ||||

2 Years | 5 Years | 10 Years | 2 Years | 5 Years | 10 Years | |

Brockville PCC | −0.23 | −0.48 | −0.63 | 0.12 | 0.32 | 0.44 |

Kemptville CS | 0.79 | 0.92 | 0.96 | 0.62 | 0.80 | 0.93 |

Kingston Pumping Station | 0.01 | 0.27 | 0.43 | −0.04 | 0.11 | 0.20 |

Ottawa CDA RCS | −0.73 | −0.69 | −0.64 | −0.50 | −0.37 | −0.28 |

Ottawa Macdonald−Cartier Int’l A | −1.29* | −1.09 | −0.94 | −1.18* | −0.95 | −0.81 |

Belleville | 1.15* | 1.29 | 1.38 | 1.02* | 1.06* | 1.10 |

Trenton A | 1.01 | 1.81 | 2.31 | 0.71 | 0.89 | 1.00 |

Bowmanville Mostert | 1.81* | 1.35 | 1.08 | 1.88* | 1.41 | 1.17 |

Burketon McLaughlin | 0.51 | 1.29 | 1.86 | 0.83 | 0.51 | 0.34 |

Oshawa WPCP | 1.22 | 0.37 | −0.15 | 0.94 | 0.61 | 0.42 |

Toronto City | −0.56 | −0.43 | −0.34 | −0.53 | −0.61 | −0.65 |

Toronto Island A | −2.11 | −1.69 | −1.39 | −1.79* | −1.82 | −1.84 |

Toronto Lester B. Pearson Int’l A | −0.66 | −1.05* | −1.26* | −0.54 | −0.62 | −0.66 |

Peterborough A | −0.64 | −1.02 | −1.20 | −0.43 | −0.73 | −0.89 |

Sarnia Airport | 0.64 | 1.02 | 1.36 | 0.43 | 0.35 | 0.27 |

Chatham WPCP | −0.50 | −1.10 | −1.52 | −0.17 | −0.25 | −0.31 |

Harrow CDS Auto | 0.66 | 0.49 | 0.29 | 0.42 | 0.36 | 0.30 |

St Thomas WPCP | 0.20 | 0.06 | −0.08 | 0.18 | 0.16 | 0.13 |

Windsor A | −0.61 | −0.52 | −0.43 | −0.54 | −0.52 | −0.52 |

Fergus Shand Dam | 0.03 | −0.15 | −0.30 | 0.05 | −0.02 | −0.08 |

London CS | 0.02 | −0.16 | −0.31 | 0.12 | 0.04 | −0.02 |

Preston WPCP | 2.62 | 2.81 | 2.78 | 2.14 | 1.93 | 1.77 |

Stratford Moe | 0.35 | 0.60 | 0.89 | 0.42 | 0.28 | 0.17 |

Waterloo Wellington A | 1.63 | 2.00 | 2.30 | 1.47* | 1.23 | 1.07 |

Delhi CS | 0.97* | 1.75* | 2.38* | 0.46 | 0.52 | 0.58 |

Port Colborne | −0.06 | 0.44 | 0.96 | 0.01 | −0.02 | −0.01 |

St Catharines A | −0.46 | −1.19 | −1.72 | −0.46 | −0.47 | −0.45 |

Vineland Station RCS | −0.62 | −0.67 | −0.72 | −0.62 | −0.54 | −0.48 |

Brantford Moe | −0.07 | −0.27 | −0.32 | 0.08 | 0.04 | 0.05 |

Hamilton A | −0.43 | 0.40 | 1.41 | 0.05 | 0.06 | 0.09 |

Hamilton RBG CS | 0.00 | −0.37 | −0.63 | 0.05 | 0.01 | 0.02 |

Guelph Turfgrass CS | 0.65 | −0.45 | −1.26 | 0.65 | −0.45 | −1.26 |

Five gauges showed statistically significant changes in the regional model, as shown in Table 7. The differences compared to the at-site model results are: no changes identified at Delhi or Toronto Pearson International Airport gauges; changes identified at the Toronto Island Airport gauge (2 y return storm, decreased), Waterloo–Wellington Airport (2 y return storm, increased), and Belleville (5 y return storm, increased).

Table 7 L-CVs and rainfall intensity estimates of gauges identified as showing different changes from the regional and the at-site model.

Gauge | T | Model | 1960–1983 | 1984–2007 | ||||||

L-CV | Q(F) | 0.05 | 0.95 | L-CV | Q(F) | 0.05 | 0.95 | |||

Delhi CS | 5yr | regional | 0.11 | 29.71 | 26.66 | 33.36 | 0.11 | 33.89 | 29.69 | 38.59 |

at-site* | 0.12 | 27.17 | 25.10 | 28.99 | 0.20 | 36.70 | 31.64 | 42.35 | ||

10yr | regional | 0.11 | 33.20 | 29.08 | 38.05 | 0.11 | 39.10 | 33.43 | 46.23 | |

at-site* | 0.12 | 28.95 | 26.47 | 31.43 | 0.20 | 42.73 | 35.19 | 51.01 | ||

Toronto Lester B. Pearson Int’l Airport | 5yr | regional | 0.21 | 34.96 | 30.68 | 39.94 | 0.18 | 30.01 | 26.93 | 33.68 |

at-site* | 0.20 | 36.59 | 31.25 | 42.14 | 0.12 | 27.77 | 25.35 | 30.70 | ||

10yr | regional | 0.21 | 40.71 | 34.78 | 47.72 | 0.18 | 34.53 | 30.14 | 39.72 | |

at-site* | 0.20 | 43.21 | 35.87 | 51.04 | 0.12 | 30.65 | 27.44 | 34.63 | ||

Toronto Island Airport | 2yr | regional* | 0.19 | 28.09 | 24.84 | 31.66 | 0.17 | 22.05 | 19.63 | 24.76 |

at-site | 0.20 | 28.26 | 24.02 | 33.01 | 0.20 | 21.11 | 17.08 | 25.96 | ||

Waterloo–Wellington Airport | 2yr | regional* | 0.24 | 24.77 | 21.64 | 27.47 | 0.23 | 30.41 | 27.93 | 33.61 |

at-site | 0.17 | 24.25 | 20.97 | 27.77 | 0.22 | 30.36 | 26.28 | 34.93 | ||

Belleville | 5yr | regional* | 0.19 | 25.51 | 22.62 | 28.63 | 0.20 | 31.60 | 28.67 | 35.02 |

at-site | 0.16 | 24.63 | 21.43 | 28.13 | 0.19 | 31.78 | 27.99 | 35.93 |

## 6 Discussion

The rainfall intensities at gauges in the same region were assumed to be statistically identically distributed apart from scale factors (defined as the sample means at each gauge). This assumption implies that storms observed at one gauge are likely to be observed at other gauges in the same region, both historically and in the future. Therefore, the between-site variations in the occurrences of heavy storms, which substantially determine the at-site model estimates, are averaged over all gauges in the regional model (except for the variation in scale factors). For regions that are not acceptably homogeneous, the heterogeneities between gauges are preserved in the Monte Carlo simulations by means of assigning L-CV values for replications of rainfall records. The reassigned L-CV values vary over a range that is not the same as the original L-CV values. In consequence, new variations between gauges may not be the same amounts as in the original rainfall records.

For the gauges at Delhi CS and Toronto Pearson International Airport, the rearrangement of sample variation changed the design rainfall intensity estimate in a direction opposite to the direction of changes identified in the at-site model (as listed in Table 6); therefore the differences between estimates from the two time periods were substantially reduced, and no statistically significant differences could be identified in the regional model. The gauge at Delhi CS had an original L-CV = 0.20 over the period 1984–2007 but was simulated with a value of 0.11 in the regional model. The reassignment of L-CV values resulted in reduced rainfall intensity estimates in the regional model. The Toronto Pearson International Airport gauge had an original L-CV = 0.12 and a regional simulated L-CV = 0.18 over the second period; therefore the estimates were substantially increased for 5 y and 10 y return events. The 2 y event estimates at Toronto Island Airport were changed as the L-CV values also changed. The error bounds at the Waterloo–Wellington Airport and Belleville gauges only slightly overlapped in the at-site model, and the changes in L-CV values led to completely unconnected error bounds. The correlation coefficients between *Q**̂*(*F*)* _{regional}*/

*Q*

*̂*(

*F*)

*and*

_{at-site}*L-CV*/

_{simulated}*L-CV*were calculated for all gauges, where

_{original}*F*is the non-exceedance frequency of 2 y, 5 y, and 10 y return events. The correlation coefficients were all statistically significant (0.95, using the Student t-test) for both time periods, which supports the relationship between intensity estimates and simulated L-CV values.

The regional L-moment algorithm was expected to improve the accuracy of the rainfall intensity estimates, especially for gauges with limited rainfall records. This expectation is partially supported by the results of accuracy improvement. It shows that for region 1, which is statistically homogeneous, the rRMSE values of the 10 y storm estimates were reduced on average by 26% (ranging from 11% to 53%) in the regional model for the first time period and by 35% (ranging from 20% to 63%) for the second time period. Given that the record lengths are ~20 y, and therefore may be inadequate to estimate 10 y return storms, the reduction in rRMSE is very impressive. The performance of the regional model on the 5 y return storm estimates is encouraging: an average 22% (first period) and 31% (second period) reduction in rRMSE. For regions with statistical heterogeneity (e.g. regions 2 and 4) many gauges with record lengths <15 y show reduced rRMSE values in the regional model, while gauges with relatively longer records have increased rRMSE values in the regional model. This is as expected since RFA is more suitable to reduce the uncertainty of estimates for sites with shorter records. For sites with longer records, RFA may not reduce the uncertainty. The improvement in accuracy of rainfall intensity return period estimates for the regional model has considerable merit.

## 7 Conclusion

Traditionally, regional frequency analysis (RFA) is performed using the annual maximum series. This has been demonstrated to be less accurate than partial duration series (PDS) for shorter return periods. We updated the RFA using PDS and implemented the R algorithm based on lmomRFA. The application to rain gauges in Southern Ontario showed that:

- The regional L-moment algorithm can reduce the uncertainty involved in rainfall intensity estimates, especially at gauges with limited records and within statistically homogeneous regions.
- The rainfall intensity estimates for the regional model were sensitive to the simulated L-CV values, which determine the between-site variation in a region. The design rainfall intensity changes identified in the regional model can be different from the at-site PDS model, due to the reassignment of between-site variations (L-CV values) in the regional model.
- The regional L-moment algorithm using PDS data, as introduced in this study, can develop a regional model with respect to PDS extracted with different annual arrival rates (
*λ*). This allows capture of the benefits of using PDS in heavy rainfall estimates, including improved design rainfall estimates, in contrast to an AMS model. - Among the 32 rain gauges in the case study, both the regional models and the at-site models identified increasing trends for 11 sites. The regional models provided improved accuracy over the at-site models (i.e. reduced root mean square error for most sites and return periods). Accordingly, the regional models provided corrected increasing or decreasing trend estimates for seven sites and significant changes for five sites compared to those provided by the at-site models.

## Acknowledgments

The research funding provided by the Canada Research Chair program and the Ontario Research Foundation, and the historical rainfall records provided by Environment Canada, are gratefully acknowledged.

## References

- Adamowski, J., Adamowski, K. and Bougadis, J. 2010. “Influence of trend on short duration design storms.”
*Water Resources Management*24:401–13. - Adamowski, K., Alila, Y. and Pilon, P. J. 1996. “Regional rainfall distribution for Canada.”
*Atmospheric Research*42:1–4, 75–88. - Adamowski, K. and Bougadis, J. 2003. “Detection of trends in annual extreme rainfall.”
*Hydrological Processes*17 (18): 3547–60. - Alexander, L. V., Zhang, X., Peterson, T. C., Caesar, J., Gleason, B., Klein Tank, A. M. G., Haylock, M., Collins, D., Trewin, B., Rahimzadeh, F., Tagipour, A., Rupa Kumar, K., Revadekar, J., Griffiths, G., Vincent, L., Stephenson, D. B., Burn, J., Aguilar, E., Brunet, M., Taylor, M., New, M., Zhai, P., Rusticucci, M. and Vazquez-Aguirre, J. L. 2006. “Global observed changes in daily climate extremes of temperature and precipitation.”
*Journal of Geophysical Research: Atmospheres*111:22. - Bonnin, G., Martin, D., Lin, B., Parzybok, T., Yekta, M. and Riley, D. 2006.
*Precipitation–frequency atlas of the United States*. NOAA Atlas 14. Silver Spring, MD: National Oceanic and Atmospheric Administration. - Brissette, F., Khalili, M. and Leconte, R. 2007. “Efficient stochastic generation of multisite synthetic precipitation data.”
*Journal of Hydrology*345 (3–4): 121–33. - Burn, D. H. and Taleghani, A. 2012. “Estimates of changes in design rainfall values for Canada.”
*Hydrological Processes*27 (11): 1590–9. - Coles, S. 2001.
*An introduction to statistical modeling of extreme values*. London: Springer-Verlag. - Dad, S. and Benabdesselam, T. 2018. “Regional frequency analysis of extreme precipitation in northeastern Algeria.”
*Journal of Water and Land Development*39:27-37, DOI: 10.2478/jwld-2018-0056 - Darwish, M., Fowler, H., Blenkinsop, S. and Tye, M. 2018. “A regional frequency analysis of UK sub-daily extreme precipitation and assessment of their seasonality.”
*International Journal of Climatology*38 (13): 4758–76. https://doi.org/10.1002/joc.5694 - Forestieri, A., Conti, F., Blenkinsop, S., Cannarozzo, M., Fowler, H., Noto, L. 2018. “Regional frequency analysis of extreme rainfall in Sicily (Italy).”
*International Journal of Climatology*38:e698-e716 - Gringorten, I. I. 1963. “A plotting rule for extreme probability paper.”
*Journal of Geophysical Research*68 (3): 813–4. - Groisman, P. Y., Karl, T. R., Easterling, D. R., Knight, R. W., Jamason, P. F., Hennessy, K. J., Suppiah, R., Page, C. M., Wibig, J., Fortuniak, K., Razuvaev, V. N., Douglas, A., Førland, E. and Zhai, P.-M. 1999. “Changes in the probability of heavy precipitation: Important indicators of climatic change.”
*Climatic Change*42 (1): 243–83. - Hao, W., Hao, Z., Yuan, F., Ju, Q. and Hao, J. 2019. “Regional frequency analysis of precipitation extremes and its spatio-temporal patterns in the Hanjiang River basin, China.”
*Atmosphere*10 (3): 130. https://doi.org/10.3390/atmos10030130 - Hare, F. K. and Thomas, M. K. 1979.
*Climate Canada*, 2nd ed. Toronto: John Wiley & Sons. - Hosking, J. R. M. 2012.
*R package, version 2.4: Regional frequency analysis using L-moments*. Vienna: R Foundation for Statistical Computing. - Hosking, J. and Wallis, J. 1997.
*Regional frequency analysis: An approach based on L-Moments*. Cambridge: Cambridge University Press. - Khan, S. A., Hussain, I., Hussain, T., Faisal, M., Muhammad, Y. and Shoukry, A. 2017. “Regional frequency analysis of extreme precipitation events using L-moments and Partial L-moments, 2017.”
*Advances in Meteorology*2017:6954902 - Laurenson, E. M. 1987. “Back to basics on flood frequency analysis.”
*Transactions of the Institution of Engineers*29 (2): 47–53. - Madsen, H., Mikkelsen, P., Rosbjerg, D. and Harremoës, P. 1998. “Estimation of regional intensity–duration–frequency curves for extreme precipitation.”
*Water Science and Technology*37 (11): 29–36. - Madsen, H., Mikkelsen, P. S., Rosbjerg, D., and Harremoës, P. 2002. “Regional estimation of rainfall intensity–duration–frequency curves using generalized least squares regression of partial duration series statistics.”
*Water Resources Research*38 (11). https://doi.org/10.1029/2001WR001125 - Mekis, E. and Hogg, W. D. 1999. “Rehabilitation and analysis of Canadian daily precipitation time series.”
*Atmosphere–Ocean*37 (1): 53–85. - Ngongondo, C., Xu, C.-Y., Tallaksen, L., Alemaw, B. and Chirwa, T. 2011. “Regional frequency analysis of rainfall extremes in southern Malawi using the index rainfall and l-moments approaches.”
*Stochastic Environmental Research and Risk Assessment*25:939–55. - R Core Team. 2013. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing.
- Stone, D. A., Weaver, A. J. and Zwiers, F. W. 2000. “Trends in Canadian precipitation intensity.”
*Atmosphere–Ocean*38 (2): 321–47. - Sveinsson, O., Salas, J. and Boes, D. 2002. “Regional frequency analysis of extreme precipitation in northeastern Colorado and Fort Collins flood of 1997.”
*Journal of Hydrologic Engineering*7 (1): 49–63. - Trefry, C. M., Watkins, D. W. and Johnson, D. 2005. “Regional rainfall frequency analysis for the state of Michigan.”
*Journal of Hydrologic Engineering*10 (6): 437–49. - Vincent, L. A. and Mekis, É. 2006. “Changes in daily and extreme temperature and precipitation indices for Canada over the twentieth century.”
*Atmosphere–Ocean*44 (2): 177–93. - Vrban, S., Wang, Y., McBean, E., Binns, A. and Gharabaghi, B. 2018. “Evaluation of stormwater infrastructure design storms developed using partial duration and annual maxima series models.”
*ASCE Journal of Hydrologic Engineering*. https://doi.org/10.1061/(ASCE)HE.1943-5584.0001712 - Wang, Y., McBean, E. A. and Jarrett, P. 2015. “Identification of changes in heavy rainfall events in Ontario, Canada.”
*Stochastic Environmental Research and Risk Assessment*29 (8): 1949–62. - WMO (World Meteorological Organization). 2009.
*Guidelines on analysis of extremes in a changing climate in support of informed decisions for adaptation*. Geneva: WMO. - Zhang, X., Hogg, W. D. and Mekis, É. 2001. “Spatial and temporal characteristics of heavy precipitation events over Canada.”
*Journal of Climate*14 (9): 1923–36.