<figure><img src="image_1.png" /></figure>

HAL
open science

# Integrating non-stationarity and uncertainty in design life levels based on climatological time series

Occitane Barbaux, Philippe Naveau, Nathalie Bertrand, Aurélien Ribes

► To cite this version:

Occitane Barbaux, Philippe Naveau, Nathalie Bertrand, Aurélien Ribes. Integrating non-stationarity and uncertainty in design life levels based on climatological time series. Weather and Climate Extremes, 2025, 50, pp.100807. 10.1016/j.wace.2025.100807. hal-05267687v2

HAL Id: hal-05267687

https://hal.science/hal-05267687v2

Submitted on 17 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Integrating non-stationarity and uncertainty in design life levels based on climatological time series

Occitane Barbaux<sup>a,b,c,●,*</sup>, Philippe Naveau<sup>d,e</sup>, Nathalie Bertrand<sup>b</sup>, Aurélien Ribes<sup>a,c,e,f</sup>

## ARTICLE INFO

Dataset link: https://github.com/Occitane-Barbaux/Integrating-non-stationarity-and-uncertain ty-in-design-life-levels

Keywords:
Uncertainty
Non-stationarity
Risk
Climate change

## ABSTRACT

This work focuses on inferring design life levels for extreme events under non-stationary conditions. Its objectives are twofold. The first one is to provide a single indicator that summarizes relevant and interpretable information about large values in time series, even when stationarity cannot be assumed. Classical risk indicators such as the 100-year return level become difficult to interpret in a non-stationary framework. To address this, we leverage the existing concept of the equivalent reliability (ER) level. Under stationarity, the ER level coincides with the classical return level, but it differs otherwise. More precisely, the ER level ensures that the probability of having all observations below the ER level during a specified design period is controlled. This definition ensures interpretability in terms of safety or failure risk. A second objective is to capture stochastic and estimation uncertainty, a key aspect in any risk analysis, as uncertainties due to inference schemes can grow with extreme intensities. We incorporate both by using the Bayesian predictive distribution. Although well known in Bayesian statistics, the predictive distribution has rarely been applied to climatological time series risk analysis.

Our approach is demonstrated on simulated data and on a case study of annual maxima of temperatures at a site in Southern France. To do so, a non-stationary Bayesian hierarchical extreme value model is used to combine data from 26 CMIP6 general circulation model simulations (SSP2-4.5, 1850-2100) with observations. The resulting predictive ER levels clearly indicate that non-stationarity over a design period of interest, as well as sampling and estimation uncertainty, have to be taken into account for risk assessment. For example, the 1000-year posterior predictive ER level for 2050-2100 is higher than any non-stationary 1000-year return level median estimate over the same period, reflecting the increasing risk due to the non-stationarity of the SSP 2-4.5 pathway.

## 1. Introduction

The effects of climate change on natural hazards have been heavily documented, relying on multiple sources of evidence such as observations and general circulation models (Intergovernmental Panel On Climate Change (Ipcc), 2023). Climatological variables such as temperatures are expected to see their distribution change significantly in the future. This information is crucial for projecting future trajectories and communicating on upcoming risks. At the same time, many components of climatological risk analysis are subject to uncertainty (see, e.g., Hawkins and Sutton, 2009; Lehner et al., 2020). Such uncertainty arises from various sources (for example, measurement

errors and approximations in physical models) and generally increases when the underlying climate processes are non-stationary (see, e.g., Naveau et al., 2020). Various studies have reported and explored the possibility of low-likelihood high-risk scenarios that should be considered in decision making (Intergovernmental Panel On Climate Change (Ipcc), 2023; Hawkins and Sutton, 2009). These include cases where climate change rates may exceed expected ranges or where record-shattering extremes occur (Fischer et al., 2021; Naveau and Thao, 2022). In terms of statistical analysis, specific probability models have been designed to represent uncertainty in records (see, e.g., Robin and Ribes, 2020; Gonzalez et al., 2025). In particular, current and future

<footnote><sup>a</sup> Météo France, 42 Av. Gaspard Coriolis, Toulouse, 31100, France</footnote>

<footnote><sup>b</sup> Autorité de Radioprotection et de Sureté Nucléaire, PSE-ENV/SCAN/BEHRIG, F-92260, Fontenay-aux-Roses, France</footnote>

<footnote><sup>c</sup> Centre National de Recherches Météorologiques, 42 Av. Gaspard Coriolis, Toulouse, 31100, France</footnote>

<footnote><sup>d</sup> Laboratoire des Sciences du Climat et de l'Environnement, Orme des Merisiers, Bat 714, Saint-Aubin, 91190, France</footnote>

<footnote><sup>e</sup> Centre national de la recherche scientifique, 16 Av. Edouard Belin, Toulouse, 31400, France</footnote>

<footnote><sup>f</sup> Université de Toulouse, 41 All. Jules Guesde, Toulouse, 31000, France</footnote>

<footnote>* Corresponding author at: Autorité de Radioprotection et de Sureté Nucléaire, PSE-ENV/SCAN/BEHRIG, F-92260, Fontenay-aux-Roses, France. E-mail address: occitane.barbaux@umr-cnrm.fr (O. Barbaux).</footnote>
warming trends pose challenges not only in terms of constructing complex non-stationary statistical models, but also because the non-stationary elements themselves make it difficult to condense model outputs into simple and comprehensive summaries for risk managers. In practice, risk analysts often seek a single, interpretable number to represent risk. Under stationary assumptions, a common metric is the T-year return level, a level expected to be exceeded on average once every T years. This index's simplicity and interpretability, however, break down in non-stationary settings. A few alternatives for non-stationary cases exist (see, e.g., Rootzén and Katz, 2013), but a consensus has yet to emerge in the statistical climatology literature.

In this context, we define the design period as the time period over which the risk assessment applies. We focus on problems where the decision-making context naturally provides a design period. For example, building codes for critical infrastructure (such as dams or power plants) specify risk levels over the structure's lifetime (Wigley, 2009). In such cases, the design period might be legally set from an initial year $t_1$ to a final year $t_2$. Non-stationarity implies that any relevant risk measure must account for both the lifetime length $t_2 - t_1$ and the initial year.

To address these challenges, we set two objectives. First, we seek a single indicator for risk over a design period that remains interpretable even when stationarity is not assumed. Our approach builds on the existing concept of equivalent reliability (ER) level. In stationary cases, the ER level aligns with the classical return level for the chosen probability of excess. More generally, the ER level is defined to directly bound the risk of any exceedance during the design period. The details will be explained in Section 2.1.

Our second objective is to combine both stochastic and estimation uncertainty into this indicator. The question of how to account for uncertainty in decision-making has received a large amount of attention in the statistical literature. The most common solution relies on confidence or credibility intervals. Their spread depends on the chosen confidence level and can face challenges of over-estimation (Fawcett and Green, 2018) or under-estimation (Zeder et al., 2023). Practically, a high-level decision maker, as already mentioned, may not be comfortable with interpreting confidence intervals and may prefer a single number that combines all uncertainty about the quantity of interest. This is true of construction codes, where a single numerical value is the reference for a large area. For example, Markova et al. (2019) compares reference maps for various natural variables such as maximum air temperature for several European countries. As advocated by Davison and Huser (2015) and Coles and Tawn (2005), the Bayesian predictive distribution naturally combines these uncertainty sources. While predictive distributions are common in Bayesian analysis, they have seldom been used in climatological risk analysis. This framework is presented in Section 2.2. We adapt the predictive framework to ER levels, yielding a “predictive ER level” that combines stochastic and estimation uncertainties. Details are given in Section 2.3

Taken together, these developments allow us to produce a single predictive risk level that reflects non-stationarity and inference uncertainty. This methodology is illustrated in Section 3. The approach is applied to annual maxima of temperatures in southern France fitted with the Bayesian non-stationary GEV (generalized extreme value) model of Robin and Ribes (2020). Ensembles from 26 CMIP6 general circulation models (GCMs) (SSP2-4.5) are combined with the observational record (CRUTEM5 1850–2020 and local meteorological data) at one location in the South of France. While our case study will focus on high temperatures, the presented dual approach of ER and predictive distributions has broad applications, such as in biology and ecology (see, e.g., Pörtner et al., 2006; Asseng et al., 2021; Locosselli et al., 2020).

## 2. Life design levels

### 2.1. Dealing with non-stationary: Equivalent reliability

Before treating non-stationary cases, we begin by reviewing foundational concepts and notation in the simpler setting of independent and identically distributed (i.i.d.) random variables.

Let $Z$ be a real-valued continuous random variable with cumulative distribution function (cdf) $F(z) = P(Z \le z)$. For a return period $T = 1/p$, the return level $z_p$ is the quantile of level $1 - p$ of the cdf:

$$z_p = F^{-1}(1-p). \quad (1)$$

In the i.i.d. case, the return level can be simply interpreted as the level exceeded on average once over the period $T$, and the return period corresponds to the mean interval between exceedances above $z_{1/T}$.

However, in a wide range of applications, and particularly for climatological time series, the assumption of stationarity no longer holds (e.g., due to trends in climate variables). Accordingly, it becomes necessary to explicitly account for temporal dependence by incorporating a temporal index<fnref n="1" /> $t$ into the definition of the random variable, then denoted as $Z_t$. To extend the concept of return levels and periods to non-stationary cases, two natural avenues can be explored, either by fixing the probability $p$ or by imposing a constant level $z$. The former strategy leads to a time-dependent return level defined by

$$z_p(t) = F_t^{-1}(1-p), \quad (2)$$

