Intensity–duration–frequency (IDF) rainfall curves in Senegal

Urbanization resulting from sharply increasing demographic pressure and infrastructure development has made the populations of many tropical areas more vulnerable to extreme rainfall hazards. Characterizing extreme rainfall distribution in a coherent way in space and time is thus becoming an overarching need that requires using appropriate models of intensity–duration–frequency (IDF) curves. Using a 14 series of 5 min rainfall records collected in Senegal, a comparison of two generalized extreme value (GEV) and scaling models is carried out, resulting in the selection of the more parsimonious one (four parameters), as the recommended model for use. A bootstrap approach is proposed to compute the uncertainty associated with the estimation of these four parameters and of the related rainfall return levels for durations ranging from 1 to 24 h. This study confirms previous works showing that simple scaling holds for characterizing the temporal scaling of extreme rainfall in tropical regions such as sub-Saharan Africa. It further provides confidence intervals for the parameter estimates and shows that the uncertainty linked to the estimation of the GEV parameters is 3 to 4 times larger than the uncertainty linked to the inference of the scaling parameter. From this model, maps of IDF parameters over Senegal are produced, providing a spatial vision of their organization over the country, with a north to south gradient for the location and scale parameters of the GEV. An influence of the distance from the ocean was found for the scaling parameter. It is acknowledged in conclusion that climate change renders the inference of IDF curves sensitive to increasing non-stationarity effects, which requires warning end-users that such tools should be used with care and discernment.

1 Introduction

The fast-growing pressure of mankind on planet Earth causes populations to be increasingly exposed to hydrometeorological hazards such as torrential rains and floods (IPCC, 2012; Mechler and Bouwer, 2015). Hydrologists are thus more compelled than ever to deal with the problem of assessing the probability of extreme rainfall events at different timescales and for various return periods, depending on the area of the target catchment and the issue at stake, most notably human life protection and infrastructure dimensioning. A classical way of synthesizing the results of such studies is the production of so-called rainfall intensity–duration–frequency (IDF) curves, which provide estimates of rainfall return levels over a range of durations. In doing so, scientists face two sets of difficulties: one related to data availability and the other to the necessity of a proper methodological framework.

On the data side, the frequency analysis of extremes requires long and continuous records of rainfall at the same location, something fairly common at a daily time step albeit unavailable in some regions. Moreover, a complicating factor is that, in many cases, it is necessary to consider sub-daily time steps. However, long-term records of sub-daily rainfall are much less numerous or much less reliable and accurate than daily series.

The methodological challenge arises from the complex combination of factors that cause rainfall to be strongly variable at all scales (from the microphysics droplet scale to synoptic scale), as a result of the nonlinear interaction of different atmospheric processes (e.g., Schertzer and Lovejoy, 1987). This implies that it is not at all obvious to find a proper theoretical framework to compute IDF curves in a way that ensures coherency between timescales. Early works on IDF proposed empirical methods consisting of first adjusting a frequency distribution model fitted to rainfall series {R(D)} for each duration D of interest and then fitting the IDF formula {i_T(D)} fitted independently to each series of quantiles derived from the first step and corresponding to a given return period T (see, e.g., Miller et al., 1973; NERC, 1975). This has the advantage of being easily implementable and is thus commonly used by hydrological engineers and operational climate and hydrological services. However, because of uncertainties in the computation of the quantiles derived for the different durations, the scaling formulation may be physically inconsistent and may lead to gross errors such as parasitic oscillations or intersections between IDF curves computed for two different durations (see Koutsoyiannis et al., 1998, for more details). As a remedy to such inconsistencies, Koutsoyiannis et al. (1998) were the first to propose a general IDF formulation that remains consistent with both the foundations of the probabilistic theories and the physical constraints of scaling across durations. Another notable advance was provided by Menabde et al. (1999), who demonstrated that the changes in rainfall distribution with duration formulated by Koutsoyiannis et al. (1998) can be expressed as a simple scaling relationship, opening the path for using the fractal framework in order to describe the temporal scaling between IDF curves established over a range of durations in various regions of the world (see, e.g., Yu et al., 2004; Borga et al., 2005; Gerold and Watkins, 2005; Nhat et al., 2008; Bara et al., 2009; Blanchet et al., 2016; Rodríguez-Solà et al., 2016; Yilmaz et al., 2016).

However, having a consistent scaling framework does not eliminate the crucial sampling issues associated with the estimation of the parameters of the IDF model. This involves significant uncertainties in the final determination of rainfall return levels, a question rarely addressed in the literature; on that subject, see the pioneering work of Mélèse et al. (2017), which presents and compares different methods for computing IDF confidence intervals (on a GEV and scaling model) over the Mediterranean region. It is especially important to investigate this issue in tropical regions such as sub-Saharan Africa, where a number of new infrastructure projects are in the pipeline while at the same time a significant increase of flood risks and related human casualties has been reported over the past two decades (Di-Baldassarre et al., 2010; Tschakert et al., 2010; IPCC, 2014).

