Intensity-duration-frequency curves at the global scale

Intensity-duration-frequency (IDF) curves usefully quantify extreme precipitation over various durations and return periods for engineering design. Unfortunately, sparse, infrequent, or short observations hinder the creation of robust IDF curves in many locations. This paper presents the first global, multi-temporal (1–360 h) dataset of generalized extreme value (GEV) parameters at 31 km resolution dubbed PXR-2 (Parametrized eXtreme Rain). Using these data we generalize site-specific studies to show that that GEV parameters typically scale robustly with event duration (r2 > 0.88). Thus, we propose a universal IDF formula that allows estimates of rainfall intensity for a continuous range of durations (PXR-4). This parameter scaling property opens the door to estimating sub-daily IDF from daily records. We evaluate this characteristic for selected global cities and a high-density rain gauge network in the United Kingdom. We find that intensities estimated with PXR-4 are within ±20% of PXR-2 for durations ranging between 2 and 360 h. PXR is immediately usable by earth scientists studying global precipitation extremes and a promising proof-of-concept for engineers designing infrastructure in data-scarce regions.


Introduction
Historical precipitation records are widely employed by civil engineers to compute intensity-durationfrequency (IDF) curves, which are essential for the design of infrastructure like highways (e.g. Brown et al 2013, NYS DoT 2018, urban drainage networks (e.g. Battaglia et al 2003, Brown et al 2013 and dams (e.g. NYS DoEC 1989). Indeed, IDF curves are used to create synthetic rainfalls that permit the sizing of structures for a given return period, often required by local regulations.
However, not all countries have historical rain gauge records that are long or dense enough to compute reliable IDF curves (e.g. Lumbroso et al 2011). The lack of observational data for IDF analysis is particularly true in continents such as Africa (van de Giesen et al 2014) and Asia, where most of the world's urbanization is expected to take place over coming decades (UN DESA 2018). As a result, much new infrastructure is being built in regions where the historical record of precipitation is scarce or uncertain, hindering adequate sizing of water-related works.
The first limitation of the observational data records is the scarcity of spatial coverage. The classical solution to this problem is to interpolate rainfall between weather stations. However, interpolation performs poorly in locations where pluviometers are sparse (e.g. Xu et al 2015, Kumari et al 2017). A more sophisticated approach consists of analyzing regional precipitation patterns to estimate local characteristics such as IDF curves at the location of interest (e.g. Roux and Desbordes 1996, Fowler and Kilsby 2003, Domínguez et al 2018. Most recently, IDF curves have been derived over the continental US (Ombadi et al 2018) using the PERSIANN-CDR satellite-based precipitation dataset (Ashouri et al 2015). However, a global, consistent IDF dataset is still lacking.
The increasing resolution and reliability of global or near-global gridded precipitation datasets represents a key opportunity to develop alternative approaches to tackle engineering challenges such as the Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. correct sizing of flood infrastructure. Global gridded precipitation estimates are typically obtained by meteorological reanalysis (e.g. Uppala et al 2005, Gelaro et al 2017 where weather observations are assimilated by numerical weather prediction models to estimate precipitation. Advanced data fusion schemes have been developed by merging gauge-, satellite-, and reanalysis-based data to generate enhanced global precipitation estimates (e.g. Beck et al 2018, Sun et al 2018. Although global weather data products are widely employed by the Earth science community, their use in engineering is still limited to applications such as wind power generation (e.g. Staffell and Pfenninger 2016, Olauson 2018) or drought monitoring (e.g. Hao et al 2014). This paper is the first attempt to study global IDF relationships using gridded precipitation datasets, in an effort to help address the issue of precipitation data scarcity.
A second issue that is common with precipitation data is that of temporal resolution. In many cases, subdaily IDF records are required for engineering uses because small catchments that are sensitive to brief rainfall events often require appropriate storm-water drainage structures. There have been recent efforts to compile high-resolution rainfall datasets (Blenkinsop et al 2018), however the vast majority of historical precipitation data are still collected at a daily resolution in most parts of the world. Such low temporal resolution presents a major challenge for engineers designing urban water infrastructure, where catchments are commonly tens of hectares with lag times of less than an hour (Berne et al 2004).
This temporal limitation might be overcome by using a duration scaling property of IDF curves to estimate sub-daily IDF from daily precipitation. Extreme precipitation intensities for a given event duration d typically follow a generalized extreme value (GEV) distribution, and it has been shown that the location and scale parameters of the GEV scale with d (e.g. Menabde et al 1999, Veneziano and Furcolo 2002, Bougadis and Adamows 2006, Overeem et al 2008. For instance, it has been demonstrated that the IDF characteristic of a given site could be described by a Gumbel distribution where the parameters follow a power law for durations between 30 min and 24 h (Menabde et al 1999). However, those studies analyzed only a few sites: 2 in South Africa and Australia (Menabde et al 1999), 5 in Canada (Bougadis and Adamows 2006), and 12 in the Netherlands (Overeem et al 2008). Therefore it has not been possible to determine whether this scaling property holds at a global scale. This would be of practical interest in many parts of the world where daily rainfall data are more widely available than sub-daily records.
This paper first uses the ERA5 reanalysis to generate global IDF relationships modeled with a GEV distribution, then investigates the extent to which these relationships scale with the event duration d at a global level. The analysis led to the creation of two datasets. The Parametrized eXtreme Rain-2 (PXR-2) dataset compiles the GEV parameters for 19 event durations spanning 1-360 h, whereas the Parametrized eXtreme Rain-4 (PXR-4) represents the global distribution of the four parameters of a generalized IDF formula.
2. Methodology 2.1. Input data Precipitation data were obtained from the ERA5 deterministic reanalysis (Hersbach and Dick 2016, Copernicus Climate Change Service 2018) at a spatial resolution of 0.25°(∼31 km) and temporal resolution of 1 h. We chose the ERA5 dataset for its high spatial and temporal resolution, its performance (Beck et al 2018), and its permissive usage license. We employ all the complete calendar years available at the time of writing (i.e. 1979-2018).
We use hourly rain gauge records from the MIDAS database of the UK Meteorological Office (Met Office 2012) as a reference dataset for comparison with the reanalysis data. The original dataset contains 650 stations with variable record lengths. We use only the stations where geographical coordinates are provided. Following Blenkinsop et al (2017), we keep only the observations that do not exceed by more than 20% the 1h and 24 h precipitation historical maxima for the UK, defined respectively as 92 mm and 279 mm by Met Office (2018). After this quality control, we omit the years with <90% of observations remaining, and retain only those stations with 90% of years remaining over the entire time period. This results in a subset of 35 stations for analysis (see figure S6 available online at stacks.iop.org/ERL/14/084045/mmedia for their location).

Estimation of distribution parameters
The GEV distribution is widely used to represent annual maxima series (AMS) (e.g. Fowler and Kilsby 2003, Overeem et al 2008, Papalexiou and Koutsoyiannis 2013 with the cumulative distribution function (CDF): where I is the rainfall intensity, μ the location parameter, σ the scale parameter and κ the shape parameter. Parameters μ and σ have the same unit as I, here mm h −1 . This formulation of the GEV implies that the distribution is bounded from below, corresponding to the Fréchet distribution, when κ<0 (Hosking et al 1985, Overeem et al 2008; other authors use a formulation of equation (1) that implies the opposite sign of κ (e.g. Papalexiou and Koutsoyiannis 2013, Ragulina and Reitan 2017). When κ=0, the GEV is the Gumbel distribution.
The estimation of κ is notoriously difficult due to its sensitivity to record length (Overeem et al 2008, Papalexiou and Koutsoyiannis 2013, Ragulina and Reitan 2017. Indeed, it has been demonstrated that for sample sizes less than 50, using a two parameter Gumbel distribution results in a smaller error than the three parameter GEV (Lu and Stedinger 1992). However, using κ=0 could considerably underestimate the precipitation intensity, especially for long return periods (Koutsoyiannis 2004a(Koutsoyiannis , 2004b. This represents an important safety concern when designing engineering structures like dams. On the other hand, studies have found that κ tends to −0.114 irrespective of geographical location (Papalexiou and Koutsoyiannis 2013) and event duration d (Overeem et al 2008). Considering the difficulty of robustly estimating κ, and that the value tends to −0.114 independently of duration or geographical location, we decide to globally set the GEV shape parameter to κ=−0.114.
The parameters of location (μ) and scale (σ) are estimated using the probability-weighted moments (PWMs) method (Hosking et al 1985). The confidence intervals of the parameter estimates are obtained via bootstrapping with 1000 samples. We assess the goodness of fit of the distribution using the Filliben test (Wilks 2011). This test is based on a Q-Q plot between the AMS and the quantile estimate (mm h −1 , see equation (2)). The probability of exceedance of the AMS is estimated with the Cunnane plotting position (Wilks 2011), and the quantile estimate is calculated at the same return period. The Filliben test statistic consists of the Pearson's correlation coefficient r on that Q-Q plot. The critical value for the test is estimated using the formula given by Heo et al (2008). If rr crit , the null hypothesis that the AMS follows the GEV distribution cannot be rejected.
After fitting the GEV, intensity for a given probability is estimated by the inverse of the CDF, also called the quantile function: where y is expressed relative to the return period T in years with T=1/(1−F): With k = -0.114, the variance of the estimated rainfall intensity i can be obtained with the formula proposed by Lu and Stedinger (1992): where n is the length of the AMS in years. The confidence interval of the precipitation intensity I can then be estimated as: where * -t n 2 is the quantile of the Student's t distribution with n−2 degrees of freedom. For a sample size of 40 years and a 95% confidence interval, this value is 2.024.

Scaling of distribution parameters
To assess the scaling of the distribution parameters (μ, σ) relative to the event duration, we find the annual maxima for 19 event durations d by using a rolling mean. Window sizes are chosen to reflect a relatively regular spacing on a logarithmic scale and to present an equal number of durations for sub-and super-daily events. The selected durations are 1, 2, 3, 4, 6, 8, 10, 12, 18, 24, 48, 72, 96, 120, 144, 192, 240, 288 and 360 h. Subsequently, the GEV parameters are estimated for each duration and ERA5 cell. Considering the spatial resolution of 0.25°, the 19 durations and the 1000 samples bootstrap, in total the GEV is fitted 1.97×10 10 times to the AMS. Global maps of the estimated GEV parameters for each duration are compiled in the PXR-2 dataset (Courty et al 2019), alongside their uncertainties.
Following Menabde et al (1999), we assume that μ and σ scale with d according to a power law, but where they assert a single scaling gradient for both parameters we allow each to scale independently. This independent scaling of the two parameters appears typical for ERA5 data (figures S4 and S5). The scaling is therefore expressed as where d is the event duration in hours and a and α are the scaling parameters for μ, and b and β those for σ. Substituting μ d and σ d into equation (2) with their scaled form from equation (6) we obtain a general IDF formula equation (7) that takes the parameters a, α and b, β specific to a given geographical location.

Global GEV parameters
The PXR-2 dataset comprises worldwide GEV parameters estimated from the ERA5 data for all 19 durations (1-360 h), along with their uncertainties. This set of parameters is made freely available to accompany this paper (Courty et al 2019).
The Filliben goodness of fit test indicates that at the 5% significance level the null hypothesis that the annual maxima follow a GEV distribution with κ=−0.114 can be rejected in 5.7% of the fitted cells (geographical location and duration). Considering the number of tests involved (19.7×10 6 ), ≈5% of rejection is expected at the 5% significance level (Bland and Altman 1995). Here, the lower goodness of fit occurs in spatially coherent areas, located mostly in the drier regions of the Atlantic and Pacific oceans, the Sahara, the south of the Arabian peninsula, tropical Africa and the upper Amazon (see figure S1). Similar goodness of fit estimates are obtained with the Lilliefors test (see figure S2). The values of the test statistic for goodness of fit are included in the PXR dataset (Courty et al 2019).
The GEV parameter maps shown in figure 1 clearly display regional rainfall patterns, such as tropical rainfall and monsoon (e.g. south Asia, Kripalani et al 2007), orographic rainfall over mountainous regions (e.g. central Andes, Viale et al 2011), and desert areas (e.g. Antartica, Vaughan et al 1999).

Scaling of the GEV parameters
The fit of the power law relationship between μ or σ and d is quantified by calculating Pearson's r 2 for data presented on a log-log scale. In 99% of the ERA5 cells r 2 exceeds 0.91 for μ and 0.88 for σ. Thus, the GEV parameters (μ, σ) both scale linearly on a log-log scale and this property appears to be robust and consistent at the global scale. The spatial distribution of r 2 values is also more consistent for μ than σ, with greater spatial variability for the latter (see figure S3).
The global distribution of the scaling gradients (α, β) is shown in figure 2. The spatial distribution of α seems to follow patterns with steeper gradients in deserts (e.g. Baja California, Patagonia, Sahara, Namib, Arabian peninsula, Taklamakan, Gobi) and smaller gradients in mountain ranges (e.g. Andes, Sierra Nevada, East African Rift, Scandes, Alps, Ural, Alborz, Caucasus, Himalaya, Kamchatka). β appears to follow similar patterns, but associated with a greater spatial variability than α. Figure 3 illustrates how this scaling applies for selected global cities. The goodness of fit varies depending on the city, and the fitted regression lines tend to overestimate both GEV parameters at the shortest durations. Additionally, the scale parameter σ displays a weaker linear scaling and higher uncertainty than μ, a property that is in accordance with the r 2 values. Figure 4 shows the differences between the parameter estimate (PXR-2) and the values derived from the scaling relationship (PXR-4) for both the whole world using ERA5, and at MIDAS rain gauges in the UK (see section 2.1). In the case of the value of the whole world estimated from ERA5 (i.e. the PXR dataset), the application of the scaling relationship induces an overestimation of both μ and σ when < d 3 h, of μ when > d 50 h, and of σ when > d 100 h. When limited to the selected MIDAS gauges the differences for PXR follow the same shape. The differences due to scaling at the MIDAS stations are smaller when using the gauge data, especially for σ, and most notably for < d 3 h. Those differences in parameter estimates between PXR-2 and PXR-4 translate to the intensity estimates, with little variation due to the return period (see figure S7).

Global GEV parameters
Our analysis (section 3.1, figure S1) concurs with previous work (e.g. Papalexiou and Koutsoyiannis 2013) suggesting that the AMS of precipitation intensities are usefully described by a GEV distribution. However, using the newly-compiled PXR-2 dataset (Courty et al 2019) we show this applies both on a global scale and with gridded precipitation data. PXR provides a useful simplified description of global extreme precipitation. By describing the entire intensity-frequency distribution for a given d with only two parameters (i.e. not mean, median, mode, range etc), more meaningful inter-comparison between areas is facilitated, as has been done in other fields in analogous situations (e.g. Hillier et al 2013). The utility of PXR is enhanced by the relative ease with which the GEV parameters and their spatial distribution (e.g. figure 1) can be interpreted. Higher μ indicates greater typical precipitation intensities (i.e. the entire distribution becomes more intense), whilst higher σ values indicate more extreme events in the 'tail' of the distribution. Additionally, we showed in section 3.1 that the parameter maps constituting PXR-2 represent qualitatively the expected geographical patterns of extreme precipitation, such as monsoon . Another possible application is as diagnostics for climate and weather models to assess their capacity to reflect the same scalings as those observed in nature.
Our results suggest that μ is broadly more robust than σ. Indeed, the estimates of σ reveal more variability than those of σ in both space (figure 1), duration (figure 3), and the scaling property (figures 2, 4, and S3). This higher variability of σ might be explained by the fact that the scale parameter is related to the intensity of less probable events (i.e. the tail of the probability density function). We acknowledge also that this work employs a relatively short AMS of 40 years that may well miss more extreme events. Indeed, using a longer series of annual maxima is key to improving estimates of GEV parameters (Papalexiou and Koutsoyiannis 2013), although at the risk of overlooking the non-stationary nature of precipitation distributions (Westra et al 2014).
Longer AMS could also allow a better estimate of κ. The fixed value of κ=−0.114 used in this study, based on point rainfall, might not be optimum for areal rainfall and could impact the intensity estimates for long return periods. We acknowledge that the precipitation events considered in this study might not be entirely independent, either in duration or in space, and that the actual confidence interval could be wider than reported here (e.g. Ouarda et al 2019, Wilks 2016).

Scaling in duration of the GEV parameters
Our study confirms that the GEV parameters μ and σ scale robustly with the duration d (e.g. Menabde et al 1999, Overeem et al 2008, and demonstrate that this relationship applies globally. However, in contrast to previous work there is strong evidence that μ and σ scale with different gradients (see section 3.2). As a caveat, we note that the relationship between the parameters and d may be multi-scale (as denoted by breaks in slope of the log-log plots), and that more sophisticated scaling laws may be specified (Clauset et al 2009). This multi-scaling property has been observed in a similar analysis of radar rainfall (Overeem et al 2009), but is not obvious in previous work using gauges (Menabde et al 1999, Overeem et al 2008. We suspect that the scaling differences between ERA5 and the rain gauges (see figure 4) are due to two main factors. First, the weather model used to generate ERA5 might underestimate the actual precipitation intensities of events of shorter durations, which are likely to be convective in nature and of limited spatial scale (Prein et al 2015). This is especially important because the precipitation in ERA5 is modeled, not directly assimilated. Second, the actual scaling property of a gridded product might be inherently different to the scaling of precipitation measured at a point, due to the averaging of rainfall occurring on gridded data. Therefore those differences in scaling might not be due to an inadequacy of the scaling hypothesis, but to an under-reporting of precipitation for events of < d 3 h in the ERA5 dataset.
Therefore, in addition to providing sub-daily IDF information in parts of the world where no such data are readily available, PXR-4 also gives an insight about the feasibility of using daily rainfall records from pluviometers to estimate sub-daily IDF. Indeed, daily records are more common than data from automatic sub-daily gauges, and the lack of the latter is a challenge for engineers (e.g. Lumbroso et al 2011). Additionally, the parsimonious representation in PXR-4 permits the generation of IDF curves for a continuous range of durations rather than discrete d in the case of PXR-2.
The PXR datasets represent areal precipitation which, as would be expected, results in lower intensities than gauges. For some applications an areal representation is actually preferred, as many hydrological processes of interest take place at the catchment scale. Decades of research have generated insights into the relationship between point and areal rainfall, and the estimation of the areal reduction factor (ARF) (e.g. To illustrate the practical impacts of using PXR-4 or PXR-2, we compare the sizing of a culvert in an hypothetical 80 ha catchment in Jakarta with a time of concentration = T 2 h c . Compared to the use of PXR-2, the use of PXR-4 results in an increase in the catchment outflow by 9.5% for the 10 year rainfall and 14.4% for the 100 year rainfall, which in turn induces a modification in the culvert sizing from 1 to 1.2 m for the 10 year rainfall and 1.2-1.5 m for the 100 year event. In this case, PXR-4 yields a more precautionary and potentially costly design than PXR-2. However, as discussed previously, more research is needed to identify whether those differences are the result of an underestimation of short rainfall intensities from ERA5, or an overestimation due to the scaling law used in PXR-4. For further information on this worked example, see the sizing calculations in section S1.1 and the IDF curves for Jakarta from PXR-2 in figure S8.

Concluding remarks
Our results demonstrate the promising applicability of (1) reanalysis data to estimate IDF relationships, and (2) daily rainfall records to estimate sub-daily IDF curves. Our findings may be of notable interest for engineers working in data scarce regions and earth scientists studying extreme precipitation variations. For durations between 2 h and 360 h, the precipitation intensities estimated from PXR-4 are within ±20% of those estimated from PXR-2 (see figure S7). . Median percentage error (MdPE) between the GEV parameters μ, σ estimated using the scaling relationship ad α , bd β (as in PXR-4) and the actual μ d , σ d estimated by fitting the AMS (as in PXR-2). The greater the deviation from zero, the less accurate are α and β at estimating μ and σ. Values above zero indicate an overestimation of the GEV parameter.
We acknowledge that our research is only a first step in the analysis of global IDF and scaling relationships. More work is needed on the subject, including: • Studying the range of validity of both PXR-2 and PXR-4 relative to catchment sizes and when compared to other options in data-scarce areas.
• Evaluating the nature of the scaling property, including subhourly precipitation.
• Determining the physical causes of the scaling.
• Comparing the scaling property in both empirical observations and in climate model output.
• Investigating the seasonal variations in the occurrence and scaling of extreme precipitation.
• Evaluating the impact of climate change on the scaling property (e.g. Ouarda et al 2019).
• Assessing the use of other probability laws to describe precipitation extremes (e.g. Marani and Ignaccolo 2015, Ouarda et al 2019).
• Evaluating the implications of using a fixed κ to fit the GEV, and assessing alternative values or methodologies.
• Studying the impact of precipitation types (snow, hail etc) and processes (orographic, convective etc) on the scaling relationship.
• Examining the underlying factors of the spatial variability observed in both the GEV parameters and their scaling gradients.
In the mean time, PXR-2 and PXR-4 are shared with the research and engineering communities to foster further work on global precipitation extremes.