where $F_t(\cdot)$ is the cumulative distribution function of $Z_t$ at time $t$. This quantity has different names; for example, it is called the ‘effective return level’ in Katz et al. (2002). This provides an annual quantile history and is often used to explain how future levels with a fixed probability risk change with time. In contrast, when communicating on the risk of exceeding a fixed level $z$ (say, a dam height), we can define an instantaneous probability:

$$p_t = 1 - F_t(z).$$

In this case, the instantaneous return period $T_t = \frac{1}{p_t}$ no longer represents the average interval between two events, and $z_p(t)$ loses its interpretation as the level exceeded once on average over the period $T$. This highlights a key challenge in non-stationary settings: the intuitive link between return levels and periods breaks down, as year-to-year conditions can no longer be assumed identical.

To address this, several alternative risk metrics have been proposed (see, e.g., Cooley, 2013). One possibility proposed by Laurent and Parey (2007) was to define the return period as the interval containing, on average, a single exceedance. Another proposition was to view the return period as the average waiting period between two events, leading Salas and Obeysekera (2014) and Olsen et al. (1998) to define the expected waiting time. However, these approaches face similar criticisms; see Rootzén and Katz (2013) and Read and Vogel (2015) for in-depth reviews. They do not address the distinction between the instantaneous return period, as a proxy of the instantaneous risk of excess, and the lifetime period of the plant building. Moreover, for projections over long horizons — such as a 1000-year event with datasets spanning only a few centuries — extrapolation far over the design period is needed (Kim et al., 2020; Yan et al., 2017; Read and Vogel, 2015; Serinaldi, 2015). In some cases, trends in the data can even yield infinite return periods for fixed levels (Cooley, 2013; Yan et al., 2017; Liang et al., 2016). As noted by Read and Vogel (2015), return periods are “not explicitly tied to a planning horizon and thus are unable to characterize the likelihood of an event occurring during a project lifetime”.

<footnote n="1">¹ Throughout this work, the temporal increment is assumed to be one year; however, other time units could be used if appropriate.</footnote>
To overcome these issues, various studies have advocated moving focus from instantaneous risk to design-period-integrated risk measures (Salas and Obeysekera, 2014; Serinaldi, 2015; Yan et al., 2017; Liang et al., 2016). In particular, Rootzén and Katz (2013) highlighted the concept of reliability, defined as the probability that a given level z is not exceeded at any time during the design period starting at time t₁ and ending at time t₂:

$$R_{t_1:t_2}(z) = P[\max(Z_{t_1}, Z_{t_1+1}, \dots, Z_{t_2}) \le z]. \quad (3)$$

Reliability was first studied in a stationary framework (see, e.g., Wigley, 2009) when researchers began to question the use of the average return period as a design metric. Unlike instantaneous measures, reliability (3) is defined over a full period, making it especially suitable for long-term infrastructure planning. The survival function $1 - R_{t_1:t_2}(z)$ simply provides the probability that, given level z, at least one exceedance ($Z_t > z$) occurs during the period $t_1 : t_2$.

In the i.i.d. case, reliability has an explicit form<fnref n="2" />:

$$R_{t_1:t_2}(z_p) = (1-p)^{t_2-t_1+1}, \quad (4)$$

where $z_p$ is the return level of probability p. For example, the reliability over 100 years for a 100-year return level is approximately 0.37, implying a 63% chance of at least one exceedance. This is important information when buildings need to be protected for such a period. The right-hand side of Eq. (4) provides a solid benchmark to link the non-stationary cases with the i.i.d. setup.

Extending this idea to the non-stationary setting, Liang et al. (2016) and Hu et al. (2018) introduced the equivalent reliability (ER) level $z_p^{ER}$ as the level satisfying:

$$R_{t_1:t_2}(z_p^{ER}) = (1-p)^{t_2-t_1+1}. \quad (5)$$

The reliability function $R_{t_1:t_2}(\cdot)$ in the left-hand side of (5) corresponds to (3) but its right-hand side comes from (4), meaning it is set to the same reliability as a sequence of $t_2 - t_1 + 1$ independent years each with an exceedance probability of p. This definition ensures that the interpretation of $z_p^{ER}$ remains consistent with the stationary framework: it is the level exceeded with probability $1 - (1-p)^{t_2-t_1+1}$ over the period from $t_1$ to $t_2$, regardless of underlying non-stationarity. In the i.i.d. (stationary) setting $z_p^{ER}$ is exactly equal, by construction, to the stationary quantile $z_p$. Hence, the interpretation of $z_p^{ER}$ remains simple with the advantage of integrating information over the design period of interest.

Moreover, $z_p^{ER}$ satisfies the inequalities:

$$\min_{t \in t_1:t_2} z_p(t) \le z_p^{ER} \le \max_{t \in t_1:t_2} z_p(t),$$

where $z_p(t)$ was defined by (2). Hence, $z_p^{ER}$ represents a middle ground over all possible values of time-dependent return levels. Applications of ER levels have included extreme precipitation analysis in Germany (Laux et al., 2023) and flood risk studies (Hu et al., 2018).

ER levels become particularly interesting when the design and return periods are large and of similar magnitude. In such cases, the right-hand side of (5) can be approximated (see Appendix C) as:

$$(1-p)^{t_2-t_1+1} \sim \exp(-(t_2-t_1+1)/T),$$

where $T = 1/p$. For example, a return period 1.44 times longer than the design period corresponds to a 50% exceedance probability. Similarly, a ratio of 19.5 corresponds to a 5% exceedance probability, illustrating how proportional scaling of return and design periods maintains constant reliability.

To visualize how $z_p^{ER}$ compares to $z_p(t)$ with time, we assume in this paragraph that the random variable $Z_t$ follows a generalized

<figure><img src="image_5.png" /><figcaption>Fig. 1. Comparison of time-dependent return levels $z_p(t)$ and equivalent reliability (ER) levels under stationary and non-stationary assumptions. The left panel represents the stationary case with constant location parameter $\mu(t) = 0.7$, while the right panel represents a non-stationary case with a time-varying location parameter $\mu(t) = 0.7(t-2050)$. In both cases, the time-dependent return level corresponding to exceedance probability $p = 10^{-3}$ is shown in blue and the ER level $z_p^{ER}$ in orange. The shaded region represents the 2050–2100 design period. The figure highlights how ER levels provide a single summary for the design period in a non-stationary framework.</figcaption></figure>

extreme value (GEV) distribution. This family of distributions approximates the limiting behavior of the maximum of a series of i.i.d. observations (Coles, 2013; Katz et al., 2002). Its cdf is given by:

$$P(Z \le z | \theta) = \begin{cases} \exp[-(1 + \xi(\frac{z-\mu}{\sigma}))^{-1/\xi}] & \text{if } \xi \neq 0, \\ \exp[-\exp(\frac{z-\mu}{\sigma})] & \text{if } \xi = 0, \end{cases} \quad (6)$$

where $\theta = (\mu, \sigma, \xi)$ are the location, scale, and shape parameters, respectively.

Assuming $\sigma = 1.7$, $\xi = -0.2$ and a time-varying location $\mu(t) = 0.7(t-2050)$ for $t = 2000, \dots, 2100$, we simulate the effect of non-stationarity on return levels. These values are chosen to mimic high temperatures, see Section 3. We select 2050–2100 as the building design period and $p = 10^{-3}$ as the risk probability. Various illustrations with differing types of non-stationarity and inclusion in the GEV parameters are available in Appendix D. Since reliability is representative of a maximum over the design period, the corresponding level is mostly driven by the highest years of the period, even when the highest years are not at the end of the period.

The left panel of Fig. 1 displays the stationary case, where both the ER level $z_p^{ER}$ and the time-dependent return levels $z_p(t)$ defined by (2) remain constant at 6.3. In contrast, the right panel highlights that under non-stationarity, $z_p(t)$ increases from 6 to 10 over the design period, while $z_p^{ER}$ provides a time-aggregated summary value of 8.7.

This shows the difficulty of choosing between $z_p(2050) = 6$ or $z_p(2100) = 10$ for design, while $z_p^{ER}$ offers a clear interpretation as the level exceeded with 5% probability over the entire 51-year design period. This interpretation is straightforward: given the probability level $(1-p)^{t_2-t_1+1}$, it does not differ from a classical risk analysis with the advantage of one single summary for the design period in a non-stationary framework.

## 2.2. Accounting for estimation uncertainty: Predictive distribution

Quantifying uncertainty is crucial for informed decision-making, particularly when decisions rely on estimated values. Davison and Huser (2015) identifies several sources of uncertainty affecting the statistics of extreme events, among which two types are of particular relevance in our context. First, stochastic uncertainty arises from the inherent variability in future values, even when the underlying probability model is perfectly known. It corresponds to random variations

<footnote n="2">² Since $P[\max(Z_{t_1}, Z_{t_1+1}, \dots, Z_{t_2}) \le z_p] = F_{t_1}(z_p) \times \dots \times F_{t_2}(z_p) = F(z_p)^{t_2-t_1+1}$.</footnote>
of a draw within a given probability distribution. Secondly, estimation uncertainty arises from the inherent discrepancy between a fitted model and the true underlying model, even within the correct model family. It relates to the uncertainty in parameter estimates. Estimation uncertainty can be substantial, especially when models are complex or when the data available for model fitting is limited, as is typical of environmental data, such as temperature extremes. Failing to account for estimation uncertainty, or underestimating it, can foster overconfidence in the resulting estimates and confidence ranges and may lead to a dangerous underestimation of risk. For example, in June 2021, temperatures in Canada exceeded the upper bound derived from best estimates of GEV parameters inferred from historical data (Bartusek et al., 2022; Philip et al., 2022).