Several recent studies have dealt with the question of IDF calculation for different West African countries. Some focused on analyzing the behavior of the extreme rainfall distribution at a given location (such as Soro et al., 2008, 2010, for Côte d'Ivoire) while others (such as Mohymont and Demarée, 2006, for Congo and Oyegoke and Oyebande, 2008, for Nigeria) looked at the scaling behavior over durations. Van-De-Vyver and Demarée (2010) also analyzed the scaling properties of rainfall over a range of durations for a couple of stations in Congo, finding the value of the main scaling parameter to be larger than the one obtained for Uccle in Belgium, implying that the small durations are heavily driving the behavior of extreme rainfall at larger durations for a tropical climate. De Paola et al. (2014) have also inferred IDF curves from disaggregated daily rainfalls for three African cities (Addis Ababa, Ethiopia; Dar Es Salaam, Tanzania; and Douala, Cameroon).

More recently, Panthou et al. (2014b) and Agbazo et al. (2016) showed that the generalized extreme value (GEV) and simple scaling framework is well-suited to estimate rainfall return levels in a coherent way at various durations for an array of stations covering a mesoscale area of typically a dozen thousands km² in southwest Niger and in northern Benin, respectively. While Agbazo et al. (2016) assumed a Gumbel distribution of the annual maxima, Panthou et al. (2014b) used the approach in its broader formulation, showing that the annual maxima distribution was heavy-tailed (positive value of the shape parameter of the GEV). Indeed, such heavy-tailed behavior in daily rainfall samples is generally found: in the African region (e.g., Panthou et al., 2012; Giugni et al., 2015) but also all around the world (e.g., Koutsoyiannis, 2004b; Papalexiou and Koutsoyiannis, 2013).

It is worth noting that both Panthou et al. (2014b) and Agbazo et al. (2016) made use of the high-quality and fine timescale resolution data collected by the AMMA-CATCH research observatory (Lebel et al., 2009). This data set homogeneously covering a wide range of time steps (from 5 min upward) over more than 20 years is unique in tropical Africa. This means than in every other area, the parameters of the scaling relationship will have to be inferred from a very limited number of sub-daily rainfall series, not all of them being of equal length, thus raising the question of which parameters have the largest influences on the final uncertainties of rainfall return levels. This issue is extremely important when dealing with large regions (such as a whole country) over which the scaling parameters may vary spatially, making it not straightforward to infer rainfall return levels for sub-daily durations when only daily data are available.

Focusing on Senegal, a region of contrasted coastal to inland semi-arid climate, our paper's ambitions are both to address the uncertainty issue not dealt with in the above-mentioned papers and to provide IDF curves for a region located at the western edge of the Sahel, evaluating the spatial variability generated by the transition from the coast to inland. In addition to its methodological bearing, the paper aims at making these IDF curves widely accessible to a large range of end-users in the whole country by mapping the values of the scaling parameters and of the rainfall return levels. Furthermore, selecting an IDF model that is the least sensitive possible to data sampling effects and computing the associated IDF confidence intervals facilitates updating of the IDF curves when new data are available.

2 Data and region

2.1 Senegal climatological context

Senegal is located at the western edge of the African continent between latitudes 12 and 17^∘ N (Fig. 1a). The climate of Senegal is governed by the West African monsoon (Lafore et al., 2011; Janicot et al., 2011; Nicholson, 2013), resulting in a two-season annual cycle: a dry season marked by the predominance of maritime and continental trade winds in winter, and a rainy season marked by the progressive invasion of the West African monsoon (Fig. 2a to c) during the summer. The length of the rainy season varies by latitude and ranges roughly from 5 months (early June to the end of October) in the south to 3 months (mid-July to mid-October) in the northern part of Senegal. Rainfall amounts peak in August and September, coinciding with the period when the Intertropical Convergence Zone (ITCZ) reaches its northernmost position over Senegal.

Figure 1 (a) The location of Senegal is indicated by the black area. (b) A map of Senegal. White circles represent tipping bucket rain gauges (a cross and a label are added for stations having more than 10 years of valid data).

There is a strong north–south gradient of the mean annual rainfall (Fig. 2d) ranging from 300 mm in the north to more than 1000 mm in the south (Diop et al., 2016). This gradient is mainly explained by the number of rainy days (in average between 20 and 80 from north to south) and to a lesser extent by the mean intensity of rainy days (in average between 10 and 15 mm day⁻¹); see Fig. 2e and f.

Figure 2 Rainfall regime statistics obtained from daily rain gauges over the 1950–2015 period: (a) mean seasonal cycle at four stations, (b) start of the rainy season, (c) end of the rainy season, (d) mean annual rainfall, (e) mean intensity of wet days, and (f) mean number of wet days.

The rains are mainly caused by mesoscale convective systems sweeping the country from east to west (Laurent et al., 1998; Mathon et al., 2002; Diongue et al., 2002). Occasionally, cyclonic circulations off the Senegalese coast direct moisture-laden air flow over the western part of the country, dumping heavy rains that often cause floods in coastal cities. Due to its convective nature, rainfall in West Africa is strongly variable in space and time, especially at the event scale for which large differences in rainfall amounts are frequently observed at two nearby points (Sane et al., 2012). The rain durations are also generally short, except in the rare case of stationary convective systems (blocked situation).

Senegal regularly undergoes heavy damaging downpours. A recent example is the rainfall event that occurred in Dakar in the morning of 26 August 2012, causing the largest flood in the last 20 years in the city. An amount of 160 mm was recorded at the Dakar-Yoff station, a large quantity but not a historical record for daily rainfall at this station. Rather, this event was exceptional because of its intensities at short durations (54 mm were recorded in 15 min and 144 mm in 50 min), far exceeding the previous records in Dakar-Yoff. Such rainfall intensities and their associated disasters justify the importance of better documenting extreme rainfall distributions at short timescales.

2.2 Rainfall data

The archives of climate and hydrological services of West African countries sometimes contain large amounts of sub-daily rainfall records. However, most of the time these records are stored in paper strip chart formats, requiring the tedious task of digitization in order to use them in numerical applications.

The present study has been made possible thanks to the important work of analyzing and digitizing rain gauge charts carried out for the main synoptic stations of Senegal. This process was undertaken by the hydro-morphology laboratory of the Geography Department at the University Cheikh Anta Diop of Dakar (UCAD) in collaboration with the National Agency of Civil Aviation and Meteorology (ANACIM) that provided the rainfall paper charts.

Senegalese synoptic stations are equipped with tipping bucket rain gauges; the receiving ring is 400 cm² and a bucket corresponds to 0.5 mm of rain. The roll rotation is daily. The chart analysis has been performed with the software “Pluvio” developed by Vauchel (1992) allowing the computation of 5 min time step digitized rainfall series from the paper diagrams. It is a long and laborious task, which has the advantage of allowing a careful “chart by chart” checking of the quality of the records before digitization. For more information on the digitization process, the reader may refer to the publications of Laaroubi (2007) and Bodian et al. (2016).

A total of 23 tipping bucket rain gauges were analyzed, with data going back to 1955 for the oldest and to 2005 for the most recent. As the assessment of extreme rainfall distributions is known for being highly sensitive to sampling effect and erroneous data (Blanchet et al., 2009; Panthou et al., 2012, 2014b), particular attention was paid to check and select the most appropriate series. The data selection had to reconcile two constraints: (i) keeping the data set as large as possible and (ii) eliminating series that contain too much missing data.

The procedure for classifying 1 station year as valid or not is the following: (i) the annual number of 5 min data and the annual amount of rain are computed, (ii) the mean interannual values of these two statistics are computed on the whole series, (iii) a year is classified as valid if either the number of 5 min rain data or the amount of rainfall is comprised between 1∕2.5 and 2.5 times their mean interannual values, and (iv) other years are classified as missing and removed from the whole series. Since missing years influence the mean interannual values, step (ii), (iii) and (iv) are repeated until all remaining years are classified as valid (note that, in fact, no station year had to be excluded after the initial step). All valid years for all series are plotted in Fig. 3. In order to keep the IDF fitting robust, only series with at least 10 years of valid data have been used. This led us to retain 14 stations with record length varying from 10 (Fatick station) to 44 years (Ziguinchor station) with a median of 28 years. This data set has the advantage of spatially representing the entire country, but as the length of the series varies, the quality of the IDF estimates might differ from one station to another. This effect will be analyzed more precisely in Sect. 5.2.1.

Figure 3 Years of network operation. The valid years are represented by black squares. The total number of valid years is displayed for each station on the right y axis.

3 Theoretical background

3.1 General IDF formulations

3.1.1 Empirical IDF formulations

IDF curves provide estimates of rainfall intensity for a range of durations {D} and for several frequencies of occurrence (usually expressed as a return period T). Each curve corresponds to the evolution of a return level (i_T) as a function of rainfall duration D. Historically, several empirical formulations of IDF curves have been proposed. All can be described by the following general equation (Koutsoyiannis et al., 1998):

$$\begin{array}{}\text{(1)}& i_T\left(D\right)=w\left(T\right)\times [D+\mathit{\theta }(T){]}^{\mathit{\eta }\left(T\right)},\end{array}$$

where w, θ, and η are parameters to be calibrated from rainfall observations.

Figure 4 Kolmogorov–Smirnov GOF p values for both IDF models: (a) in calibration mode; (b) in validation mode; left panel for IDF_SiSca; right panel for IDF_Koutso.

Figure 5 Quantile–quantile plots for both IDF models for the different durations and for all stations (global scores in the legend): (a) in calibration mode; (b) in validation mode; left panel for IDF_SiSca; right panel for IDF_Koutso. The x axis and y axis units are mm h⁻¹.

Figure 6 Comparison between IDF_SiSca and IDF_Koutso in calibration mode. Top panels: location parameter – μ(D) [mm h⁻¹]; scale parameter – σ(D) [mm h⁻¹]; and shape parameter – ξ [–]. Bottom panels: return levels – i_T [mm h⁻¹] – obtained for different return periods T ranging from 2 to 100 years displayed.

Figure 7 IDF_SiSca return level plots obtained at four emblematic stations over Senegal (a–d: Saint-Louis, Dakar, Tambacounda, and Ziguinchor). Shaded area represents the 90 % confidence interval. The markers represent the observed annual maxima intensities i(D).

3.1.2 Koutsoyiannis scaling relationship

Koutsoyiannis et al. (1998) have demonstrated that the empirical formulations (Eq. 1) can be expressed as follows:

$$\begin{array}{}\text{(2)}& i_T\left(D\right)=a\left(T\right)\times b_{\mathrm{Koutso}}\left(D\right),\end{array}$$

where b_Koutso(D) is the scaling function:

$$\begin{array}{}\text{(3)}& b_{\mathrm{Koutso}}\left(D\right)=(D+\mathit{\theta }{)}^{\mathit{\eta }}.\end{array}$$

The advantage of Eq. (2), as compared to Eq. (1), is to separate the dependency on T (return period) from the dependency on D (duration): a(T) only depends on T, and b_Koutso(D) only depends on D. A consequence is that for the particular case of $D_{\mathrm{0}}=\mathrm{1}-\mathit{\theta }$:

$$\begin{array}{}\text{(4)}& i_T\left(D_{\mathrm{0}}\right)=a\left(T\right).\end{array}$$

Then, it becomes a classical frequency analysis of the random variable I(D₀) to estimate the return levels i_T(D₀) – i.e., evaluate ℙ[I(D₀)≤i(D₀)]. Then, Eq. (2) can be reformulated as an equality of distribution of random variables I:

$$\begin{array}{}\text{(5)}& I\left(D\right)\stackrel{d}{=}I\left(D_{\mathrm{0}}\right)\times b_{\mathrm{Koutso}}\left(D\right).\end{array}$$

3.1.3 Simple scaling relationship

In the particular case of θ=0, Eq. (5) becomes

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }I\left(D\right)\stackrel{d}{=}I\left(D_{\mathrm{0}}\right)\times b_{\mathrm{SiSca}}\left(D\right),\text{(7)}& {\displaystyle }& {\displaystyle }{b}_{\mathrm{SiSca}}\left(D\right)=D^{\mathit{\eta }},\end{array}$$

