Evaluating the Performance of a Max-Stable Process for Estimating Intensity-Duration-Frequency Curves

Abstract

To explicitly account for asymptotic dependence between rainfall intensity maxima of different accumulation duration, a recent development for estimating Intensity-Duration-Frequency (IDF) curves involves the use of a max-stable process. In our study, we aimed to estimate the impact on the performance of the return levels resulting from an IDF model that accounts for such asymptotical dependence. To investigate this impact, we compared the performance of the return level estimates of two IDF models using the quantile skill index (QSI). One IDF model is based on a max-stable process assuming asymptotic dependence; the other is a simplified (or reduced) duration-dependent GEV model assuming asymptotic independence. The resulting QSI shows that the overall performance of the two models is very similar, with the max-stable model slightly outperforming the other model for short durations ($d\le 10\mathrm{h}$). From a simulation study, we conclude that max-stable processes are worth considering for IDF curve estimation when focusing on short durations if the model’s asymptotic dependence can be assumed to be properly captured.

1. Introduction

Much research has been recently done on the application of multivariate methods to estimate Intensity-Duration-Frequency (IDF) curves. IDF curves are a popular tool among hydrologists to estimate exceedance probabilities of extreme rainfall events with different durations. In broad terms, IDF curves model a relationship between intensities of extreme rainfall events and their frequencies (i.e., return periods) as a function of event duration. A challenge in estimating IDF curves is how to deal with the simultaneous modeling of intensities for different durations, in particular, how to account for the possible dependence that could arise between intensities of different durations.

Initially, the conventional approach to model IDF curves was based on univariate extreme value theory (EVT) models. Early work on the topic estimates extreme value distributions individually for several fixed durations and subsequently fits an empirical relation to quantiles (return levels) as a function of duration [1,2,3,4]. This approach is prone to inconsistencies as the natural ordering of quantiles is not guaranteed to be preserved over all durations (in other words: quantiles cross). To address this problem, Koutsoyiannis et al. [5] suggested a consistent extreme value model for intensities as a function of duration with location and scale being functions of duration. Later on, a couple of studies implemented methods from Bayesian statistics for the univariate relationship between intensity and duration. Lehmann et al. [6] formulated a Bayesian Hierarchical Model (BHM) based on the model from [5], and Van de Vyver [7] presented a multi-scale model using Bayesian inference. More recently, Ritschel et al. [8] used this model to characterize stochastic precipitation models, and Ulrich et al. [9] proposed the addition of spatial covariates to the model from Koutsoyiannis et al. [5], extending the work of Fischer et al. [10] using spatial covariates to model daily precipitation maxima. These studies make an assumption of stationarity, which may not be valid under a changing climate. Some recent studies have focused on tackling this issue with univariate methods to construct consistent IDF curves in a nonstationary setting. Some examples include the work of Padulano et al. [11] using the storm index method, as well as those of Ganguli and Coulibaly [12], Ganguli and Coulibaly [13] comparing the estimates from a stationary and nonstationary method.

Many of the previous univariate models assume that rainfall intensities are independent for different durations, thus simplifying the modeling efforts. However, the estimation of an extreme value distribution as a function of durations brings along the problem of dependence of extremes associated with different durations, as longer duration series are always aggregated from series of shorter durations. An important consequence of this way of aggregating is that there exists no “single” event for any given duration. As an example, the values of the 15 min duration events are not events that lasted exactly 15 min, but rather the “largest” 15 min long average values from longer events.

In recent years there has been widespread use of multivariate EVT methods for modeling IDF curves, which allow the explicit modeling of dependence structures that could not be captured by the univariate approach. Essentially, a multivariate extreme value distribution (MEVD) is fitted to extreme precipitation data, with the marginal distributions being frequently modeled by a univariate extreme value distribution. Simultaneously, the dependence structure is described with methods such as max-stable processes or copulas. An early example of this was proposed by Muller et al. [14] who, alongside the independence likelihood (analogous to the univariate approach), proposed the use of a so-called trivariate likelihood for three durations: 1, 24, and 72 h. In this study, the dependence between 24 and 72 h was modeled with a bivariate extreme distribution from the logistic family. Afterward, Van de Vyver [15] investigated the use of the trivariate likelihood by calculating the parameters’ uncertainty using a Bayesian approach. He found that, while the resulting posterior distributions from the trivariate likelihood were narrower than the independence likelihood ones, its limitations were too strict to recommend its use over the independence likelihood. An early example of modeling IDF curves using copulas is the study of Singh and Zhang [16], who used a Frank Archimedean copula to estimate the IDF relationship in a bivariate setting.

Some of the most recent advances have been a result of studies that attempted to model IDF curves within a spatial setting. These approaches take advantage of the methods developed for modeling so-called spatial extremes [17]. For example, Stephenson et al. [18] implemented a spatial max-stable process for IDF estimation; they estimated IDF curves in a spatial setting by incorporating a Bayesian Hierarchical Model in every station with a max-stable process. While their approach was able to capture the spatial dependence, they had to limit their scope to assume that the rainfall maxima were independent for different durations. Later on, Tyralis and Langousis [19] proposed the use of a max-stable process to estimate IDF curves for a single station in a way that the asymptotic dependence between rainfall intensity maxima of different durations was explicitly modeled, extending the spatial methodology by proposing a so-called duration space instead of geographical coordinates. Remarkably, their proposed max-stable process was able to explicitly account for the asymptotic dependence between intensities for different durations.