In structure design, standard practice involves evaluating high-order quantiles such as return levels corresponding to specified return periods, derived using frequentist maximum likelihood or Bayesian posterior estimates. These estimates are typically accompanied by confidence or credibility intervals (Davison and Huser, 2015) (also called a credible interval (Eberly and Casella, 2003)) to describe the associated estimation uncertainty. Despite underlying uncertainty, regulators and risk managers frequently demand a single, actionable value for design purposes (Mateus and Coonan, 2023; Rust et al., 2011). Limited guidance exists on how to incorporate estimation uncertainty to derive a single value for structural design. A common approach is to consider the upper bound of a confidence or credibility interval at a chosen confidence level (Bousquet and Bernardara, 2021; Rust et al., 2011). However, this approach poses several challenges. Confidence interval may be excessively wide, resulting in bounds that exceed physical plausibility (Fawcett and Green, 2018). On the contrary, short observational records with no or few low-likelihood events can lead to sampling bias and overly narrow intervals (Zeder et al., 2023), failing to adequately encompass future values.

In practice, the upper limit of a return level's confidence interval is influenced by two distinct user-specified inputs: the return period (or annual exceedance probability) and the confidence level of the confidence interval. The selection of these two criteria is often arbitrary, resulting in a wide range of levels being employed in practice. In applications, confidence levels range from 70% (Hamdi et al., 2018) to 95% (Yan et al., 2017; van den Brink and Können, 2011) or 99% (Toimil et al., 2021). Return periods also vary across applications, from 100-year rainfall (Gentilucci et al., 2023) to up to 50,000-year events (van den Brink and Können, 2011). These diverse combinations often lead to overlapping estimates, potentially confusing the true risk level under consideration. This ambiguity is illustrated in Fig. 2. In this example, the upper bound of the 99% credibility interval for a return level with $p = 10^{-2}$ slightly exceeds the upper bound of the 70% interval for $p = 2 \times 10^{-3}$. Given that these two quantities almost coincide, if this value is used for structural design, which risk level is truly being considered: is it $10^{-2}$ or $2 \times 10^{-3}$? Moreover, Coles and Tawn (2005) raised concerns regarding the exclusive use of confidence intervals of return levels for decision making, noting that “Framing return level questions in predictive terms is the most natural for decision making, while alternative approaches based on confidence intervals of classical return level estimates can be misleading due to their lack of a formal probabilistic interpretation”. Consequently, using their upper bound as a design point is debatable, as it does not provide a clear exceedance probability and may create confusion about the actual level of risk.

An alternative approach, well known in Bayesian statistics, relies on the predictive distribution (Coles and Tawn, 2005). This idea has a frequentist analog using likelihood (DAVISON, 1986). In a Bayesian framework, model parameters are treated as random variables whose probability distribution function (pdf) reflects prior parameter uncertainty. The predictive distribution is the marginal distribution of the outcome variable, averaged over the distributions of the model parameters:

$$P(Z \le z) = \int_{\theta} F(z|\theta)p(\theta)d\theta, \qquad (7)$$

<figure><img src="image_6.png" /><figcaption>Fig. 2. Comparison between predictive, median, and upper-bound estimates of return levels for exceedance probabilities p. The figure represents the value of various estimates (y-axis) for various p (x-axis): median estimate $z_p^{0.5}$ (blue), credibility intervals of levels 70%, 95% and 99% (decreasing shades of blue), and the predictive estimate $z_p^{\text{pred}}$ (orange). A gray line represents the value 7, to highlight an example where the upper bound of the 99% credibility interval for a return level with $p = 10^{-2}$ slightly exceeds the upper bound of the 70% interval for $p = 2 \times 10^{-3}$. The figure highlights how the difference between $z_p^{\text{pred}}$ and $z_p^{0.5}$ increases with p while credibility intervals potentially overlap.</figcaption></figure>

where $F(z|\theta)$ is the cdf of Z given a vector of parameters $\theta$ following a prior pdf $p(\theta)$. (7) represents where a new draw of z may lie, not conditioned on a specific value of $\theta$, but rather considering the current knowledge and uncertainty surrounding $\theta$. This directly integrates estimation uncertainty into the distribution of future outcomes. Similarly to (1), the predictive return level $z_p^{\text{pred}}$ for a return period $T = 1/p$ is then the $1-p$ quantile of the predictive cdf.

Prior distributions aim at encapsulating expert knowledge that will be updated using the sample $z^{\text{obs}}$ to create the posterior cdf. When necessary, the predictive distribution can be used with the posterior cdf to define and approximate a posterior predictive distribution.

Given a set of observations $z^{\text{obs}}$, the posterior predictive cdf is given by:

$$F(z|z^{\text{obs}}) = \int_{\theta} F(z|\theta)p(\theta|z^{\text{obs}})d\theta. \qquad (8)$$

A parametric form for $p(\theta|z^{\text{obs}})$ is rarely explicit, and methods like the Markov chain Monte Carlo (MCMC) are typically used to approximate it through sampling $\theta_{i \in \{1, N\}}$:

$$F(z|z^{\text{obs}}) \approx \frac{1}{N} \sum_{i=1}^{N} F(z|\theta_i). \qquad (9)$$

This produces a numerical estimate of the predictive distribution. Quantiles of this predictive distribution provide point summary estimates that incorporate parameter uncertainty.

Fawcett and Green (2018) and Coles and Tawn (2005) advocated for the use of predictive return levels as “the most natural for decision making” due to their intuitive interpretation as the expected level of a future event. Predictive distributions have been applied within a stationary framework to various environmental variables (Davison and Huser, 2015; Fawcett and Walshaw, 2006, 2008; Coles and Tawn, 1996a; Smith, 1999), including wind speed maxima (Fawcett and Green, 2018), sea level surges (Coles and Tawn, 2005; Fawcett and Walshaw, 2016), and solar flares (Griffiths et al., 2022). A key advantage of the predictive approach is its capacity to incorporate unfavorable, low-likelihood but plausible parameter values, making it particularly suited for evaluating low-likelihood, high-impact outcomes. For instance, Coles (2003) showed that a predictive return level could have anticipated the 1999 Venezuela record rainfall. Lastly and importantly, unlike interval-based methods, which require both a return period and
a confidence level, the predictive return level depends only on a single user-specified parameter: the acceptable risk level p.

To illustrate the computation of predictive return levels, we consider a GEV-distributed variable (Coles, 2013; Katz et al., 2002) with cdf defined in (6) and parameter priors given by:

$$ \theta = \begin{cases} \mu & \sim \mathcal{N}(0, 0.5), \\ \log(\sigma) & \sim \mathcal{N}(0.5, 0.05), \\ \xi & \sim \mathcal{N}(-0.2, 0.05), \end{cases} \qquad (10) $$

where $\mathcal{N}(a, b)$ denotes a Gaussian random variable with mean a and standard deviation b, and $\sigma$ has a strictly positive support. Various illustrations with differing distributions for $\mu$ and $\sigma$ are available in Appendix E.

From $N = 1000$ draws $\{\theta_1, \dots, \theta_N\}$, we compute return levels $z_p(\theta_i)$ such that $F(z_p(\theta_i)|\theta_i) = 1-p$. The median estimate $z_p^{0.5}$ and the credibility interval, defined by the $\alpha$ and $1-\alpha$ quantiles, are derived from this distribution. For comparison with the parameter-estimate-centered setting, we define $z_p(\theta^{\text{BE}})$ as the return level based on the Bayes estimator $\theta^{\text{BE}}$, taken as the multivariate median of $\{\theta_1, \dots, \theta_N\}$ (Vardi and Zhang, 2000). An approximated predictive distribution is obtained from (9) by sampling from the distributions GEV($\theta_1$),...,GEV($\theta_N$). $z_p^{\text{pred}}$ is defined as the $1-p$ quantile from this distribution.

Fig. 3 illustrates how $z_p^{\text{pred}}$ encompasses a broader range of outcomes, particularly in the tail, than $z_p(\theta^{\text{BE}})$ or even $z_p^{0.5}$. p is set at $10^{-5}$. The first panel illustrates the general interval approach. It displays a large number of GEV pdf from GEV($\theta_1$),...,GEV($\theta_N$), the median estimate $z_p^{0.5}$ and the credibility interval ($z_p^{0.05}$ and $z_p^{0.95}$). The second panel illustrates the difference between the pdf based on the Bayes estimator $\theta^{\text{BE}}$ (blue) and the approximated predictive distribution (orange). The third panel displays the cdf of the previous distributions, and the last panel is zoomed in on the tail. While all the $\xi_i$ sampled here are negative, the tail of the predictive distribution is thicker than any of the individual GEV distributions. This reflects the predictive distribution’s inclusion of rare but plausible outcomes not captured by a single parameter set.

Fig. 2 illustrates how the difference between $z_p^{\text{pred}}$ and $z_p^{0.5}$ increases with p while credibility intervals potentially overlap. For various values of p, it displays the median estimate $z_p^{0.5}$, credibility intervals of levels 70%, 95% and 99% and the predictive estimate $z_p^{\text{pred}}$. The value 7 is an example of overlap where the upper bound of the 99% credibility interval for a return level with $p = 10^{-2}$ slightly exceeds the upper bound of the 70% interval for $p = 2 \times 10^{-3}$. It illustrates how the choice of the confidence level can influence the result, and partly offset or confuse the choice of the return level. For higher exceedance probabilities such as $p = 0.02$, the predictive and median-based quantiles align closely. However, for rare events such as $p = 10^{-5}$, $z_p^{\text{pred}}$ may even exceed the upper bound of a high-level credibility interval. This divergence highlights the increased importance of accounting for parameter uncertainty in extreme quantiles.