where b_SiSca(D) is a simple scaling formulation of b.

3.2 IDF scaling formulations in the frame of the extreme value theory

In the scaling approach described above, the estimation of rainfall return levels requires a statistical model of rainfall intensity distribution since Eqs. (5) and (6) take the form of an equality of distributions. The extreme value theory (EVT, Coles, 2001) is the most commonly used framework for deriving these models.

3.2.1 Block maxima framework in extreme value theory

EVT proposes two methods to extract samples of extreme values from a time series (Coles, 2001): block maxima analysis (BMA), which consists of defining blocks of equal lengths (often one year in hydrology) and extracting the maximum value within each block; and peak over threshold (POT), which consists of extracting all the values exceeding a given threshold.

Figure 8 Return level plot for the 1 h duration – GEV(1 h) – at two emblematic stations over Senegal: Dakar and Ziguinchor. The estimates come from the GEV estimates at the 24 h duration – GEV(24 h) – downscaled to the 1 h duration (dashed line). Shaded areas represent the 90 % confidence interval of uncertainty due to the fitting of GEV(24 h) and of uncertainty due to the downscaling. Obtained GEV(1 h) from the fitted IDF_SiSca at these stations are also displayed (solid line). The circles represent the observed annual maxima intensities i(D).

Figure 9 Evolution of the spread of the 90 % confidence interval of return levels depending on the return period. The red color of the spread is due to the uncertainty of GEV(24 h) fitting, while the blue color of the spread is due to the uncertainty of the scaling. Results for all stations have been combined. The mean, min, max, and 50 % confidence interval of the spread obtained at the different stations are also shown.

Figure 10 Maps of IDF_SiSca parameters: μ₀ [mm h⁻¹], σ₀ [mm h⁻¹], ξ [–], and η [–].

Figure 11 Maps of return levels intensities (from IDF_SiSca) for different return periods (T=2, 10, and 100 years a–c) for the reference duration (1 h): i_T(D=1) [mm h⁻¹].

Compared to BMA, POT has the advantage of allowing the selection of more than one value per year, thus increasing the sample size used for inferring the model. However, the choice of an appropriate threshold is often difficult (Frigessi et al., 2002). Here the BMA approach was preferred as it is more straightforward to implement.

In BMA, when the block is large enough (which is ensured for annual maxima), the EVT states that the generalized extreme value (GEV) distribution is the appropriate model for block maxima samples (Coles, 2001). The GEV distribution is fully described by three parameters: the location (μ), the scale (σ), and the shape (ξ), which are respectively related to the position, the spread and the asymptotic behavior of the tail of the distribution:

$$\begin{array}{ll}{\displaystyle }{F}_{\mathrm{GEV}}(i;\mathit{\mu },\mathit{\sigma },\mathit{\xi })=& {\displaystyle }\exp \left\{-{\left[\mathrm{1}+\mathit{\xi }\left({\displaystyle \frac{i-\mathit{\mu }}{\mathit{\sigma }}}\right)\right]}^{-\frac{\mathrm{1}}{\mathit{\xi }}}\right\}\\ \text{(8)}& {\displaystyle }& {\displaystyle }\mathrm{for}\mathrm{1}+\mathit{\xi }\left({\displaystyle \frac{i-\mathit{\mu }}{\mathit{\sigma }}}\right)>\mathrm{0}.\end{array}$$

A positive (negative) shape corresponds to a heavy-tailed (bounded in the upper tail) distribution. When ξ tends to 0, the GEV reduces to the Gumbel distribution (light-tailed distribution):

$$\begin{array}{}\text{(9)}& F_{\mathrm{GUM}}(i;\mathit{\mu },\mathit{\sigma })=\exp \left\{-\exp \left[-\left({\displaystyle \frac{i-\mathit{\mu }}{\mathit{\sigma }}}\right)\right]\right\}.\end{array}$$

3.2.2 GEV parameter formulation in a scaling framework

Menabde et al. (1999) have derived the equations merging the scaling formulations presented above (both b_SiSca and b_Koutso) with the extreme value distributions (see also Panthou et al., 2014b; Blanchet et al., 2016). In this approach, the I(D) samples are modeled by a GEV distribution for which the location and scale parameters are parameterized as a function of D as follows:

$$\begin{array}{}\text{(10a)}& {\displaystyle }& {\displaystyle }I\left(D\right)\sim \text{GEV}\left\{\mathit{\mu }\left(D\right);\mathit{\sigma }\left(D\right);\mathit{\xi }\right\},\text{(10b)}& {\displaystyle }& {\displaystyle }\mathit{\mu }\left(D\right)={\mathit{\mu }}_{\mathrm{0}}\times b\left(D\right),\text{(10c)}& {\displaystyle }& {\displaystyle }\mathit{\sigma }\left(D\right)={\mathit{\sigma }}_{\mathrm{0}}\times b\left(D\right).\end{array}$$

The return levels are easily obtained at all durations D as follows:

$$\begin{array}{}\text{(11)}& i_T\left(D\right)=F_{\mathrm{GEV}}^{-\mathrm{1}}\left(D,\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{T}}\right).\end{array}$$

This formulation is equivalent to the following:

$$\begin{array}{}\text{(12)}& i_T\left(D\right)=F_{\mathrm{GEV}}^{-\mathrm{1}}\left(D_{\mathrm{0}},\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{T}}\right)\times b\left(D\right)\end{array}$$

with $D_{\mathrm{0}}=\mathrm{1}-\mathit{\theta }$.

Note that Eqs. (10) to (12) are valid for both b_SiSca and b_Koutso. In log–log space, the IDF_SiSca return levels have a linear shape, indicating a single temporal scaling regime, while those of IDF_Koutso could present a more or less pronounced curvature, indicating a transition between two temporal scaling regimes.

4 Methodology: inference, evaluation, and uncertainty of IDF models

In this study, two IDF models are compared: the IDF_Koutso obtained from the Koutsoyiannis scaling b_Koutso, and the IDF_SiSca obtained from the simple scaling b_SiSca. Both models describe the distribution of extreme rainfall intensities across durations but they differ in their formulation and in their number of parameters: IDF_SiSca has four parameters $\mathit{\{}{\mathit{\mu }}_{\mathrm{0}},{\mathit{\sigma }}_{\mathrm{0}},\mathit{\xi },\mathit{\eta }\mathit{\}}$, while IDF_Koutso has five parameters $\mathit{\{}{\mathit{\mu }}_{\mathrm{0}},{\mathit{\sigma }}_{\mathrm{0}},\mathit{\xi },\mathit{\eta },\mathit{\theta }\mathit{\}}$. BMA samples from which the two scaling models are inferred and evaluated are built by using 1-year block lengths, in order to ensure independence between the elements of the sample. At each station, the extreme rainfall sample thus consists of annual maximum intensities i(D) with D ranging from 1 to 24 h: $\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{6},\mathrm{8},\mathrm{10},\mathrm{12},\mathrm{15},\mathrm{18},\mathrm{24}\mathit{\}}$ h. The lower bound of this range (1 h) was selected in order to limit the risk of underestimating the true annual maximum intensity when evaluating at shorter durations (close to the 5 min fixed window of the raw series). The upper bound of the range (24 h) was chosen because it is a standard duration for hydrological applications and climate studies, but also because it is much more frequently recorded (by daily rain gauges). Note that for each duration D, a rolling mean of length D is applied to the 5 min rainfall series before extracting the maxima. This ensures that the extracted maxima are not underestimated, which is the case when using a fixed window.

4.1 IDF model inference

Different fitting methods have been tested to adjust the IDF model parameters to rainfall data. One of them (the two-step method) is applicable to both IDF_Koutso and IDF_SiSca models.

Note that two other methods specifically dedicated to the IDF_SiSca model were also tested: one based on the moment scaling function (as in e.g., Borga et al., 2005; Nhat et al., 2008; Panthou et al., 2014b), and one based on the global maximum likelihood estimation (as in Blanchet et al., 2016). As they did not perform better than the two-step method, they are not presented here.

The fitting of the scaling b(D) is based on the equality of distribution given in Eq. (5) for IDF_Koutso and Eq. (6) for IDF_SiSca. If these equations hold, the scaled random variables I(D)∕b(D) have the same distribution as the random variable I(D₀) for all durations D. This means that the observed scaled samples i(D)∕b(D) have similar statistical properties for each duration D. Based on this property, the parameters of b(D) are calibrated in order to minimize a statistical distance between the different scaled samples i(D)∕b(D). As suggested by Koutsoyiannis et al. (1998), the difference in medians computed by the Kruskal–Wallis statistic applied on multi-samples (Kruskal and Wallis, 1952) was chosen to characterize the distance between the scaled samples i(D)∕b(D).