Although Tyralis and Langousis [19] demonstrated the feasibility of a duration-dependent max-stable process to estimate IDF curves, they did not investigate their performance compared to, e.g., univariate EVT methods. Their results showed that both the univariate and multivariate approaches adequately approximated the empirical quantiles, where the max-stable approach resulted in more conservative (i.e., higher) intensities for large quantiles (i.e., longer return periods). In our study, we aimed to build upon the results of Tyralis and Langousis [19] by estimating the impact on performance when accounting for the asymptotical dependence between rainfall intensity aggregated over different durations. We expected that, whenever the asymptotical dependence between durations is high and the max-stable approach is able to capture the “strength” of such dependence for the estimated intensities, the model performance should be significantly higher than for a model assuming independence.

The present paper introduces a scheme to evaluate the impact of the asymptotic dependence between rainfall intensities aggregated over different durations on the performance of EVT-based IDF models. This involves a comparison of skill between two IDF models: one accounting for asymptotic dependence between durations, the other assuming independence. The comparison shows that accounting for the dependence between rainfall intensity aggregated over different durations slightly improves the point estimates for long return periods, and particularly for “short” durations $(d\le 10\mathrm{h})$ usually associated with convective phenomena. However, this comes at the price of increased complexity of modeling the asymptotic dependence.

2. Methods and Data

Our study involves numerical experiments to estimate the relative performance of IDF curves modeled with two approaches: one based on a max-stable process to account for asymptotic dependence (henceforth named as the MS-GEV approach), and another one based on the assumption of independence using the reduced d-GEV model (henceforth named as the rd-GEV approach). By doing so, we aimed to estimate how the performance of IDF curves is affected by considering (or ignoring) the asymptotic dependence between rainfall intensity maxima for different durations.

We performed the study in two broad steps. First, we conducted a case study using data from 6 rain-gauge stations. This data was used to estimate the respective IDF-model parameters from both the MS-GEV and rd-GEV approaches. We compared the performance of both approaches using a measure of skill. In the second step, we introduced synthetic data with known levels of dependence to estimate the effect that the level of dependence has on the resulting estimations. We used the synthetic data to estimate and compare their performance again. Finally, we compared the results from both steps to determine how the performance was affected by the asymptotic dependence between durations. The two different methods used for estimating IDF curves are explained in more detail in the following section. Subsequently, the methods used for verification are described, and, finally, the observation data and the synthetic data are presented.

2.1. Estimation of IDF Curves

Let ${\zeta }_d\left(t\right)$ be the instantaneous rainfall intensity series integrated over a time window of length d, where d is an arbitrary time duration (which commonly ranges from the measurement interval to 72–120 h). Given ${\zeta }_d\left(t\right)$, we obtain the series of the maximum annual average rainfall intensity for each value of d as

$$i\left(d\right)=\underset{y^{-}<t<y^{+}}{max}\left\{{\zeta }_d\left(t\right)\right\}y=(1,\dots ,n),$$

(1)

where n is the total number of observation years, and $y^{-},y^{+}$ are the beginning and end of the yth year, respectively. As a rule, k durations $d_j,j=1,\dots ,k$ are simultaneously used when constructing $i\left(d\right)$, with values from the measurement interval to 120 h (depending on the application). The resulting k series can be thought of as a realization of a random variable $I\left(d\right)$. Notice that, in Equation (1), d is not a random variable but a parameter for the intensity, as noted by Koutsoyiannis et al. [5].

The construction of $i\left(d\right)$ from a single duration series (e.g., hourly precipitation sums) generates a statistical dependence between $i\left(d_1\right)$ and $i\left(d_2\right)$ corresponding to different aggregation durations $d_1\ne d_2$. For example, $i(d=2\mathrm{h})$ and $i(d=3\mathrm{h})$ show a very high dependence (that is, when one of them has what is considered to be a high value, the other intensity also has a high value). However, as the gap in aggregation duration between the values grows, this dependence diminishes: $i(d=2)$ and $i(d=24)$ show almost complete independence. Nadarajah et al. [20] proposed a scheme to work with this type of random variable, which they denoted as ordered random variables. Some authors have linked this concept to the different physical processes that result in different time scales for precipitation events. For example, Muller et al. [14] claimed that the independence between the 1 h and 24-h events were due to the 1-h events being a result of convective (local) motions, while the 24-h event was related to synoptic phenomena. In this paper, we take a closer look at how this dependence affects the estimation of IDF curves down the road.

Following the block-maxima approach, e.g., [21], a GEV can be fitted to each k-series of $i\left(d_j\right)$. Then, the return levels $z_{d,T}$ associated with return periods $T=1/p$ can be calculated for each duration used for the fit $d_j$, where p is denoting the non-exceedance probability with values usually in the range corresponding to upper extreme quantiles. The idea behind IDF curves is to describe the return level $z_{d,T}$ for arbitrary durations d based on the sample $i\left(d_j\right)$ in a meaningful way. This typically involves a parametric form of $z_{d,T}$ as a function of d. The choice of the model used for estimating IDF curves depends on the choice of parametric form of $z_{d,T}$.