Consistent with the recommendations of Davison and Huser (2015), Coles and Tawn (1996b), we advocate using predictive quantiles as point estimates for structural design. Predictive quantiles provide a compelling alternative to conventional interval-based summaries to account for both stochastic and estimation uncertainty.

## 2.3. Combining both concepts: Predictive equivalent reliability

We now propose a unified framework that combines a non-stationary framework and parameter uncertainty through the concept of Predictive Equivalent Reliability (PER). This framework builds on the predictive distribution (Section 2.2) and the equivalent reliability (ER) framework (Section 2.1), and is designed to assess the probability of exceedance over a design period, while fully accounting for uncertainty in model parameters.

Let $Z_t$ be a random variable of interest at year t, characterized by a vector of parameters $\theta$ with distribution $p(\theta)$. We consider a design

<figure><img src="image_7.png" /><figcaption>Fig. 3. Comparison between the predictive estimate, the median estimate, and the parameter-estimate-centered estimate of the 10$^5$-year return level. Let $\{\theta_1, \dots, \theta_N\}$ be a sample of parameter values drawn from the parameter distribution $p(\theta)$. p is set at $10^{-5}$. Panel a represents various probability distribution functions (pdfs) for GEV($\theta_1$),...,GEV($\theta_N$). Vertical bars show the median estimate of the return level $z_p^{0.5}$ (blue), and the corresponding 5%–95% credibility interval ranging from $z_p^{0.05}$ to $z_p^{0.95}$ (black). Panel b represents the pdf given by the Bayes estimator $\theta^{\text{BE}}$ (blue), and the approximated predictive distribution from (9) (orange). Panel c represents the cumulative distribution functions (CDFs) corresponding to the previous distributions shown in Panel b. Panel d is a zoom on the right tail of the CDFs from Panel c. For panels b, c, and d, vertical bars show return levels calculated using either a predictive approach (orange) or the Bayes estimator. This figure highlights how the predictive distribution covers rare but plausible outcomes not captured by a single parameter set.</figcaption></figure>

period from year $t_1$ : $t_2$. The predictive reliability associated with a threshold z over this period is defined as:

$$ \mathrm{PR}_{t_1:t_2}(z|z^{\mathrm{obs}}) = P[\max(Z_{t_1}, Z_{t_1+1}, \dots, Z_{t_2}) \le z|z^{\mathrm{obs}}], \qquad (11) $$

$$ = \int_{\theta} P[\max(Z_{t_1}, Z_{t_1+1}, \dots, Z_{t_2}) \le z|\theta] p(\theta|z^{\mathrm{obs}}) d\theta. \qquad (12) $$

The predictive equivalent reliability level (PER) $z_p$ is defined as the solution to:

$$ \mathrm{PR}_{t_1:t_2}(z_p) = (1-p)^{t_2-t_1+1} \qquad (13) $$
Given a set of observations $z^{\text{obs}}$, the posterior predictive reliability can be approximated and estimated numerically by:

$$PR_{t_1:t_2}(z | z^{\text{obs}}) \approx \frac{1}{N} \sum_{i=1}^{N} P[\max(Z_{t_1}, Z_{t_1+1}, \dots, Z_{t_2}) \le z | \theta_i], \quad (14)$$

where $\theta_i$ are samples from $p(\theta | z^{\text{obs}})$. An approximation of the PER of annual probability $p$ is given by the quantile $(1-p)^{t_2-t_1+1}$ of the approximated predictive reliability.

As a framework combining multiple concepts, the PER level is particularly suited to situations where a single, interpretable value is needed to summarize risk over a multi-year period under uncertainty. Examples of the use of a single index appear in a large range of applications. In structural design (Mateus and Coonan, 2023; Rust et al., 2011), only one value may be kept as a reference. In a context of adaptation, the viability of different projects could be evaluated using a PER level. This concept is also of interest in decision-making. For instance, a decision maker may wish to invest in a forest replanting project only if the risk of failure is below 5% (corresponding to an annual risk of failure of 0.1%, if a 50-year project is considered). Only tree species able to survive over the corresponding PER level would be selected in such a case.

The PER level offers two major advantages over classical return levels. First, it provides a single value for the entire design period. Classical return levels were originally defined in a stationary context and are often used to assess the risk of excess either for a single year or over the associated return period, assuming identical yearly conditions. In contrast, PER levels build upon equivalent reliability levels as defined in Section 2.1. ER levels inherently assess exceedance risk over a period, even if the distribution of events varies. It also avoids the conceptual pitfall of using a single-year metric to characterize multi-year risk, a limitation previously noted in the literature Salas and Obeysekera (2014), Serinaldi (2015), Yan et al. (2017) and Liang et al. (2016). Second, it integrates multiple sources of uncertainty in a single value. The typical use case of return levels often includes a confidence interval to express parameter uncertainty, which requires a subjectively chosen confidence level. The upper bound is sometimes used as a conservative estimate (Bousquet and Bernardara, 2021; Rust et al., 2011), but its probabilistic interpretation is unclear. In contrast, the PER level relies on the predictive distribution as defined in Section 2.2. It encompasses similar sources of uncertainty as confidence intervals while retaining the original index's consistent probabilistic interpretation.

# 3. Application

## 3.1. Data set

Our approach is applied to the study of annual temperature maxima in the city of Pierrelatte (44.33°N 4.73°E), in the South of France. Local observational data are available from 1985 to 2021 and sourced from the Pierrelatte meteorological station. We focus on annual maxima of daily maximum temperatures (TXx) to characterize extreme heat events and denote the corresponding time series $z^{\text{obs}}$. Simulated TXx time series are derived from 26 General Circulation Models (GCMs) from the CMIP6 ensemble, covering the period from 1850 to 2100. We keep both historical simulations and future projections under the SSP1-1.9, SSP2-4.5, and SSP5-8.5 scenarios. A visual overview of these simulations is provided in Fig. A.6.

To model non-stationarity, we use the annual mean temperature over continental Europe (35–70°N, 10 W–30°E) as a covariate. Observational data for this covariate is taken from the CRUTEM5 (Osborn et al., 2021) dataset and spans 1850–2020, noted $x^{\text{obs}}$. Simulated covariate series from the previous GCMs are also used and shown in Fig. A.5.

## 3.2. Statistical model and (Bayesian) estimation

Following the framework of Robin and Ribes (2020), we model the annual maxima $Z_t$ using a non-stationary Generalized Extreme Value (GEV) distribution (Coles, 2013), with time-varying parameters driven by the covariate $X_t$.

The parameterization of the GEV distribution is given by:

$$\theta(t) = \begin{cases} \mu(t) = \mu_0 + \mu_1 X_t, \\ \sigma(t) = \exp(\sigma_0 + \sigma_1 X_t), \\ \xi(t) = \xi. \end{cases} \quad (15)$$

$\mu(t)$ and $\sigma(t)$ are the location and scale parameters, with the log-link ensuring positivity of $\sigma(t)$. The shape parameter $\xi$ is assumed constant. This specific parameterization has been selected by Robin and Ribes (2020) using a likelihood ratio test on simulated data from the CMIP5 dataset.

The vector of parameters $\theta = (\mu_0, \mu_1, \sigma_0, \sigma_1, \xi)$ is estimated using Robin and Ribes (2020)'s Bayesian approach. In the next part, we briefly review the key concepts of Robin and Ribes (2020) and refer to it for full details. First, non-stationary GEV parameters and covariate values are estimated for each GCM via maximum likelihood. The resulting sample of model estimates is used to fit a multivariate Gaussian distribution, which is taken as a prior for $\phi = (\theta, X_{t_1:t_2})$.

Posterior inference is done in two steps:

* The prior is updated with observed covariate data using a conjugate distribution (Ribes et al., 2020) to obtain $p(\phi|x^{\text{obs}})$.

* This intermediate posterior is then updated using local data $z^{\text{obs}}$ to obtain $p(\phi|x^{\text{obs}}, z^{\text{obs}})$. Since no conjugate is available in this case, sampling from the posterior distribution is done using the NUTS (No-U-Turn Sampler) MCMC algorithm (Homan and Gelman, 2014). 10000 posterior samples are obtained and presented in Appendix G.

Calculations are done using Robin and Ribes (2020)'s NSSEA and Pystan (Stan Development Team, 2024) Python packages.

This approach allows us to sample from the posterior distribution of the parameters, enabling estimation of quantities such as instantaneous return levels and predictive equivalent reliability (PER) levels.

## 3.3. Results

Results of this application are presented in Fig. 4, with additional summaries including credibility interval of the ER levels and predictive time-dependent return levels shown in Appendix B. Classical goodness-of-fit diagnostics to check our GEV hypothesis (Coles, 2013) are presented in Appendix F. Under each emission scenario, quantile-quantile plots and Kolmogorov–Smirnov tests do not contradict our GEV hypothesis, suggesting our non-stationary GEV model is appropriate.

For scenario SSP2-4.5 (yellow) and an annual probability of exceedance $p = 10^{-3}$, the PER level for 2050–2100 is 47.4 °C. As a comparison, the median time-dependent return level of annual probability $p = 10^{-3}$ increases from 44.7 °C to 46.5 °C over the same period. Similar but exacerbated behavior is visible for scenario SSP5-8.5 (red), where the PER level is 53.5 °C while median time-dependent return levels increase from 45.5 °C to 51.6 °C. For scenario SSP1-1.9 (blue), return levels increase up to 43.8 °C in 2062, and then stabilize. For each scenario, the PER level's total exceedance probability over 2050–2100 is 0.05, whereas the corresponding annual exceedance probabilities vary from year to year. This illustrates the benefit of using the PER level as a single value representing exceedance risk over the design period.