Once the scale relationship is identified ($\widehat{b}$), scaled samples $i\left(D\right)/\widehat{b}\left(D\right)$ are computed and pooled in a single sample since they are expected to follow the same GEV distribution (see Eq. (12). The GEV parameters are estimated on this aggregated scaled sample by using the L-moments method. This method was retained as it is more suitable for small samples (Hosking and Wallis, 1997) than the maximum likelihood estimation algorithm, which sometimes fails in optimizing the likelihood for small samples.

4.2 Model evaluation and selection

With the aim of selecting the best IDF model from the two compared IDF formulations (IDF_Koutso and IDF_SiSca), a process of model evaluation and comparison is proposed here by looking at both their flexibility (the models are fitted on a calibration sample) and their robustness (the models are fitted in a predictive mode on samples not used for calibration).

The flexibility characterizes the capacity of a model to fit the observed data that are used to calibrate its parameters. To evaluate flexibility, the IDF models are fitted at each station, then different scores are computed to assess the fitting performances.

The robustness, on the other hand, aims at evaluating whether or not the IDF model is too flexible due to the model having too many parameters with respect to the number of observations. It thus depends on the sensitivity of the IDF model parameters to sampling effects: the less the model parameters are sensitive to sampling effects, the more robust the model. As the two models tested here have a different number of parameters (4 for IDF_SiSca, 5 for IDF_Koutso), there is particular interest in comparing how the goodness of fit for each model is degraded when shifting from the calibration mode to the predictive mode. The predictive capacity of the IDF models is assessed by using a classical calibration/validation process. At each station, a subset of data is used to fit the IDF model; a second independent subset is used to validate it. The same scores used in the calibration mode are computed for the validation subset. Rather than using two consecutive sub-periods, one for the calibration sample and one for the predictive sample, a year-to-year separation was used to build the two subsets. This limits the risk of obtaining samples made of years belonging predominantly to a dry period or to a wet period.

The flexibility and the predictive capacity of the IDF models are quantified based on two types of scores: global and quantile–quantile.

The two global scores used are the statistics returned by two goodness of fit (GOF) tests: Kolmogorov–Smirnov (KS) and Anderson–Darling (AD). Each test computes a statistic based on the differences between a theoretical cumulative distribution function (CDF) and the empirical CDF. The null hypothesis is that the sample is drawn from the fitted model. The test also returns the corresponding p value (probability of type 1 error). The p value is used as an acceptation/rejection criterion by fixing a threshold (here 1, 5, and 10 %). These tests and p values were computed for each rainfall duration at each station.

GOF tests allow for evaluating the entire distribution but do not guarantee that all quantiles are correctly estimated. Thus, as a complement, quantile-based scores are also computed. They characterize the relationship between theoretical (obtained from the fitted CDF) and empirical (obtained from the empirical CDF) quantiles. The root-mean-square error (RMSE), the mean error (ME), and the mean absolute error (MAE) quantile-based scores are computed. The full presentation of these scores can be found in Panthou et al. (2012). A weighted version of these scores is also used in order to assign greater weight to unusual quantiles, as proposed by Begueria and Vicente-Serrano (2006) and also presented in Panthou et al. (2012).

4.3 Uncertainty assessment

From a methodological point of view, the central contribution of this paper is its attempt at quantifying the uncertainty associated with IDF calculation in a scaling framework. This involves two distinct aspects. One is the uncertainty linked to the estimation of the scaling parameters. The other is the uncertainty linked to the inference of the GEV parameters. This second component is especially important to consider when applying a scaling model to a location where only daily rainfall series are available, which is the ultimate purpose of regional IDF models. Indeed, in some regional studies, the scaling parameters will have to be inferred from the very few stations where rainfall is recorded at sub-daily time steps; if they display variations in space, then they will have to be spatially interpolated so as to provide scaling parameter at any location of interest, notably at the location of daily rainfall stations. At these stations, the scaled GEV distribution is thus estimated from the daily observations only, making the inference far less robust than when using a richer scaled sample obtained from observations ranging from 1 h (or less) to 1 day.

Therefore, in the following, the uncertainty assessment at a given location will be addressed separately for the two situations: (i) firstly, when observations at this location are available over a whole range of time steps; (ii) secondly, when only daily observations are available.

4.3.1 Global IDF model uncertainty when multi-timescale samples are available

Confidence intervals for IDF parameters and return levels are estimated using a non-parametric bootstrap (Efron and Tibshirani, 1994). For each station, it consists of fitting IDF curves to bootstrap samples (i(D)_boot) obtained from the original i(D) samples. The entire process consists of four steps:

1. The vector of years is resampled with replacement (Monte-Carlo resampling) until its length equals the length of the original vector.

2. Once a year y is drawn in the bootstrap sample of years and the annual maximum for that year is retained for each duration D in order to build the bootstrap sample of rainfall intensities i(D)_boot. This guarantees the coherence between the samples at different durations.

3. The IDF model is fitted on the bootstrap sample i(D)_boot.

4. The obtained parameters – $\mathit{\{}{\mathit{\mu }}_{\mathrm{0}},{\mathit{\sigma }}_{\mathrm{0}},\mathit{\xi },\mathit{\eta },\mathit{\theta }{\mathit{\}}}_{\mathrm{boot}}$ for IDF_Koutso and $\mathit{\{}{\mathit{\mu }}_{\mathrm{0}},{\mathit{\sigma }}_{\mathrm{0}},\mathit{\xi },\mathit{\eta }{\mathit{\}}}_{\mathrm{boot}}$ for IDF_SiSca – and the associated return level i_T(D)_boot are stored.

These four steps are repeated 1000 times leading to generate 1000 i(D)_boot samples and corresponding vectors of length 1000 for the different IDF parameters and IDF return levels stored in step 4. Confidence intervals are computed on these vectors. It is important to emphasize that these confidence intervals are a measure of the global uncertainty associated with the inference IDF model (uncertainty due to the inference of the scaling relationship and uncertainty generated by the inference of the parameters of the scaled GEV).

4.3.2 Scaling versus GEV related uncertainty when only daily samples are available

When only daily observations are available, the GEV parameters are inferred on the corresponding annual block maxima sample of daily data, which contains far less information that the scaled samples used for fitting a scaled GEV when multi-timescale observations are available. The GEV parameters for the sub-daily time steps are then deduced from the daily GEV parameters using scaling parameters that must be inferred from nearby multi-timescale observations. In some cases this might generate a GEV model that differs significantly from the GEV model that would have been fitted directly on the observations at the proper time steps if they were available. This effect is studied here by assuming that only the daily data were available for fitting the GEV at our 14 stations and by implementing the bootstrap approach in a way that allows separating the uncertainty linked to the GEV parameter inference and the uncertainty linked to the inference of the scaling parameters. Analyzing the uncertainty involves two independent bootstrap resampling processes.

The first bootstrapping method used consists in resampling i(24 h) based on 1000 bootstrap drawings and fitting 1000 GEV(24 h) to these bootstrap samples. These 1000 GEV(24 h) are then downscaled to a target duration D using Eq. (10), yielding 1000 different GEV(D) (only the results obtained for the 1 h duration are presented here). The scaling parameters used to inform Eq. (10) are those computed from the complete multi-timescale samples as explained in Sect. 4.1. This process yields a sample of 1000 GEV at 1 h duration – denoted {GEV(1 h)}_GEV. The dispersion of these 1000 GEV(1 h) is linked to the sole sampling effect underlying the adjustment of the initial GEV(24 h), assuming the scaling parameters to be perfectly known.

In a parallel way, the uncertainty associated with scaling is evaluated by generating 1000 downscaled samples. The reference GEV(24 h) fitted on the original sample i(24 h) is downscaled using 1000 scaling parameters from the bootstrap procedure described in the previous section (Sect. 4.3.1). This produces a sample of 1000 GEV(1 h) denoted {GEV(1 h)}_Scal, whose internal dispersion is only influenced by the uncertainty in inferring the scaling parameters, assuming the reference GEV(24 h) to be perfectly known.

5 Results

5.1 Model evaluation and selection

The model evaluation results are presented in Figs. 4, 5 and Table 1. In these figures and table, the Panels (a) and (b) relate to the calibration and validation results, respectively.

Table 1 Global quantile–quantile scores results for the different IDF models: (a) calibration mode and (b) validation mode. All scores are expressed in mm h⁻¹.

Figure 4 presents the GOF p value of the KS test obtained for both models (IDF_SiSca and IDF_Koutso) in calibration and validation mode at each station (the AD test gives similar results, not shown). In Fig. 5, all stations are gathered in one single q–q plot from which global scores are computed. All global results (non-weighted and weighted q–q scores) are reported in Table 1.

5.1.1 Flexibility and robustness

Figure 4a shows that for all stations and durations, the KS p values are higher than 10 % (i.e., the risk of being wrong by rejecting the null hypothesis “observations are drawn from the models” is greater than 10 %). This means that both IDF models fit the observed data with a reasonable level of confidence in calibration and have thus good flexibility skills. The global scores reported in Fig. 5a and Table 1 show that in calibration, IDF_Koutso slightly outperforms IDF_SiSca. This result was expected as IDF_Koutso has an additional degree of freedom (θ parameter) compared to IDF_SiSca.

With regards to the validation mode, four stations display p values below 10 % at almost each duration (Fig. 4b); globally, both models display a similar number of occurrences of p values below 10 % (37 for IDF_SiSca and 35 for IDF_Koutso) as well as below 5 % (21 for IDF_SiSca and 20 for IDF_Koutso) and below 1 % (1 for IDF_SiSca and 2 for IDF_Koutso).

The global q–q plots in Fig. 5 and the statistics summarized in Table 1 confirm that the two IDF models perform very similarly in validation. IDF_SiSca has slightly smaller biases (mean errors) while RMSE and MAE are slightly better for IDF_Koutso.

5.1.2 Model selection

In addition to performing closely to each other in both calibration and validation modes, the two models yield very similar parameters and return levels, as may be seen from Fig. 6. It is worth noting that the fitted values of the additional parameter θ of the IDF_Koutso model range from −0.02 to 0.39, which is relatively close to zero compared to the [1–24 h] range of durations considered here. This means that the IDF_Koutso model is de facto very close to the IDF_SiSca model, which is a simplification of the IDF_Koutso model assuming θ being equal to zero.

Consequently, while there is no factual reason for considering one of the models to be better than the other, the IDF_SiSca model will be retained, according to the following considerations.

1. The model is more parsimonious, with no clear advantage brought by the fifth parameter of the IDF_Koutso model.

2. The model is easier to implement, especially from the perspective of regional studies involving the mapping of the scaling parameters.

3. The model a straightforward link between the formulation of the IDF_SiSca model and that of the Montana formula (see Appendix) commonly used in national or regional agencies; this makes the formulation of the final IDF product easier to grasp by end-users, thus facilitating its adoption and use.

5.2 Assessing uncertainties

5.2.1 Global IDF model uncertainty when multi-timescale samples are available

The bootstrap approach presented in Sect. 4.3.1 yields confidence intervals representing the global uncertainty linked to sampling in a situation where several samples at different time steps are available at the same location. More precisely, it makes a Monte-Carlo exploration of how the aggregated scaled sample built from the multi-timescale initial samples may vary depending on the random variations of each initial sample. The results are presented in Fig. 7 for four major cities spread over Senegal. Three of them have all of their GOF p values above 0.1, in both calibration and validation modes (Fig. 4), while the fourth (Dakar) has its GOF p values mostly below 0.1 in validation mode, a few of them being even below 0.05 (meaning that, at that particular station, the model is less skillful).

The 90 % confidence intervals of the IDF curves are displayed as colored stripes in Fig. 7. As intuitively expected for a given station, the higher the return periods considered, the larger the confidence intervals. Equally in agreement with knowledge and practice is the fact that, for a given parameter, the largest uncertainty intervals are usually obtained for the shortest series (Fatick, Podor, and Thies), while the longest series (Dakar-Yoff, Tambacounda, Kaolack, and Ziguinchor) display the narrowest intervals (Table 2). However, this relation weakens when considering higher moments or higher return periods: the coefficients of correlation between the confidence interval width and the sample size (available number of years) are r²=0.80 for μ, r²=0.88 for σ, r²=0.69 for ξ, r²=0.55 for i_T=2 and i_T=10, and r²=0.04 only for i_T=100. The presence of very rare events in an observed sample is another factor widening the confidence intervals because some bootstrap samples will include these values, while others will not.

Table 2 IDF_SiSca fitted parameter values and 90 % confidence interval estimated by bootstrap (in brackets).

N corresponds to the number of available years, thus the number of annual maxima.

When comparing the confidence intervals computed for each parameter of the scaled GEV, it appears that their width is well correlated between μ and σ (r²=0.82) and much less so between μ or σ and ξ (r²=0.32 between σ and ξ). The widths of the confidence intervals are quite large for both σ and ξ, which was expected since 8 stations out of 14 have a sample size smaller than 30. The uncertainty on ξ is a sensitive issue, since the confidence interval may include negative values, implying a bounded behavior (Weibull domain of attraction), whereas a light (zero shape value – Gumbel domain of attraction) or heavy (positive shape value – Fréchet domain of attraction) behavior is usually expected for rainfall extremes. However, it is slightly positive on average (+0.046), which is in agreement with the results obtained by Panthou et al. (2012, 2013, 2014b) that point to a dominantly heavy-tailed behavior in the central Sahel region.

5.2.2 Scaling versus GEV related uncertainty when only daily samples are available

As previously explained, at stations where only daily data are available, the sub-daily GEV distributions have to be estimated from this limited set of 24 h values. This significantly increases the uncertainty, as seen in Fig. 8. In this figure, the total uncertainty on the 1 h GEV distribution is divided into (i) the uncertainty linked to the initial fitting of the 24 h distribution – GEV(24 h) uncertainty – and (ii) the uncertainty generated by using the scaling relationships of Eqs. (10b) and (10c) in order to downscale to 1 h distribution GEV(1 h) – scaling uncertainty. This decomposition is carried out by following the procedure presented in Sect. 4.3.2. The results are given for the two longest series of our data set (Dakar Yoff, 38 years; Ziguinchor, 44 years), which display two different behaviors. At Dakar Yoff, the GEV(24 h) uncertainty becomes distinctly larger than the scaling uncertainty from the 10-year return period onwards; at Ziguinchor, this occurs only from the 100-year return period onwards. Associated with this difference is the fact that the downscaled GEV model (dots in Fig. 8) diverges from the reference scaled model (continuous line in Fig. 8) for Dakar Yoff while they are almost identical for Ziguinchor. At Dakar, the width of the 90 % confidence interval associated with the estimation of GEV(24 h) reaches 130 mm h⁻¹ for a return period of 500 years, compared to 30 mm h⁻¹ for the confidence interval associated with the scaling uncertainty. At Ziguinchor the values are respectively 50 and 20 mm h⁻¹.

Figure 9 synthesizes the results obtained at all stations, essentially confirming that the inference of the daily scale GEV(24 h) is a far more important source of uncertainty than the inference of the scaling relationship when it comes to estimating the GEV(1 h). Figure 9 displays the minimum, mean, and maximum uncertainty spread obtained on the 14 stations for GEV(24 h) (red) and the scaling relationship (blue); the 50 % shaded interval contains the seven central values. In order to be able to compute these spreads, the values are expressed as a percentage of the rainfall value given by the GEV(1 h) for each station at a given return level. It was found that the spreads due to the GEV(24 h) fit using daily samples are 3 to 4 times larger than those due to the scaling estimate for the 100-year return level and 5 times larger for the 500-year return level.

5.3 IDF products

5.3.1 IDF curves

A typical representation of IDF curves is given in Fig. 7. As a result of the IDF model formulations and the fitting on a unique scaled sample (for both IDF_Koutso and IDF_SiSca), the return level curves are parallel (they do not cross) and the intensities decrease as the duration increases. The log–log linearity between return levels and durations comes from the simple scaling formulation (the curves would be bent but still parallel for the IDF_Koutso model). Rainfall return levels are of a similar order of magnitude for the four stations, although a north–south gradient is apparent, with rainfall intensities gradually increasing from Saint-Louis to Dakar and from Dakar to Ziguinchor. At the 2-year return period, rainfall intensities vary from roughly 40 mm h⁻¹ (between 33 and 60 mm h⁻¹ when considering all 14 stations) for the 1 h duration to approximately 3 mm h⁻¹ (between 2 and 5 mm h⁻¹) for the 24 h duration. For any station, the return levels for the 10-year and 100-year return periods are approximately 1.5 and 2 times higher than the 2-year return levels, respectively; these ratios hold at all timescales (from 1 to 24 h) as a result of the log–log linearity of the intensity versus the duration. As already discussed in Sect. 5.2, the novelty of these IDF curves is the fact that they are provided with their confidence intervals, allowing the user to get a representation of the uncertainty surrounding the estimated intensity return levels, which is linked to both the sample size and by the quality of the whole GEV and scaling model.

5.3.2 IDF mapping for Senegal

Maps of the four IDF parameters (GEV + scaling) over all of Senegal are plotted in Fig. 10. They have been produced by kriging the parameters inferred at each of our 14 stations. Two of these parameters (ξ and η) are independent of the duration D, while μ and σ are functions of D; these two parameters are thus mapped for the reference duration of 1 h only (corresponding thus to μ₀ and σ₀). They both display a clear north–south increasing gradient, a feature already found by Panthou et al. (2012) for the central Sahel: the location and scale parameter ranges from around 30 mm h⁻¹ and 10 mm h⁻¹ in the north, respectively to around 50 mm h⁻¹ and 15 mm h⁻¹ in the south, respectively. While there are different factors that may explain this gradient, it is clearly coherent with the similar gradient of the mean number of wet days (Fig. 2) that makes the occurrence of a rainfall intensity less frequent in the north than in the south, simply because there are fewer rainfall events there (as evidenced for the whole region by Le Barbé et al., 2002).

Regarding the two non-duration dependent parameters (ξ and η), the shape parameter ξ does not display any clear spatial organization while the scaling parameter η displays a southwest–northeast gradient (with values ranging from −0.8 to −0.9). This suggests that, in addition to the latitudinal effect, the distance from the ocean might also influence the temporal structure of rainfall events. The values of the scaling parameter are very close to those observed by Panthou et al. (2014b) over the AMMA-CATCH Niger network located near Niamey.

The general pattern of the maps of 2-year and 10-year return levels given in Fig. 11 is almost entirely driven by the north–south rainfall gradient. The pattern of the 100-year return level is a bit less regular, with the distance to the ocean seeming to play a role in the western part of the country and a higher patchiness that is certainly largely due to the sampling uncertainty at such a low frequency of occurrence.

6 Conclusions and discussion

6.1 Main results

This study of extreme rainfall over Senegal for durations ranging from 1 to 24 h confirms previous research reporting that a single temporal regime of scale invariance (simple scaling) seems to hold in tropical Africa for this range of timescales (Panthou et al., 2014b; Agbazo et al., 2016; Ghanmi et al., 2016). Whether this range could be extended to sub-hourly and/or sup-daily rainfall intensities is an open research question, out of the scope of this paper, but that can be addressed using the recent methodology developed in Innocenti et al. (2017). The simplified GEV and scaling formulation proposed by Panthou et al. (2014b) with four parameters (three for the GEV and one for the scaling) performs similarly to the five parameter formulation of Koutsoyiannis et al. (1998). This simplified formulation permits an easier study of the sampling uncertainties associated with the inference of the 4 parameters, carried out by a bootstrap resampling in the observed sample of extreme rainfall at 14 stations. Thus in addition to more solidly establishing that scaling is an appropriate hypothesis for this region of the world, our study provides for the first time a comprehensive assessment of the different uncertainties affecting the IDF curves produced by the model (other studies dealing with uncertainty focus on the whole IDF uncertainty, such as Mélèse et al., 2017).

The key advantage of the GEV and scaling approach for computing IDF curves is twofold: (i) it ensures timescale coherency (for the range of explored durations) when working at a regional scale, thus allowing for a coherent spatial interpolation of the IDF model parameters over the region of interest; and (ii) it offers the possibility of deducing GEV distributions for shorter durations at locations where only 24 h data are available, thanks to this spatial interpolation. Both properties have been used in this paper. First, a one at a time simulation approach was used to explore the partition of the overall uncertainty between the GEV inference uncertainty and the scaling model inference uncertainty. One important result in this respect is that the uncertainty produced by the inference of the GEV parameters is 3 to 4 times larger than the uncertainty associated with the inference of the scaling relationship. This means that the scaling relationship requires far fewer data to be inferred correctly than the GEV model. Secondly, maps of the four IDF model parameters and associated intensity return levels have been computed, allowing for the retrieval of the general spatial pattern of these parameters over Senegal. The location (μ) and scale (σ) parameters of the GEV distribution, as well as the rainfall intensity levels for the 2-year and 10-year return periods, display a clear increasing gradient from north to south in line with the climatological gradient of the mean annual rainfall and of the occurrence of wet days. By contrast, for the temporal scaling parameter η, the increasing gradient is rather oriented from northeast to southwest, reflecting the influence of both the occurrence of wet days and of the distance to the ocean. The map of ξ is somewhat patchy, reflecting the fact that this parameter is usually difficult to estimate, but another important result of this study is that its average value is slightly positive, suggesting that the rainfall distribution is heavy-tailed as often observed in several regions in the world (Koutsoyiannis, 2004a, b). Also worth noting is the fact that the value of η is close to −1 (ranging from roughly −0.9 and −0.8) indicating a steep decrease of intensities as the duration increases. This is a common feature of short and intense rainfalls such as those produced by convective storms. These values are comparable to those found by Mohymont et al. (2004) in the tropical area of central Africa, and to those obtained in the Sahelian region of Niamey by Panthou et al. (2014b), close to −0.9 in both cases. They are larger in absolute value than those found in mid-latitude regions, as already underlined by Van-De-Vyver and Demarée (2010).

A final consideration relates to the implementation of such IDF models in operational services. While the theoretical framework of coupling the GEV and scaling models might be considered difficult to handle outside the world of academic research, implementing them for producing IDF curves is relatively easy, especially when using the simplified approach tested here. This approach has the additional advantage of producing relationships between rainfall return levels that are formally equivalent to the so-called Montana relationship (see Appendix), widely used in operational services, facilitating the implementation and usage of our IDF model in meteorological/climatological services and hydrological agencies.

6.2 Points of discussion and perspectives

In the perspective of extending this work to other tropical regions of the world where sub-daily rainfall data might be rare, it remains to explore the effect of using a fixed window to extract the daily rainfall annual maxima, whereas a moving window was used for all durations (including 24 h) in this study. As a matter of fact, daily records of rainfall are carried out at a given hour of the day (usually 06:00 GMT or local time), producing smaller totals than when a mobile window is used to extract the daily rainfall maximum of a given year. Since the scaling relationships that are used to deduce sub-daily statistics from these fixed-window 24 h maxima are tuned on multi-temporal maxima extracted using mobile windows, there is a potential underestimation bias of the sub-daily statistics inferred at 24 h stations that merits further study.

Another critical question relates to using statistical inferences that presuppose stationarity in time in a context of a changing climate. Warming is already attested in the Sahel and is bound to increase, involving possible changes in annual rainfall patterns induced by changes in the positioning of the Bermuda–Azores High and of the Saharan Heat Low. Indeed, rainfall intensification in this region has already been reported by Panthou et al. (2014a) and by Taylor et al. (2017), likely in connection with an average regional warming of about 0.18 K decade⁻¹ over the past 60 years. While dealing with this question was far beyond the scope of this paper, it is a major challenge for both end-users and researchers. It requires developing non-stationary IDF curves, one possible solution in this respect being to use both long historical rainfall series and the information that can be extracted from future climate model projections (see, e.g.,  Cheng and AghaKouchak, 2014).

At the same time it is important to emphasize that stationarity is an elusive concept whose reality is never guaranteed in nature, even without climate change. The Sahelian rainfall regime, for instance, is known for its strong decadal variability (Le Barbé et al., 2002) with potentially great impacts on most extreme rainfall events (Panthou et al., 2013). The use of long (multi-decadal) rainfall series to fit IDF curves can thus reduce the sampling effects and reduce the IDF uncertainties but they can also introduce some hidden biases linked to this decadal-scale non-stationarity. This happened with the dams built on the Volta River in the 1970s. The dams were dimensioned based on the rainfall information of the previous three decades, which included two abnormally wet decades. The reservoirs never filled up in the 1980s and 1990s. Therefore, while IDF curves are intended to be disseminated to a large community of end-users, users must be warned that they are nothing other than a decision-making support tool to be used with care and to be updated regularly.

Appendix A: Simple scaling IDF to Montana IDF

The IDF Montana formulation is as follows:

$$\begin{array}{}\text{(A1)}& i_T\left(D\right)=a\left(T\right)\times D^{b_{\mathrm{m}}}.\end{array}$$

The underscript _m is used to differentiate the Montana formulation from the scaling expression b in the main paper (_m stands for Montana). In our case, the scaling function is simple scaling (Eq. 12), thus Eq. (A1) becomes

$$\begin{array}{}\text{(A2)}& i_T\left(D\right)=F_{\mathrm{GEV}}^{-\mathrm{1}}\left(D_{\mathrm{0}},\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{T}}\right)\times D^{\mathit{\eta }}.\end{array}$$

The two Montana parameters a and b_m can be derived by using the equality between the two formulations:

$$\begin{array}{}\text{(A3)}& {\displaystyle }a\left(T\right)=& {\displaystyle }{F}_{\mathrm{GEV}}^{-\mathrm{1}}\left(D_{\mathrm{0}},\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{T}}\right),\text{(A4)}& {\displaystyle }{b}_{\mathrm{m}}\left(T\right)=& {\displaystyle }\mathit{\eta }.\end{array}$$

Note that when the simple scaling is verified then: (i) D₀ is equal to 1, and depends only on the unit chosen to expressed the intensity of rainfall; and (ii) the assumption of the dependence of b_m on the return period T in the Montana formulation is no longer valid (b_m is equal to η for all return periods).