For this study, we employed two different approaches for the parametric form of $z_{d,T}$ as a function of d. For the model that assumes asymptotical independence for $i\left(d\right)$ for different d, we follow the duration dependent GEV model (d-GEV) of Ritschel et al. [8]. Based on Koutsoyiannis et al. [5], the d-GEV model of Reference Ritschel et al. [8] estimates a GEV simultaneously from annual maxima associated with various durations, thus conceiving the GEV as a function of duration. The d-GEV yields consistent quantiles $z_{d,T}$, which cannot cross by definition. The other approach is based on a max-stable process for modeling the relationship of $z_{d,T}$ that accounts for the asymptotic dependence between durations.

2.1.1. Using the Duration-Dependent GEV

Following Ritschel et al. [8] and Ulrich et al. [9], we used the duration dependent GEV (d-GEV) to model $i\left(d\right)$ with the distribution

$$G\left(x\right)=exp\left[-{\left(1+\xi \left(\frac{x}{\sigma \left(d\right)}-\tilde{\mu }\right)\right)}^{-1/\xi }\right],$$

(2)

where $\sigma \left(d\right)={\sigma }_0/{(d+\nu )}^{\eta }$ is the duration dependent scale parameter, and $\tilde{\mu }=\mu \left(d\right)/\sigma \left(d\right)$ is the modified location parameter. Given the estimated $\left(\mu \right(d),\sigma (d),\xi )$ parameters, it is straightforward to calculate the return level $z_{d,T}$ for any arbitrary duration using

$$z_{d,T}=\mu \left(d\right)+\frac{\sigma \left(d\right)}{\xi }\left[{\left(-log\left(1-\frac{1}{T}\right)\right)}^{-\xi }-1\right].$$

(3)

To compare the resulting $z_{d,T}$ of this approach with those estimated using the MS-GEV approach, we set the parameter $\nu =0$. This parameter is related to sub-hourly duration values ($0<d<1$), which we do not consider here. Therefore, the dependence of location and scale parameter on duration follows

$$\begin{array}{c}\hfill \mu \left(d\right)=\tilde{\mu }{\sigma }_0{d}^{-\eta },\end{array}$$

(4)

$$\begin{array}{c}\hfill \sigma \left(d\right)={\sigma }_0{d}^{-\eta }.\end{array}$$

(5)

This results in a model with four parameters to be estimated, namely $\{\tilde{\mu },{\sigma }_0,\xi ,\eta \}$. We call this modified distribution the reduced d-GEV or rd-GEV for short. The parameters of the d-GEV distribution are estimated by maximizing the likelihood as implemented in the R-package IDF [22]. Equation (3) is used to get intensities for arbitrary durations $(d\ge 1)$.

2.1.2. Using a Max-Stable Process

Max-stable processes are extensions to infinite dimensions of finite-dimensional extreme value theory models (i.e., extremes of random variables or vectors). They arise as “the pointwise maxima taken over an infinite number of (appropriately rescaled) stochastic processes” [23]. Let $\left\{X\right(x):x\in \chi \}$ be a stochastic process, where $\chi$ is a compact subset of ${\mathbb{R}}^d,d\ge 1$, and $\left\{Z\right(x):x\in \chi \}$ be a max-stable stochastic process. Following de Haan [24], if there exist continuous functions $a_n\left(x\right)>0$ and $b_n\left(x\right)\in \mathbb{R}$, and provided that the limit is non-degenerate, the process $Z\left(x\right)$ can be defined as

$$Z\left(x\right)=\underset{n\to +\infty }{lim}\frac{{max}_{i=1}^n{X}_i\left(x\right)-b_n\left(x\right)}{a_n\left(x\right)},x\in \chi .$$

(6)

The max-stable process $Z(·)$ describes the limiting process of maxima from the $X_i$ IID random fields [25]. The use of this max-stable process for modeling spatial extremes is justified when, based on n independent replicates, and, if n is large enough, we assume that $Z\left(x\right)$ is a good candidate for modeling the partial maxima process $\{{max}_{i=1,\dots ,n}{X}_i\left(x\right):x\in \chi \}$ [26].

One of the main advantages of using a max-stable process is that it provides a flexible way of modeling the dependence structure between the $X_i$ IID random fields. If we assume that $\chi \subset {\mathbb{R}}^2$ represents a geographical catchment, we can think that for multiple points $(x\in \chi )$, the marginal distributions are jointly modeled via the max-stable process, resulting in continuous functions of the GEV parameters $\mu \left(x\right),\sigma \left(x\right),\xi \left(x\right)$ for each margin.

In order to implement a max-stable model for estimating IDF curves, we followed the framework proposed by Tyralis and Langousis [19]. This approach (MS-GEV) employs the Brown-Resnick process, a frequently used parametric family of max-stable processes for modeling environmental extremes [27,28,29]. Previous studies have shown the applicability of the Brown-Resnick process for extreme rainfall applications [17,30,31,32]. A central proposition of our current approach is to define a continuous variable $\overline{i}\left(d\right)$ in a one-dimensional space, where each “location” is one of the durations d. This in contrast to other applications of max-stable processes, where the variable of interest is commonly defined in a two-dimensional (e.g., latitude and longitude) space.