For both SSP2-4.5 and SSP5-8.5 scenarios, exhibiting continued warming over the 2050–2100 period, the PER level is primarily influenced by extremes occurring near the end of the design period. For scenario SSP1-1.9, the PER value is also driven by the warmest
<figure><img src="image_8.png" /><figcaption>Fig. 4. Comparison of time-dependent return levels and predictive equivalent reliability (PER) levels under three emission scenarios. This figure represents estimated time-dependent return levels $z_p(t)$ for $p = 10^{-3}$ (including credibility intervals at 5% level) and corresponding PER levels for three emission scenarios. SSP1-1.9, SSP2-4.5, and SSP5-8.5 are represented using respectively blue, yellow, and red. The shaded region represents the 2050–2100 design period. Observed annual maxima are shown with black dots. The figure highlights how PER levels provide one single summary value to represent exceedance risk over the design period for a given scenario.</figcaption></figure>

year, in the middle of the design period. In all cases, the PER value exceeds the highest median return level $z_p^{50\%}(t)$ (either 2100 or 2062), highlighting the importance of accounting for both parameter and stochastic uncertainty when estimating design levels.

## 4. Discussion

Designing structures sensitive to hot temperatures for the coming decades requires considering the effects of climate change. The observed non-stationarity in temperature extremes, driven by ongoing warming trends, is expected to intensify across all emissions scenarios (Intergovernmental Panel On Climate Change (Ipcc), 2023). Accordingly, any robust design strategy must reflect the evolving distribution of climate variables and the increased uncertainty introduced by limited observations and divergent projections from state-of-the-art general circulation models (GCMs). To address this, we extend the traditional concept of return levels using the equivalent reliability level, which remains valid in non-stationary settings. The predictive framework is adopted to encapsulate both estimation and stochastic uncertainties. The predictive equivalent reliability level (PER) integrates both concepts into a single summary index that quantifies the risk of exceedance over a design period, while accounting for temporal trends and uncertainty on the parameters. The approach was applied to annual maxima of temperatures in southern France.

In evaluating risk over a design period, one might instinctively consider the probability of exceedance in the worst year, such as the year with the strongest warming effect. Rootzén and colleagues (2013) (Rootzén and Katz, 2013) define the Minimax Design Life Level using the maximum annual exceedance probability over the design period, resulting in more conservative design levels. However, we concur with (Rootzén and Katz, 2013) that aggregated risk over the design period is both more informative and more consistent with the usual practice in a stationary context. Since Equivalent reliability (ER) levels align with this principle while maintaining consistency with the existing stationary return level framework, we regard them as the most appropriate index.

This study specifically considers the risk of exceeding a specified threshold. Depending on the application, other definitions of risk, such as the number of exceedances within a given period or their

duration, might be more appropriate. The predictive approach could be integrated into other methodologies, such as the Bayesian extension of intensity-duration-frequency (IDF) curves proposed by Van de Vyver (2015).

By design, the performance of the posterior predictive distribution also depends on prior choice. As outlined in van de Schoot et al. (2021), priors can be informed using various methodologies such as empirical data, literature meta-analyses, or expert elicitation. This application used an unbounded multivariate Gaussian prior; however, alternative prior specifications, e.g. with bounded tails, could better reflect physical constraints or domain knowledge.

Concerning the topic of verification, it is well known that the quantile score provides a clear link between some classical proper scoring rules and return levels (see, e.g., Eq. (4) and (8) in Pic et al., 2025). Proposing a scoring rule suitable for PER estimation would be an interesting research avenue if this metric had to be widely used in decision-making. This would also help to go beyond visual heuristics when comparing PER and return levels.

A notable limitation of the current method is the treatment of the scenario uncertainty. Emission scenarios are treated individually, omitting the scenario uncertainty, as a major contributor to long-term projection uncertainty (Hu et al., 2023). While PER levels are designed to combine as much information as possible and to encompass less likely scenarios, thereby accommodating uncertainty in the climate system's response, its current formulation omits the uncertainty associated with different emission trajectories. Incorporating a probabilistic description of the scenario uncertainty would be attractive to provide a more comprehensive risk picture. However, assigning probabilities to emission scenarios remains a challenge, and is beyond the scope of our paper, although previous works (Hausfather and Peters, 2020) have highlighted that all scenarios are not equally likely. The issue of scenario selection is particularly conspicuous in decision-making processes, where scenario selection often reflects both scientific judgment and normative choices related to policy or planning goals (Grant et al., 2020).

Additionally, because the predictive distribution samples from both parameter and process uncertainty, its upper tail may include values that are statistically plausible but seem physically unrealistic. While this typically occurs at extremely low probability levels, it can still pose interpretative challenges. To prevent the model from generating physically impossible values, one can use physical sources of information to specify a statistical upper bound (Noyelle et al., 2023). Further theoretical developments regarding the requirements for prior distributions to yield a bounded predictive distribution, along with their applicability to real variables, represent potential areas for future research.

For illustration, this method has been applied to a single point of interest and three separate future emission scenarios. Broader applications across regions and scenarios may be of practical interest for many users interested in designing life levels. It may also reveal new insights about future extreme temperatures in the context of a rapidly shifting climate.

## CRediT authorship contribution statement

**Occitane Barbaux:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis. **Philippe Naveau:** Writing – review & editing, Writing – original draft, Supervision, Funding acquisition, Formal analysis, Conceptualization. **Nathalie Bertrand:** Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization. **Aurélien Ribes:** Writing – review & editing, Writing – original draft, Visualization, Supervision, Methodology, Funding acquisition, Formal analysis, Conceptualization.
## Code availability

Code and data used in this work are available here: https://github.com/Occitane-Barbaux/Integrating-non-stationarity-and-uncertainty-in-design-life-levels

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The authors thank the three anonymous referees for their helpful comments that greatly improved the quality of the manuscript. This study has received funding from Autorité de sûreté nucléaire et de radioprotection (ASNR), France, Météo-France, France, Agence Nationale de la Recherche, France - France 2030 as part of the PEPR TRACCS programme under grand number ANR-22-EXTR-0005, the PEPR IRIMONT, France (France 2030 ANR-22-EXIR-0003), SHARE PEPR Maths-Vives, France (ANR-24-EXMA-0008), EXSTA, France, and from the Geolearning research chair, France, a joint initiative of Mines Paris and the French National Institute for Agricultural Research (INRAE).

## Appendix A. Data

Data from 26 GCMs was extracted and used in the final application. For some models, several runs are available, so several values are available for each year. All data is shown as anomalies with the reference period 1986–2015.

## Appendix B. Partial methodology

The calculation of a median ER Level and corresponding credibility interval is possible, as illustrated in Fig. B.7

Similarly, the calculation of a predictive effective return level is possible, as illustrated in Fig. B.8

<figure><img src="image_9.png" /><figcaption>Fig. A.5. Dots in black are the real observations. Each color is one of the GCM available, and each dot is a value of covariate X, the annual mean temperature of the daily mean temperature over Continental Europe.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. A.6. Dots in black are the real observations. Each color is one of the GCM available, and each dot is a value of the annual maxima of the daily maxima of temperature over Pierrelate.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. B.7. The median of the annual quantile $z_p^{50\%}(t)$ (given by (2)) for $p = 10^{-3}$ is in red and for $p = 0.5$ in gray. The median ER level over the 2050–2100 period $z_p^{ER}$ (given by (5)) is shown in green. Confidence intervals are all shown in similar colors at the 5% level. Existing observations of local annual maxima are shown in black.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. B.8. The median of the annual quantile $z_p^{50\%}(t)$ (given by (2)) for $p = 10^{-3}$ is in red and for $p = 0.5$ in gray. The posterior predictive quantile level $z_p^{pred}(t)$ (given by (9)) is shown in green. Confidence intervals are all shown in similar colors at the 5% level. Existing observations of local annual maxima are shown in black.</figcaption></figure>
## Appendix C. Design period and return period relationship

We define $T = \frac{1}{p}$ (the return period) and $T_{1:2} = t_2 - t_1 + 1$ (the design period). We suppose a relation between these, where $T = c \times T_{1:2}$ with $c \in R^{*+}$.

The corresponding reliability can be expressed as:

$$ (1 - \frac{1}{T})^{T_{1:2}} = (1 - \frac{1}{cT_{1:2}})^{T_{1:2}} $$

Using the following property: for $k, n \in R^{*+}$,

$$ \lim_{n \to \infty} (1 - \frac{k}{n})^n = e^{-k} $$

$e^{-k}$ is actually an upper bound: $(1 - \frac{k}{n})^n \le e^{-k} \forall n$.

We have:

$$ \lim_{T_{1:2} \to \infty} (1 - \frac{1}{cT_{1:2}})^{T_{1:2}} = e^{-\frac{1}{c}} $$

Where $e^{-\frac{1}{c}}$ is an upper bound on the probability of no event exceeding the threshold (i.e. $\forall T_{1:2}$, the probability of excess is equal or above $1 - e^{-\frac{1}{c}}$)

Examples:

Defining $p_{tot}$ as the target reliability, i.e. the probability of no observation above the threshold during the design period, we have $c = -\frac{1}{\ln(p_{tot})}$.

For $p_{tot} = 0.5$, $c \approx 1.44$.

For $p_{tot} = 0.7$, $c \approx 2.8$.

For $p_{tot} = 0.95$, $c \approx 19.5$.

For $p_{tot} = 0.99$, $c \approx 99.5$.

## Appendix D. Sensitivity of equivalent reliability levels to various non-stationary scenarios

The form of the non-stationary scenario and how non-stationarity is incorporated into the model parameters directly influence Equivalent Reliability (ER) levels and consequently PER levels.