For any max-stable process, the $X\left(s_i\right)$ marginals are generalized extreme value (GEV) distributed, with distribution function:

$$G\left(i\right)=exp\left\{-\left[{\left(1+\xi \frac{i-\mu }{\sigma }\right)}_{+}^{-1/\xi }\right]\right\},$$

(7)

where $\mu ,\sigma ,\xi$ are the location, scale and shape parameters, and $x_{+}=max(0,x)$. Following de Haan [24] (with the adaptation for $d>0$), when the limiting process $\{\overline{i}\left(d\right):d>0\}$ is non-degenerate, a simple max-stable process can be constructed by its so-called spectral characterization, which is a representation of the max-stable process in the frequency domain [23]. In the spectral characterization, the max-stable process is given by choice of the stochastic process ($X_i\left(d\right)$ in Equation (A1)).

Such max-stable processes are called simple as the margins $\overline{z}\left(d\right)$ are unit Fréchet distributed (i.e., $\mu =\sigma =\xi =1$). The use of unit Fréchet marginals is standard, as the max-stable process theory is based on the assumption that the marginals have a common, convenient max-stable distribution. There is no loss of generality in assuming that the limiting process $\{\overline{i}\left(d\right):d>0\}$ has unit Fréchet margins, as it is straightforward to transform such margins into arbitrarily GEV-distributed ones and vice versa.

For this study, we used the bivariate form of the Brown-Resnick process given by Kabluchko et al. [33] (see Appendix A) as the stochastic process $X_i\left(d\right)$. In this form, the dependence is a function of the semivariogram $\gamma$, defined as

$${\gamma }^2\left(h\right)=2{\left(\frac{h}{\rho }\right)}^{\alpha },\rho >0,0<\alpha \le 2,$$

(8)

where $\alpha$ and $\rho$ are, respectively, the smooth and range parameters of the semivariogram, and h represents a measure of the distance between two durations. Tyralis and Langousis [19] calculated this distance as the euclidean distance

$$h_e=|d_j-d_i|,$$

(9)

where the indices i and j denote different durations $d_i$ in hours, and $j\ne i$. However, this measure does not account for the non-linearity of the distance between durations for events of increasing magnitude. For example, consider that an event of 4 h compared to one of 2 h ($h_e=2$) is already twice as long, while an event of 50 h compared to one of 48 h ($h_e=2$) is only 1.04 times the second one.

To address this issue, we explored the use of a distance measure based on a logarithm following Van de Vyver and Van den Bergh [34]. This distance is defined as

$$h_l=lo{g}_2\left(d_j\right)-lo{g}_2\left(d_i\right)=lo{g}_2\left(\frac{d_j}{d_i}\right),$$

(10)

where i and j denote different durations $d_i$ in hours, and $j>i$.

We compare the resulting pairwise extremal coefficients from both distance measures to discern which one results in a more appropiate fit for the semivariogram of Equation (8).

The bivariate form of the Brown-Resnick max-stable process described in Equation (A2) is valid only for unit Fréchet marginals. Therefore, the series of yearly rainfall intensity $\overline{i}\left(d\right)$ requires an appropriate transformation (see [19] for details) to be unit Fréchet distributed, using the relationship

$$\overline{z}\left(d\right)=(1+\xi {(\overline{i}\left(d\right)-\mu \left(d\right)/\sigma \left(d\right))}_{+}^{1/\xi }.$$

(11)

To link $\overline{i}\left(d\right)$ to $\overline{z}\left(d\right)$ for all durations, we used the response surfaces for the GEV parameters [19]:

$$\begin{array}{c}\hfill \mu \left(d\right)={\mu }_0{d}^c,\end{array}$$

(12)

$$\begin{array}{c}\hfill \sigma \left(d\right)={\sigma }_0{d}^c.\end{array}$$

(13)

These response surfaces follow the constraints given by Equations (18)–(20), and describe a function $\mathsf{\Psi }\left(d\right)$ for the parameters of the marginals of $\overline{i}\left(d\right)$. Equations (12) and (13) are equivalent to Equations (4) and (5) in the rd-GEV model. The response surfaces allow to use all durations simultaneously when estimating the max-stable process parameters.

Taking the response surfaces into account, the parameters that we need to estimate for calculating IDF curves using the Brown-Resnick model are six: $[\rho ,\alpha ,{\mu }_0,{\sigma }_0,{\xi }_0,c]$. This is accomplished via the maximum likelihood estimate of the pairwise likelihood given in Equation (A3). The estimation of the parameters is done using the R package SpatialExtremes [35].

Of particular usefulness for our study, is a measure to summarize the strength of the asymptotical dependence modeled by the Brown-Resnick max-stable process. A well-known measure is the extremal coefficient $\theta$, which for the bivariate case, can take values of $1\le \theta \le 2$. The value of $\theta$ decreases when the dependence between the two margins increases. When $\theta =1$ the two margins are completely dependent, and when $\theta =2$ they are independent.

For the dependence structure of the Brown-Resnick max-stable process described by Equation (8), the extremal coefficient is a function only of the distance between durations: $\theta =\theta \left(h\right)$. Given the semivariogram $\gamma$ for a Brown-Resnick max-stable process, with $\mathsf{\Phi }$ denoting the standard normal distribution function, the extremal coefficient is given by

$$\theta {\left(h\right)}_{BR}=2\mathsf{\Phi }{(\gamma \left(h\right)/2)}^{1/2}.$$

(14)

For comparison purposes, we also calculate the extremal coefficient nonparametrically (${\widehat{\theta }}_{emp}$). For our study, we used the method proposed by Marcon et al. [36], which was found by Vettori et al. [37] to generally perform better than other nonparametric estimators.

To summarize: To estimate IDF curves with a max-stable process, we (i) transform our block-maxima data $i\left(d\right)$ into unit Fréchet using Equation (11) with the response surfaces given by Equations (12)–(13); then (ii) estimate the parameters via the maximum likelihood estimates of the pairwise likelihood given by Equation (A3); and, finally, (iii) calculate the intensity for any arbitrary duration d and return period T from Equation (3). We perform all the computations within the R language [38]. The data and code is available as supplementary information.

2.2. Verification and Model Comparison

As a performance measure, we use the quantile score (QS) [39]. This allows us to evaluate predictions of $z_{d,T}$ estimated from IDF-curves in terms of quantiles (i.e., return periods). For D durations, N years, and a given return period T the QS is defined as

$${\mathrm{QS}}_T=\frac{1}{ND}\sum _{d=1}^D\sum _{n=1}^N{\rho }_T(i_n\left(d\right)-z_{d,T}),$$

(15)

where $i_n\left(d\right)$ is the observed block maxima, $z_{d,T}$ is the corresponding intensity from the model (using Equation (3)), and ${\rho }_T\left(u\right)$ is the so-called check function

$${\rho }_T\left(u\right)=\left\{\begin{array}{cc}(1-1/T)u\hfill & u\ge 0\hfill \\ (-1/T)u\hfill & u<0;\hfill \end{array}\right.$$

(16)

thus, for this particular application, ${\rho }_T\left(u\right)={\rho }_T(i_n\left(d\right)-z_{d,T})$.

The QS is always positive and reaches an optimal value at zero. To compare the performance of a model with a reference, we use the Quantile Skill Index (QSI) [9] derived from the Quantile Skill Score $QSS=1-Q{S}_{\mathrm{model}}/Q{S}_{\mathrm{ref}}$ (see [40] for more details on skill scores), defined as

$$\mathrm{QSI}=\left\{\begin{array}{cc}1-\frac{{\mathrm{QS}}_{\mathrm{model}}}{{\mathrm{QS}}_{\mathrm{ref}}},\hfill & if{\mathrm{QS}}_{\mathrm{model}}<{\mathrm{QS}}_{\mathrm{ref}}\hfill \\ -\left(1-\frac{{\mathrm{QS}}_{\mathrm{ref}}}{{\mathrm{QS}}_{\mathrm{model}}}\right),\hfill & if{\mathrm{QS}}_{\mathrm{model}}\ge {\mathrm{QS}}_{\mathrm{ref}}\hfill \end{array}\right..$$

(17)

For our study, ${\mathrm{QS}}_{\mathrm{model}}$ is the score for the MS-GEV approach, and ${\mathrm{QS}}_{\mathrm{reference}}$ is the score for the rd-GEV approach. Positive (negative) values of the QSI indicate a gain (loss) of skill for the MS-GEV approach compared to the rd-GEV one. The advantage of using the QSI over QSS is that negative values have a more meaningful interpretation.

To get a robust estimation of the prediction error for the QSI, we applied 10-fold cross-validation to estimate the QS [41]. The QSI is then obtained using the mean cross-validated QS, averaged over all cross-validation folds, for each model. We used return periods $T=(5,10,20,40,100)$ years. However, the results from the 100 year return period should be interpreted with caution, as the data used for parameter estimation consists of much shorter series of approximately 40 years. The QSI has to be interpreted with care for return periods much larger than the length of the time series (e.g., $T>40\mathrm{years}$). The lack of observations for this region could result in a really high uncertainty of the value of the QS, and, therefore, of the QSI.

2.3. Data

Two different datasets are used for our study. The first one is a synthetic dataset generated for the simulation study. The second one is the block maxima from six rain gauge stations located in the Wupper Catchment (West Germany). We describe both datasets in the following section.

2.3.1. Synthetic Data

We generate synthetic datasets with varying levels of dependence to investigate the models’ performance in estimating IDF curves. We designed three synthetic datasets that simulate rainfall block maxima aggregated over different durations with increasing levels of dependence. For each dataset, we simulate values from a Brown-Resnick simple max-stable process with known dependence parameters using the R-package SpatialExtremes. For the marginal d-GEV distribution, we used a set of parameters characteristic of those d-GEV distributions fitted from the stations in the observational dataset. Then, to fulfill the constraints of annual rainfall maxima averaged over durations $d_i,i=1,\dots ,k$, we transformed the initial simulated data from having unit-Fréchet margins to GEV margins that follow the constraints given by [5,20]

$$\begin{array}{cc}\hfill \zeta \left(d_i\right)={\zeta }_0& ,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \sigma \left(d_i\right)(d_i/d_j)\le \sigma \left(d_j\right)\le \sigma \left(d_i\right)& ,d_i\le d_j\forall i,j,\mathrm{and}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \mu \left(d_i\right)(d_i/d_j)\le \mu \left(d_j\right)\le \mu \left(d_i\right)& ,d_i\le d_j\forall i,j.\hfill \end{array}$$

(20)

This transformation uses the response surface described in Equations (12) and (13). We simulate 40 values for each dataset (representing 40 years) for $d=(1,3,\dots ,119,120)\mathrm{h}$. Following Zheng et al. [25], three sets of dependence parameters were used: $(\rho =1,\alpha =1)$ for weak dependence, $(\rho =0.5,\alpha =0.5)$ for moderate dependence, and $(\rho =0.5,\alpha =0.2)$ for strong dependence. For each parameter set, we generate 1000 realizations. An issue we encountered is that the nature of the simulated data allowed for rainfall intensity series that did not strictly follow the constrains given by Equations (18)–(20), however, we considered this to happen at a frequency that would not affect the final result of the study.

2.3.2. Observations

We used six rain gauge stations from the Wupper Catchment in Germany (Figure 1). All the stations have hourly values of accumulated precipitation height for the period 1979–2016. The stations are all within an elevation range of 250 m, and horizontally the shortest and longest distance between each station is 7 Km and 34 Km, respectively. We chose this dataset as it has a large number of years with high-frequency (hourly) measurements. Furthermore, the stations range from the Bergisches Land to the Upper Rhine Plain and therefore represent very well the different altitudes of the catchment. Figure 1 shows the distribution of rainfall maxima for this period.

We obtain the annual block maxima $i\left(d_j\right)$ of the averaged rainfall intensity ${\zeta }_d\left(t\right)$ over the time window $d_j$ for each station using Equation (1). For estimation purposes, we used durations $d=(1,3,\dots ,119,120)\mathrm{h}$. The decision for the cut-off value of 120 h was based on previous studies on IDF curve estimation [18,19]. By visual inspection of the corresponding Quantile Quantile (QQ)-plot, we ensure good agreement of the resulting $i\left(d_j\right)$ block maxima for all 6 stations with the GEV distribution. A small subset of the QQ-plots can be seen in Figure A2.

3. Results

We present the results for the case study in the Wupper region of Germany first, followed by the results of the simulation study.

3.1. Case Study

3.1.1. Structure of the Extremal Dependence

Figure 2 and Figure 3 show a comparison of the pairwise extremal coefficient $\theta$ derived from the parameters of the MS-GEV approach (Equation (14), red line) and from a nonparametric estimate (dots) to assess how well the MS-GEV approach captures the observed asymptotic dependence. We used different distance measures h for each plot. Figure 2 uses the euclidean distance $h_e$ (Equation (9)), and Figure 3 uses the log-distance $h_l$ (Equation (10)). Additionally, the different colors show the lower distance $d_i$ used for each duration pair $(d_i,d_j)$, where $i<j$, and $(i>0,j>0)$.

For the euclidean distance $h_e$ (Figure 2), the values of ${\theta }_{\mathrm{BR}}$ are close to the binned means of ${\widehat{\theta }}_{emp}$ (black circles) with respect to the scattering of the individual empirical estimates (color circles) for all stations. Nevertheless, the empirical point clouds show a remarkably high variability of $\theta$ around the mean value. To name one example, in station Buchenhofen, the distance of $h_e=25\mathrm{h}$ has a range of very different values of $\theta$, spanning from 1.2 to 1.8. Thus, using a model that approximates the empirical mean-binned values in this case does not appear to be a meaningful representation of the overall variability of the dependence between durations.

Furthermore, each set of duration pairs with a fixed lower duration $d_i$ in Figure 2 (represented by different colors) seems to follow a different path as the distance $h_e$ grows. In particular, for duration pairs with a short lower duration ($d_i\le 10\mathrm{h}$), the value of $\theta$ grows much faster than duration pairs with longer lower durations ($d_i>10\mathrm{h}$). This suggests that several different regimes of dependence coexist, which are not only a function of the distance $h_e$, but also of the magnitude of the durations used. This is a transgression of the assumption that the extremal coefficient, and by extension, the semivariogram model of Equation (8) is isotropic (i.e., $\gamma$ should only be a function of h). It is thus not evident that using a dependence model for the MS-GEV approach based on the euclidean distance $h_e$ adequately captures the asymptotical dependence between durations seen on the data used for this study.

The resulting extremal coefficient using the log-distance $h_l$ is shown in Figure 3. The point cloud shows a remarkably lower variability around the binned means than those of Figure 2. Furthermore, the empirical values of ${\widehat{\theta }}_{emp}$ appear to follow the same regime, suggesting that $\theta$ is isotropic when using the log-distance. The shape of the point cloud seems to be appropriately captured by ${\theta }_{\mathrm{BR}}$ for duration ratios $d_j/d_i\lesssim 6$ (i.e., when the upper duration $d_j$ is around six times the lower duration $d_i$). However, the parametric model deviates from the empirical estimates as $h_l$ increases. With the exception of station Leverkusen, the parametric model consistently overestimates the strength of the asymptotical dependence for the pairs with a duration ratio $d_j/d_i>6$.

In light of the above results shown by Figure 2 and Figure 3, we decided to use the log-distance $h_l$ when estimating the semivariogram of the MS-GEV approach in all of the following calculations.

3.1.2. Estimation of IDF Curves

Table 1. Parameter estimates from the MS-GEV approach for stations in the Wupper catchment using durations $d=(1,3,\dots ,119,120)\mathrm{h}$.

Station	$\mathit{\alpha }$	$\mathit{\rho }$	${\mathit{\mu }}_0$	${\mathit{\sigma }}_0$	${\mathit{\xi }}_0$	c
Bever	1.42	2.09	13.32	2.84	0.03	−0.58
Buchenhofen	1.39	2.22	13.38	3.03	0.02	−0.63
Leverkusen	1.32	2.19	10.74	2.34	0.05	−0.64
Lindscheid	1.54	2.44	13.69	3.23	0.06	−0.64
Neumuehle	1.54	1.84	13.52	2.74	0.04	−0.60
Schwelm	1.53	1.74	13.07	2.77	0.05	−0.62

shows the parameter estimates for the MS-GEV approach for all Wupper catchment stations, estimated from durations $d=(1,3,\dots ,119,120)\mathrm{h}$. The estimates for the range parameter $\rho$ are reasonably consistent across all stations. Their value of $\sim 2$ suggests that, for all stations, the rainfall intensities of different durations become asymptotically independent when their ratio $d_j/d_i$ is larger than ($2^2=4$).

Figure 4 shows the IDF curves following from Equation (3) for the MS-GEV approach (solid lines) and compares them to the IDF curves based on the rd-GEV approach (dashed lines). Return levels from the MS-GEV are for the most part consistently higher, which is in agreement with Tyralis and Langousis [19].

To compare the results of the IDF curves using the euclidean distance $h_e$ instead of the log-distance $h_l$, a plot comparing the resulting 100-year return level of both distances is shown in Appendix A.

3.1.3. Performance Averaged over All Durations

Figure 5 shows the cross-validated QSI evaluating the performance of the MS-GEV approach compared to the rd-GEV one. Similar behavior can be seen for all stations. For short return periods, the QSI is close to zero (denoting that both models are equally good), increasing towards longer periods. Station Lindscheid profits most from the MS-GEV with a 20% increase in skill for the 100-year return level. Only for a few points is skill negative, mostly for the shorter return periods.

3.1.4. Performance for Individual Durations

For a detailed comparison, we show the QSI conditioned on duration in Figure 6. The QSI varies appreciably over different durations for a given return period. For short durations ($d<10\mathrm{h}$), the QSI is mostly positive for all return periods; for long durations $d>100\mathrm{h}$, it is mostly negative. Gauge Lindscheid exhibits positive skill for many more combinations of durations ($d<50\mathrm{h}$) and return periods $T>40\mathrm{years}$; station Bever is the only station not showing a positive skill for very short durations, showing a general loss of skill for all but the shortest durations.

3.2. Simulation Study

We studied the effect of the level of dependence on the performance of the MS-GEV approach. To this end, we used synthetic data with known dependence parameters ($\rho$ and $\alpha$ in Equation (8)) and estimate the performance for various levels of dependence. Figure 7 shows the cross-validated QSI obtained using Equation (17), using the averaged QS values over all $d=(1,2,\dots ,120)\mathrm{h}$ durations for 1000 replications as box-whisker plots. Similar to the case study, the distance h used for the semivariogram of the Brown-Resnick process was the log-distance $h_l$ (Equation (10)).

The variation of skill among the replicates increases with increasing return period. The median is consistently positive (MS-GEV superior to rd-GEV) but below 0.05, suggesting that the dependence does not impact strongly on the performance results. The strong level of dependence leads to a slightly higher median skill but also to a considerable larger variance.

Figure 8 shows the QSI calculated with the model and reference QS averaged over all replicates for individual durations. The results are similar to the pooled QSI over all durations of Figure 7: The QSI increases with return period and dependence, staying below QSI = 0.1 for a 100 year return period and the strongly dependent series. This again suggests that the level of dependence has little impact on the performance of the return levels from the IDF curves.

4. Discussion

In this study, we obtained a measure of relative performance for return level estimates of IDF curves for various durations involving a max-stable process that allows for asymptotic dependence between durations (MS-IDF), compared to a model that assumes independence (rd-GEV). To do so, we built upon the previous study of Tyralis and Langousis [19], who focused on the theoretical basis of using max-stable processes for modeling IDF-curves and did not investigate the consequences in terms of model performance.

To investigate the possible impact of the asymptotical dependence, we evaluated and compared the performance of estimating IDF curves with two different approaches: i) using a max-stable process to describe the dependence between rainfall intensity maxima for different durations and ii) assuming independent maxima. We evaluated individual performance based on a score that allows us to focus on the tail of the distribution, namely the quantile score (QS). The comparison between the MS-GEV and rd-GEV approach was carried out with the quantile score index (QSI), an index based on the QS skill score. The QSI enables us to quantify a gain/loss in performance when accounting for asymptotic dependence in return level estimation.

The results of Figure 2 showed that the resulting pairwise extremal coefficient was non-isotropic when using the euclidean distance $h_e$ as the measure h in the semivariogram of the Brown-Resnick process for the MS-GEV approach. Thus, in contrast to the approach of Tyralis and Langousis [19], we explored the use of a logarithmic distance measure instead of an euclidean one for the semivariogram of the Brown-Resnick process. Figure 3 shows that this was a reasonable choice, with the resulting parametric extremal coefficient properly capturing the variability of the empirical extremal coefficient around its binned means.

A simulation study suggests a minor advantage of the MS-GEV approach over the rd-GEV, particularly for long return periods (large quantiles). This advantage increases with the strength of the dependence, but remains low (QSI $\le 0.1$) even for the strongest level of dependence. A complementary case study for six gauges in the Wupper catchment (Germany) corroborates a general advantage for long return periods when averaging the performance measure over all durations. A detailed investigation of performance conditioned on durations shows, for our case study, that this advantage results mostly from short ($d\lesssim 10\mathrm{h}$) and, in some cases, from intermediate ($20\mathrm{h}\lesssim d\lesssim 50\mathrm{h}$) durations, depending, however, on the specific station.

The presented findings support the idea of Tyralis and Langousis [19] that max-stable processes are valuable models for IDF curve estimation. The simulation and case study results indicate that an increase in skill with the MS-GEV approach is mostly found for large quantiles and short to medium durations. This effect might be related to the fact that shorter durations have a larger number of pairs than the longer durations for the pairwise likelihood. This means that the longer durations could be underrepresented in the current likelihood expression, leading to a better fit of the model for the shorter durations for the MS-GEV approach. However, extreme events of short durations have usually a higher impact on society than those of long durations. Therefore, we believe that this underrepresentation does not necessarily detract value from using the MS-GEV approach. When focusing on large quantiles and short durations, it might be worth taking the added computational expense and model complexity in exchange for increased skill in the return levels for such events.

In addition, the simulation study demonstrated that the level of dependence had a modest impact on the overall performance and variability of the MS-GEV approach. These findings contradict those of Zheng et al. [25], who found that the dependence strength did not influence the performance of the return level estimates. However, our study focused on a different dependence structure than that of Zheng et al. [25], who used a spatial approach, in contrast to our duration space. This contradiction may be associated with the particular form of the asymptotic dependence for the different durations of $i\left(d\right)$, which stems from the nature of $i\left(d\right)$ as random ordered variables. As previously studied by [20,42], when two random variables are intrinsically ordered, each margin’s distribution is affected, something that has to be taken into account. While the response surfaces described in Equations (12) and (13) take this ordering into account, the dependence structure given in Equation (8) does not, a factor that could explain the contradiction with previous studies.

The parametric estimate of the extremal coefficient ${\theta }_{\mathrm{BR}}$ from the MS-GEV approach shown in Figure 3 appears to be a reasonable fit for the empirical values of $\theta$ when the ratio of the duration pair is around ($d_i/d_j\lesssim 6$). However, as the upper duration gets around six times larger than the lower duration, the model consistently overestimates the strength of the dependence. Figure 3 also showed that duration pairs $(d_i,d_j)$ involving short lower durations $d_i\le 10\mathrm{h}$ had an extremal coefficient that approached independence (i.e., $\theta =2$) faster than those with longer lower durations as $h_l$ increased. This suggests that the dependence between rainfall maxima of different aggregation durations is short-ranged, in particular, for the shorter durations.

For an operational use of this approach, uncertainty estimates for the quantiles (IDF curves) need to be incorporated. In this regard, Ganguli and Coulibaly [12], Mélese et al. [43] showed that a Bayesian Hierarchical Model approach resulted in reliable credibility intervals for IDF curves. For the rd-GEV approach, we investigate a bootstrap-based method for estimating the uncertainty of d-GEV based return levels in a different study [9]. We also limited our study to using data with hourly frequency, resulting in the value of the (sub-hourly) $\nu$ parameter of the d-GEV to be artificially set to zero (what we called the rd-GEV). Further studies would benefit from using sub-hourly frequencies, allowing $\nu$ to vary freely. Moreover, we assumed that the stations’ data was stationary, ignoring the possible effects of climate change. Several studies have shown that accounting for nonstationarity has a measurable effect on the return levels estimates [11,12,13].

The method presented in this study is a straightforward and practical manner of estimating IDF return level performance based on a max-stable process. The results for the case study encourage to investigate its performance within a larger geographical setting. Furthermore, it seems worth implementing more flexible functions describing the variability of GEV parameters with duration as used for the response surface, e.g., an additional parameter accounting for different behavior, particularly for short durations, as suggested in Koutsoyiannis et al. [5]. Another critical issue for future studies is to explore how different dependence structures could impact the performance of the estimation. Furthermore, a comparison with the recent developments in the use of covariates for the rd-GEV approach could result in better skill for such an approach compared with our current MS-GEV one [9].

5. Conclusions

Our findings indicate that the use of models that allow for the asymptotic dependence between rainfall maxima of different durations when estimating IDF curves can lead to moderately better return level estimates, particularly for long return periods (100 years, generally of considerable interest) and short durations ($d\le 10\mathrm{h}$). However, the former comes at the expense of the added complexity of modeling the asymptotic dependence. Furthermore, this asymptotical dependence seems to be short-ranged for the short durations. We, therefore, recommend the use of the simpler univariate-EVT methods assuming independence between durations for a single station when the main goal is obtaining return levels for a wide range of short and long durations from IDF-curves.