Here, we explore several types of non-stationarity to illustrate the sensitivity of ER through a numerical study. In practice, we reuse the toy example presented in Section 2.1:

We assume that the random variable $Z_t$ follows a generalized extreme value (GEV) distribution, with a cdf:

$$ P(Z \le z | \theta) = \begin{cases} \exp[-(1 + \xi(\frac{z-\mu}{\sigma}))^{-1/\xi}] & \text{if } \xi \neq 0, \\ \exp[-\exp(\frac{z-\mu}{\sigma})] & \text{if } \xi = 0, \end{cases} \quad (D.1) $$

with $\theta = (\mu, \sigma, \xi)$ the location, scale and shape parameters.

We choose a constant shape parameter (set to -0.2, a realistic value for extreme temperatures), while location and scale are defined as varying with time with the following parametrization:

$$ \left\{ \begin{array}{l} \mu(t) = 0.7 X(t), \\ \log(\sigma(t)) = 0.5 + \sigma_1 X(t), \\ \xi(t) = -0.2. \end{array} \right. \qquad (D.2) $$

Depending on the example, $\sigma_1$ is either worth 0 (stationary scale) or 0.05 (time-dependent scale). The covariate $X(t)$, which carries the non-stationarity, is defined for four possible scenarios:

1. A linear scenario, as presented in Section 2.1, where $X(t) = t - 2050$, for $t = 2000, ..., 2100$

2. An exponential scenario, where $X(t) = 100(1 - \exp[-0.05(t - 2050)]) - 50$, for $t = 2000, ..., 2100$

3. An hyperbolic scenario, where $X(t) = 100 * (\tanh[0.015(t - 2050)]) - 50$, for $t = 2000, ..., 2100$

4. A gaussian, where $X(t) = 100 * \exp[-10 * (\frac{t-2050}{50})^2] - 50$, for $t = 2000, ..., 2100$

We select 2000-2100 as the building design period, and the risk probability is set to $p = 10^{-3}$.

The first row of Fig. D.9 displays applications where the scale is stationary with $\sigma_1 = 0$ while the second displays applications where $\sigma_1 = 0.05$. Similarly, the first column on the left displays applications where the covariate follows a linear scenario, the second applications where it follows an exponential scenario, the third applications where it follows a hyperbolic scenario, and the fourth applications where it follows a Gaussian scenario.

The results differ between each application. There is no visible nor direct relationship between the evolution of time-dependent return levels and the value of the ER level, but the ER level is systematically quite close to, and slightly lower than, the highest return level over the period. Noticeably, ER levels are driven by the highest years of the period, particularly if very rare events are being considered. This is observed regardless of the “shape” of the non-stationarity, including when the highest years are not at the end of the period. Similarly, there is no direct relationship with the inclusion of stationarity in the scale parameter, but all level values are systematically higher.

## Appendix E. Sensitivity of predictive levels to parameter distributions

Uncertainty is included in PER levels through the use of the predictive distribution. Several factors (such as reduced sample size, increased spread of the covariate, inference issues, etc.) could lead to an increased spread of the posterior distribution of the parameters and consequently to a predictive distribution. It is useful to highlight the sensitivity of the predictive levels, and thus the PER levels, to the distribution of the parameters. To visualize the effect of various parameter distributions, we reuse the toy example presented in Section 2.2:

We assume that the random variable of interest follows a GEV distribution with its cdf defined in (6). We want to compare the effect of different standard deviations in the distributions of the location and scale parameters. We compare the results of four applications, with the GEV location, scale, and shape parameters following the probability distributions below:

$$ \theta = \begin{cases} \mu & \sim \mathcal{N}(0, 0.5) \quad \text{or} \quad \mathcal{N}(0, 1), \\ \log(\sigma) & \sim \mathcal{N}(0.5, 0.05) \quad \text{or} \quad \mathcal{N}(0.5, 0.2), \\ \xi & \sim \mathcal{N}(-0.2, 0.05), \end{cases} $$

where $\mathcal{N}(a, b)$ denotes a Gaussian random variable with mean $a$ and standard deviation $b$, and $\sigma$ has a strictly positive support. Even though the spread of the shape parameter is not explored in this illustration, the corresponding distribution has a mean of -0.2 (a realistic value for extreme temperatures, which are usually considered bounded.)

Fig. E.10 displays the distribution of all the $z_p(\theta_l)$, the median $z_p^{0.5}$ and quantiles $z_p^a$ and $z_p^{1-a}$ of the samples of $z_p(\theta_l)$, with $p = 10^{-4}$ and $\alpha = 0.05$, corresponding to a credibility interval at the 90% level. The predictive quantile $z_p^{\text{pred}}$ is also displayed for comparison.

Both the values and the proximity between each index differ between applications.

## Appendix F. Goodness-of-fit

In this study, we rely on Robin and Ribes (2020) for the statistical model selection. We also rely on classical extreme value theory distributional diagnostic tools (Coles, 2013) to assess goodness-of-fit on the residuals defined as follows. Given the estimated parameters set $\hat{\mu}_t, \hat{\sigma}_t, \hat{\xi}_t$, the standardized variables

$$ \hat{Z}_t = \frac{Z_t - \hat{\mu}_t}{\hat{\sigma}_t} $$

should approximately follow a GEV(0, 1, $\hat{\xi}_t$) stationary distribution. Using this transformation, we perform a bootstrapped Kolmogorov-Smirnov test for the median estimate of the posterior distribution of
<figure><img src="image_13.png" /><figcaption>**Fig. D.9.** For eight different applications, the annual quantile of probability $p = 10^{-3}$ (given by (2)) is in red, and the Equivalent Reliability level (given by (5)) is shown in dashed green. The blue shade represents the 2000–2100 design period. The year in which the annual quantile first exceeds the Equivalent Reliability level is given in black. The first row displays applications where the scale is stationary with $\sigma_1 = 0$, while the second displays applications where $\sigma_1 = 0.05$. Similarly, the first column on the left displays applications where the covariate follows a linear scenario, the second applications where it follows an exponential scenario, the third applications where it follows a hyperbolic scenario, and the fourth applications where it follows a Gaussian scenario.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. E.10.** For four different applications, the distribution of all the $z_p(\theta_t)$ is in blue. Vertical bars show the estimates of the return level $z_p$, for $p = 10^{-4}$; the median $z_p^{0.5}$ (blue), and the 5%–95% quantiles $z_p^{0.05}$, $z_p^{0.95}$ (black). $z_p^{\text{pred}}$, the quantile of the predictive distribution, is shown in pink. The first row displays applications where $\log(\sigma) \sim \mathcal{N}(0.5, 0.05)$ while the second displays applications where $\log(\sigma) \sim \mathcal{N}(0.5, 0.2)$. Similarly, the first column on the left displays applications where $\mu \sim \mathcal{N}(0, 0.5)$, while the second column displays applications where $\mu \sim \mathcal{N}(0, 1)$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. F.11. For three different emission scenarios, the figure shows sample quantile-quantile plots of the standardized observation using the estimated parameters.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. G.12. For the five parameters of the non-stationary GEV, the figure shows the prior distribution (blue), the sampled posterior distribution (red), and the maximum likelihood estimate based only on the available observations.</figcaption></figure>

the parameters. The null distribution of the KS test statistic has been sampled to account for the use of fitted parameters (Zeimbekakis et al., 2024). Each scenario is assessed individually and provides the following p-values:

• For scenario SSP 5–8.5, the *p*-value is 0.82.

• For scenario SSP 2–4.5, the *p*-value is 0.86.

• For scenario SSP 1–1.9, the *p*-value is 0.35.

Hence, with a *p*-value well above 0.05, each scenario is consistent with the null hypothesis ($\hat{Z}_i$ are GEV distributed). As an additional validation criterion, Fig. F.11 below displays a quantile-quantile plot for each scenario, comparing the standardized quantiles of the observed sample (y-axis) versus the fitted posterior distribution (x-axis).

Visually, points generally align along the diagonal, meaning the fit is satisfactory. Overall, our GEV distributional assumption does not seem to be rejected.

## Appendix G. MCMC sampling

The MCMC procedure is done using the STAN implementation of the NUTS (No-U-Turn Sampler) MCMC algorithm. In this implementation, hyperparameters are selected automatically. The obtained samples are then used for both the PER level and the return levels (see Fig. G.12).

## Data availability

I have shared a link to the code and the data in the tex file https://github.com/Occitane-Barbaux/Integrating-non-stationarity-and-uncretainty-in-design-life-levels.

## References

Asseng, S., Spänkuch, D., Hernandez-Ochoa, I.M., Laporta, J., 2021. The upper temperature thresholds of life. Lancet Planet Heal. 5 (6), e378–e385. http://dx.doi.org/10.1016/S2542-5196(21)00079-6.

Bartusek, S., Kornhuber, K., Ting, M., 2022. 2021 north American heatwave amplified by climate change-driven nonlinear interactions. Nat. Clim. Chang. 12 (12), 1143–1150. http://dx.doi.org/10.1038/s41558-022-01520-4, Publisher: Nature Publishing Group, URL https://www.nature.com/articles/s41558-022-01520-4.

Bousquet, N., Bernardara, P. (Eds.), 2021. Extreme Value Theory with Applications to Natural Hazards: From Statistical Theory to Industrial Practice. Springer International Publishing, Cham, http://dx.doi.org/10.1007/978-3-030-74942-2, URL https://link.springer.com/10.1007/978-3-030-74942-2.

van den Brink, H.W., Können, G.P., 2011. Estimating 10000-year return values from short time series. Int. J. Climatol. 31 (1), 115–126. http://dx.doi.org/10.1002/joc.2047, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/joc.2047.

Coles, B.S., 2003. The use and misuse of extreme value models in practice. In: Extreme Values in Finance, Telecommunications, and the Environment. Chapman and Hall/CRC, pp. 90–111, Num Pages 22.

Coles, S., 2013. An Introduction to Statistical Modeling of Extreme Values. Springer Science & Business Media, Google-Books-ID, SonbBwAAQBAJ.
Coles, S.G., Tawn, J.A., 1996a. A Bayesian analysis of extreme rainfall data. J. R. Stat. Soc. Ser. C. Appl. Stat. 45 (4), 463-478. http://dx.doi.org/10.2307/2986068.

Coles, S.G., Tawn, J.A., 1996b. Modelling extremes of the aral rainfall process. J. R. Stat. Soc. Ser. B Stat. Methodol. 58 (2), 329-347. http://dx.doi.org/10.1111/j.2517-6161.1996.tb02085.x, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2517-6161.1996.tb02085.x.

Coles, S., Tawn, J., 2005. Bayesian modelling of extreme surges on the UK east coast. Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 363 (1831), 1387-1406. http://dx.doi.org/10.1098/rsta.2005.1574, Publisher: Royal Society, URL https://royalsocietypublishing.org/doi/full/10.1098/rsta.2005.1574.

Cooley, D., 2013. Return periods and return levels under climate change. In: Aghaouchak, A., Easterling, D., Hsu, K., Schubert, S., Sorooshian, S. (Eds.), Extremes in a Changing Climate: Detection, Analysis and Uncertainty. In: Water Science and Technology Library, Springer Netherlands, Dordrecht, pp. 97-114. http://dx.doi.org/10.1007/978-94-007-4479-0-4.

DAVISON, A.C., 1986. Approximate predictive likelihood. Biometrika 73 (2), 323-332. http://dx.doi.org/10.1093/biomet/73.2.323.

Davison, A.C., Huser, R., 2015. Statistics of extremes. Annu. Rev. Stat. Appl. 2 (Volume 2, 2015), 203-235. http://dx.doi.org/10.1146/annurev-statistics-010814-020133, Publisher: Annual Reviews, URL https://www.annualreviews.org/content/journals/10.1146/annurev-statistics-010814-020133.

van de Schoot, R., Depaoli, S., King, R., Kramer, B., Märtens, K., Tadesse, M.G., Vannucci, M., Gelman, A., Veen, D., Willemsen, J., Yau, C., 2021. Bayesian statistics and modelling. Nat Rev Methods Prim. 1 (1), 1-26. http://dx.doi.org/10.1038/s43586-020-00001-2, Publisher: Nature Publishing Group, URL https://www.nature.com/articles/s43586-020-00001-2.

Eberly, L.E., Casella, G., 2003. Estimating Bayesian credible intervals. J. Statist. Plann. Inference 112 (1), 115-132. http://dx.doi.org/10.1016/S0378-3758(02)00327-0, Special issue II: Model Selection, Model Diagnostics, Empirical Bayes and Hierarchical Bayes, URL https://www.sciencedirect.com/science/article/pii/S0378375802003270.

Fawcett, L., Green, A.C., 2018. Bayesian posterior predictive return levels for environmental extremes. Stoch Env. Res Risk Assess 32 (8), 2233-2252. http://dx.doi.org/10.1007/s00477-018-1561-x.

Fawcett, L., Walshaw, D., 2006. A hierarchical model for extreme wind speeds. J. R. Stat. Soc. Ser. C. Appl. Stat. 55 (5), 631-646. http://dx.doi.org/10.1111/j.1467-9876.2006.00557.x.

Fawcett, L., Walshaw, 2008. Modelling environmental extremes (part 5). URL http://www.mas.ncl.ac.uk/~nlf8/shortcourse/part5.pdf.

Fawcett, L., Walshaw, D., 2016. Sea-surge and wind speed extremes: optimal estimation strategies for planners and engineers. Stoch. Environ. Res. Risk Assess. http://dx.doi.org/10.1007/s00477-015-1132-3, Publisher: Newcastle University, URL https://eprints.ncl.ac.uk.

Fischer, E.M., Sippel, S., Knutti, R., 2021. Increasing probability of record-shattering climate extremes. Nat. Clim. Chang. 11 (8), 689-695. http://dx.doi.org/10.1038/s41558-021-01092-9, Publisher: Nature Publishing Group, URL https://www.nature.com/articles/s41558-021-01092-9.

Gentilucci, M., Rossi, A., Pelagagge, N., Aringoli, D., Barbieri, M., Pambianchi, G., 2023. GEV analysis of extreme rainfall: Comparing different time intervals to analyse model response in terms of return levels in the study area of central Italy. Sustainability 15 (15), 11656. http://dx.doi.org/10.3390/su151511656, Number: 15 Publisher: Multidisciplinary Digital Publishing Institute, URL https://www.mdpi.com/2071-1050/15/15/11656.

Gonzalez, P., Naveau, P., Thao, S., Worms, J., 2025. A statistical method to model non-stationarity in precipitation records changes. Geophysical Research Letters 52 (6), http://dx.doi.org/10.1029/2023GL017021, e2023GL017021.

Grant, N., Hawkes, A., Napp, T., Gambhir, A., 2020. The appropriate use of reference scenarios in mitigation analysis. Nat. Clim. Chang. 10 (7), 605-610. http://dx.doi.org/10.1038/s41558-020-0826-9.

Griffiths, B., Fawcett, L., Green, A.C., 2022. Bayesian inference for solar flare extremes. Space Weather. 20 (3), http://dx.doi.org/10.1029/2021SW002886, e2021SW002886, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2021SW002886.

Hamdi, Y., Duluc, C.-M., Rebour, V., 2018. Temperature extremes: Estimation of non-stationary return levels and associated uncertainties. Atmosphere 9 (4), 129. http://dx.doi.org/10.3390/atmos9040129, Number: 4 Publisher: Multidisciplinary Digital Publishing InstituteURL https://www.mdpi.com/2073-4433/9/4/129.

Hausfather, Z., Peters, G.P., 2020. Emissions – the ‘business as usual’ story is misleading. Nature 577 (7792), 618-620. http://dx.doi.org/10.1038/d41586-020-00177-3.

Hawkins, E., Sutton, R., 2009. The potential to narrow uncertainty in regional climate predictions. Bull. Am. Meteorol. Soc. 90 (8), 1095-1108. http://dx.doi.org/10.1175/2009BAMS2607.1, Publisher: American Meteorological Society Section: Bulletin of the American Meteorological Society, URL https://journals.ametsoc.org/view/journals/bams/90/8/2009bams2607_1.xml.

Homan, M.D., Gelman, A., 2014. The no-u-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res. 15 (1), 1593-1623.

Hu, Y., Liang, Z., Peng, A., Wang, K., Wang, J., Li, B., 2023. Nonstationary hydrological frequency analysis in light of model parameters and climate projections uncertainty. J. Hydrol. 617, 129120. http://dx.doi.org/10.1016/j.jhydrol.2023.129120, URL https://www.sciencedirect.com/science/article/pii/S0022169423000628.

Hu, Y., Liang, Z., Singh, V.P., Zhang, X., Wang, J., Li, B., Wang, H., 2018. Concept of equivalent reliability for estimating the design flood under non-stationary conditions. Water Resour Manag. 32 (3), 997-1011. http://dx.doi.org/10.1007/s11269-017-1851-y.

Intergovernmental Panel On Climate Change (Ipcc), 2023. Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, first ed. Cambridge University Press, http://dx.doi.org/10.1017/9781009157896, URL https://www.cambridge.org/core/product/identifier/9781009157896/type/book.

Katz, R.W., Parlange, M.B., Naveau, P., 2002. Statistics of extremes in hydrology. Adv. Water Resour. 25 (8), 1287-1304. http://dx.doi.org/10.1016/S0309-1708(02)00056-8, URL https://www.sciencedirect.com/science/article/pii/S03091708020000568.

Kim, Y.-H., Min, S.-K., Zhang, X., Sillmann, J., Mandstad, M., 2020. Evaluation of the CMIP6 multi-model ensemble for climate extreme indices. Weather. Clim. Extrem. 29, 100269. http://dx.doi.org/10.1016/j.wace.2020.100269, URL https://www.sciencedirect.com/science/article/pii/S2212094719302439.

Laurent, C., Parey, S., 2007. Estimation of 100-year-return-period temperatures in France in a non-stationary climate: Results from observations and IPCC scenarios. Extreme Climatic Events, Glob. Planet. Change Extreme Climatic Events, 57 (1), 177-188. http://dx.doi.org/10.1016/j.gloplacha.2006.11.008. URL https://www.sciencedirect.com/science/article/pii/S0921818106002852.

Laux, P., Weber, E., Feldmann, D., Kunstmann, H., 2023. The robustness of the derived design life levels of heavy precipitation events in the pre-alpine orland region of southern Germany. Atmosphere 14 (9), 1384. http://dx.doi.org/10.3390/atmos14091384, Number: 9 Publisher: Multidisciplinary Digital Publishing Institute, URL https://www.mdpi.com/2073-4433/14/9/1384.

Lehner, F., Deser, C., Maher, N., Marotzke, J., Fischer, E.M., Brunner, L., Knutti, R., Hawkins, E., 2020. Partitioning climate projection uncertainty with multiple large ensembles and CMIP5/6. Earth Syst. Dyn. 11 (2), 491-508. http://dx.doi.org/10.5194/esd-11-491-2020, Publisher: Copernicus GmbH, URL https://esd.copernicus.org/articles/11/491/2020/.

Liang, Z., Hu, Y., Huang, H., Wang, J., Li, B., 2016. Study on the estimation of design value under non-stationary environment. South-To-North Water Transf. Water Sci Tech 14, 50-53.

Locosselli, G.M., Brienen, R.J.W., Leite, M.d., Gloor, M., Krottenthaler, S., Oliveira, A.A.d., Barichivich, J., Anhuf, D., Ceccantini, G., Schöngart, J., Buckeridge, M., 2020. Global tree-ring analysis reveals rapid decrease in tropical tree longevity with temperature. Proc. Natl. Acad. Sci. 117 (52), 33358-33364. http://dx.doi.org/10.1073/pnas.2003873117, Publisher: Proceedings of the National Academy of Sciences, URL https://www.pnas.org/doi/10.1073/pnas.2003873117.

Markova, J., Athanasopoulos, A., Sousa, M.L., Dimova, S., Pinto, A., Iannaccone, S., 2019. State of Harmonised Use of the Eurocodes. Publications Office of the European Union.

Mateus, C., Coonan, B., 2023. Return values of temperature and snow loadings for 50, 100 and 120-year return periods support building design standards in Ireland. In: Advances in Science and Research. 20, Copernicus GmbH, pp. 17-33. http://dx.doi.org/10.5194/asr-20-17-2023, ISSN: 1992-0628, URL https://asr.copernicus.org/articles/20/17/2023/.

Naveau, P., Hannart, A., Ribes, A., 2020. Statistical methods for extreme event attribution in climate science. Annu. Rev. Stat. Appl. 7 (1), 89-110. http://dx.doi.org/10.1146/annurev-statistics-031219-041314, arXiv:https://doi.org/10.1146/annurev-statistics-031219-041314.

Naveau, P., Thao, S., 2022. Multi-model errors and emergence times in climate attribution studies. J. Clim. 4791-4804. http://dx.doi.org/10.1175/JCLI-D-21-0332.1, URL https://journals.ametsoc.org/view/journals/clim/aop/JCLI-D-21-0332.1/JCLI-D-21-0332.1.xml.

Noyelle, R., Zhang, Y., Yiou, P., Faranda, D., 2023. Maximal reachable temperatures for western Europe in current climate. Environ. Res. Lett. http://dx.doi.org/10.1088/1748-9326/acf679, Publisher: IOP Publishing, URL https://hal.science/hal-04201152.

Olsen, J.R., Lambert, J.H., Haimes, Y.Y., 1998. Risk of extreme events under nonstationary conditions. Risk Anal. 18 (4), 497-510. http://dx.doi.org/10.1111/j.1539-6924.1998.tb00364.x, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1539-6924.1998.tb00364.x.

Osborn, T.J., Jones, P.D., Lister, D.H., Morice, C.P., Simpson, I.R., Winn, J.P., Hogan, E., Harris, I.C., 2021. Land surface air temperature variations across the globe updated to 2019: the CRUTEMS data set. J. Geophys. Res.: Atmospheres 126 (2), http://dx.doi.org/10.1029/2019JD032352, e2019JD032352.

Philip, S.Y., Kew, S.F., van Oldenborgh, G.J., Anslow, F.S., Seneviratne, S.I., Vautard, R., Coumou, D., Ebi, K.L., Arrighi, J., Singh, R., van Aalst, M., Pereira Marghidan, C., Wehner, M., Yang, W., Li, S., Schumacher, D.L., Hauser, M., Bonnet, R., Luu, L.N., Lehner, F., Gillett, N., Tradowsky, J.S., Vecchi, G.A., Rodell, C., Stull, R.B., Howard, R., Otto, F.E.L., 2022. Rapid attribution analysis of the extraordinary heat wave on the Pacific coast of the US and Canada in June 2021. Earth Syst. Dyn. 13 (4), 1689-1713. http://dx.doi.org/10.5194/esd-13-1689-2022, Publisher: Copernicus GmbH, URL https://esd.copernicus.org/articles/13/1689/2022/.

Pic, R., Dombry, C., Naveau, P., Taillardat, M., 2025. Proper scoring rules for multivariate probabilistic forecasts based on aggregation and transformation. Adv. Stat. Clim. Meteorol. Ocean. 11 (23-58).
Pörtner, H.O., Bennett, A.F., Bozinovic, F., Clarke, A., Lardies, M.A., Lucassen, M., Pelster, B., Schiemer, F., Stillman, J.H., 2006. Trade-offs in thermal adaptation: the need for a molecular to ecological integration. Physiol Biochem. Zool 79 (2), 295–313. http://dx.doi.org/10.1086/499986.

Read, L.K., Vogel, R.M., 2015. Reliability, return periods, and risk under non-stationarity. Water Resour. Res. 51 (8), 6381–6398. http://dx.doi.org/10.1002/2015WR017089, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/2015WR017089.

Ribes, A., Thao, S., Cattiaux, J., 2020. Describing the relationship between a weather event and climate change: A new statistical approach. J. Clim. 33 (15), 6297–6314. http://dx.doi.org/10.1175/JCLI-D-19-0217.1, Publisher: American Meteorological Society Section: Journal of Climate, URL https://journals.ametsoc.org/view/journals/clim/33/15/JCLI-D-19-0217.1.xml.

Robin, Y., Ribes, A., 2020. Nonstationary extreme value analysis for event attribution combining climate models and observations. Adv. Stat. Clim. Meteorol. Ocean. 6 (2), 205–221. http://dx.doi.org/10.5194/ascmo-6-205-2020, Publisher: Copernicus GmbH, URL https://ascmo.copernicus.org/articles/6/205/2020/.

Rootzén, H., Katz, R.W., 2013. Design life level: Quantifying risk in a changing climate. Water Resour. Res. 49 (9), 5964–5972. http://dx.doi.org/10.1002/wrcr.20425, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/wrcr.20425.

Rust, H.W., Kallache, M., Schellnhuber, H.J., Kropp, J.P., 2011. Confidence intervals for flood return level estimates assuming long-range dependence. In: Kropp, J., Schellnhuber, H.-J. (Eds.), In Extremis: Disruptive Events and Trends in Climate and Hydrology. Springer, Berlin, Heidelberg, pp. 60–88. http://dx.doi.org/10.1007/978-3-642-14863-7_3.

Salas, J.D., Obeysekera, J., 2014. Revisiting the concepts of return period and risk for nonstationary hydrologic extreme events. J. Hydrol. Eng. 19 (3), 554–568. http://dx.doi.org/10.1061/(ASCE)HE.1943-5584.0000820, Publisher: American Society of Civil Engineers, URL https://ascelibrary.org/doi/10.1061/%28ASCE%29HE.1943-5584.0000820.

Serinaldi, F., 2015. Dismissing return periods!. Stoch Env. Res Risk Assess 29 (4), 1179–1189. http://dx.doi.org/10.1007/s00477-014-0916-1.

Smith, R.L., 1999. Bayesian and frequentist approaches to parametric predictive inference. In: Bernardo, J.M., Berger, J.O., Dawid, A.P., Smith, A.F.M. (Eds.), Bayesian Statistics 6. Oxford University PressOxford, pp. 586–612. http://dx.doi.org/10.1093/oso/9780198504856.003.0026.

Stan Development Team, 2024. Stan modeling language users guide and reference manual, 2.35. URL https://mc-stan.org/users/documentation/.

Toimil, A., Camus, P., Losada, I.J., Alvarez-Cuesta, M., 2021. Visualising the uncertainty cascade in multi-ensemble probabilistic coastal erosion projections. Front. Mar. Sci. 8, http://dx.doi.org/10.3389/fmars.2021.683535, Publisher: Frontiers, URL https://www.frontiersin.org/journals/marine-science/articles/10.3389/fmars.2021.683535/full.

Van de Vyver, H., 2015. Bayesian estimation of rainfall intensity-duration-frequency relationships. J. Hydrol. 529, 1451–1463. http://dx.doi.org/10.1016/j.jhydrol.2015.08.036, URL https://www.sciencedirect.com/science/article/pii/S0022169415006083.

Vardi, Y., Zhang, C.-H., 2000. The multivariate $l_1$-median and associated data depth. Proc. Natl. Acad. Sci. 97 (4), 1423–1426. http://dx.doi.org/10.1073/pnas.97.4.1423, Company: The National Academy of Sciences Distributor: The National Academy of Sciences Institution: The National Academy of Sciences Label: The National Academy of Sciences Publisher: Proceedings of the National Academy of SciencesURL https://www.pnas.org/doi/abs/10.1073/pnas.97.4.1423.

Wigley, T.M.L., 2009. The effect of changing climate on the frequency of absolute extreme events. Clim. Change 97 (1), 67–76. http://dx.doi.org/10.1007/s10584-009-9654-7.

Yan, L., Xiong, L., Guo, S., Xu, C.-Y., Xia, J., Du, T., 2017. Comparison of four nonstationary hydrologic design methods for changing environment. Investigation of Coastal Aquifers, J. Hydrol. Investigation of Coastal Aquifers, 551, 132–150. http://dx.doi.org/10.1016/j.jhydrol.2017.06.001. URL https://www.sciencedirect.com/science/article/pii/S002216941730392X.

Zeder, J., Sippel, S., Pasche, O.C., Engelke, S., Fischer, E.M., 2023. The effect of a short observational record on the statistics of temperature extremes. Geophys. Res. Lett. 50 (16), http://dx.doi.org/10.1029/2023GL104090, e2023GL104090, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2023GL104090.

Zeimbekakis, A., Schifano, E.D., and, J.Y., 2024. On misuses of the Kolmogorov-Smirnov test for one-sample goodness-of-fit. Amer. Statist. 78 (4), 481–487. http://dx.doi.org/10.1080/00031305.2024.2356095, arXiv:https://doi.org/10.1080/00031305.2024.2356095.