Improved Correction of Extreme Precipitation Through Explicit and Continuous Nonstationarity Treatment and the Metastatistical Approach

Abstract

Climate models simulate extreme precipitation under nonstationarity due to continuous climate change. However, systematic errors in local‐scale climate projections are often corrected using stationary or quasi‐stationary methods without explicit and continuous nonstationarity treatment, like quantile mapping (QM), detrended QM, and quantile delta mapping. To bridge this gap, we introduce nonstationary QM (NS‐QM) and its simplified version for consistent nonstationarity patterns (CNS‐QM). Besides, correction approaches for extremes often rely on limited extreme‐event records. To leverage ordinary‐event information while focusing on extremes, we propose integrating the simplified Metastatistical extreme value (SMEV) distribution into NS‐QM and CNS‐QM (NS‐QM‐SMEV and CNS‐QM‐SMEV). We demonstrate the superiority of NS‐ and CNS‐QM‐SMEV over existing methods through a simulation study and show several real‐world applications using high‐resolution‐regional and coarse‐resolution‐global climate models. NS‐QM and CNS‐QM reflect nonstationarity more realistically but may encounter challenges due to data limitations like estimation errors and uncertainty, particularly for the most extreme events. These issues, shared by existing approaches, are effectively mitigated using the SMEV distribution. NS‐ and CNS‐QM‐SMEV offer lower estimation error, approximate unbiasedness, reduced uncertainty, and improved representation of the entire distribution, especially for samples of ∼70 years, and greater superiority with larger samples. We show existing methods may perform competitively for short samples but exhibit substantial biases in quantile‐quantile matching due to bypassing nonstationarity modeling. NS‐ and CNS‐QM‐SMEV avoid these biases, adhering better to their theoretical functioning. Thus, NS‐ and CNS‐QM‐SMEV enhance the correction of extremes under nonstationarity. Yet, properly identifying nonstationarity patterns is crucial for reliable implementations.

Introduction

Climate model projections are used for evaluating future precipitation extremes under climate change (Emmanouil et al., 2023; Moustakis et al., 2021; Pfahl et al., 2017; Rajczak & Schär, 2017). However, these projections often exhibit biases stemming from various sources, including model errors, areal averaging effects, and the use of parameterizations instead of resolving sub-grid processes (e.g., convection) (Bony et al., 2015; Mitovski et al., 2019; Sherwood et al., 2014). Therefore, these projections require post-processing to improve their applicability in local-scale studies, such as hazard assessments and climate change impact studies (Brunner et al., 2021; Christensen et al., 2008; Teutschbein & Seibert, 2012). However, post-processing is a key source of uncertainty in local projections of extremes (Lafferty & Sriver, 2023; Michalek et al., 2024; Wootten et al., 2017; Yuan et al., 2017), alongside internal variability (due to the chaotic Earth system), scenario uncertainty (possible greenhouse gas emissions trajectories), and climate model uncertainty (stemming from model structures and parameters). While post-processing uncertainty in precipitation varies by region, it is often substantial. Wootten et al. (2017) showed that it can account for up to ∼20% of total uncertainty for the southeastern US, or even ∼30–50% in complex terrain or when ensemble projections are limited. Lafferty and Sriver (2023) found that it often exceeds model and scenario uncertainties in general. Thus, addressing post-processing uncertainty is key to improving local climate projections.

Post-processing hydro-climatological extremes presents notable challenges (Hernanz et al., 2022), especially in establishing deterministic links between precipitation extremes and other atmospheric variables. Bias correction, an effective model output statistics approach, establishes a transfer function between the statistics of simulations and observations in a baseline period to correct simulations (Hertig et al., 2019; Lange, 2019; Maraun & Widmann, 2017; Rummukainen, 1997). In particular, quantile mapping (QM) (Panofsky & Brier, 1968) is a popular correction technique that adjusts the full distribution of simulations and is superior to other methods that only adjust the mean or variance (Teutschbein & Seibert, 2012).

Nonstationarity, which refers to the deterministic evolution of stochastic process statistics over time, challenges bias correction methods. Climate change and human activities may induce nonstationarity in hydro-climatic extremes (Milly et al., 2008, 2015). While climate models aim to depict the continuously changing Earth system and inherent nonstationarity, QM assumes stationarity, undermining its reliability for correcting climate model simulations. Despite existing variations of QM account for relative changes between baseline and prediction periods, they still assume stationarity in each period. For instance, the detrended QM (DQM) method (Bürger et al., 2013) applies QM to simulations by removing and reintroducing changes in the mean in the prediction period relative to the baseline period. Similarly, the quantile delta mapping (QDM) method (Cannon et al., 2015) preserves changes from the baseline to the prediction period across all quantiles. The cumulative distribution function transform (CDFt) (Michelangeli et al., 2009) relates stationary CDFs of simulations and observations in the baseline period to correct the simulations' CDF in the prediction period with a change between periods. Kallache et al. (2011) introduced a CDFt variant that uses concatenated Generalized Pareto (GP) distributions with stochastic covariates, leading to doubly stochastic estimates (e.g., see Serinaldi et al., 2018; Serinaldi & Kilsby, 2015). Since climate change occurs continuously over time rather than as a sudden change at a particular time, these methods are quasi-stationary and cannot capture the nonstationarity realistically. Omitting nonstationary modeling may distort the nonstationary signals and affect estimates (Cannon et al., 2015; Maraun, 2016; Maurer & Pierce, 2014). Besides, different degrees of nonstationarity between observations and simulations can lead to nonstationary biases (Hui et al., 2020; Merkenschlager et al., 2017; Van De Velde et al., 2022), introducing additional errors. Yet, existing methods cannot deal with these issues due to their lack of suitable nonstationarity treatment.

Furthermore, there are two general approaches to correct extremes using QM or its variants, each with pros and cons. One approach involves correcting all simulations and extracting extremes (e.g., annual maxima). However, this approach often employs distributions tailored for ordinary events, resulting in a poor representation of the distribution of extremes and substantial errors (Cannon et al., 2015; Heo et al., 2019; Maraun, 2013). The other approach only employs extreme-event data, such as annual maxima, and is more popular (e.g., Kim et al., 2021; Srivastav et al., 2014). This approach commonly employs distributions from Extreme Value Theory, such as the Generalized Extreme Value (GEV) (Kim et al., 2021; Srivastav et al., 2014; Um et al., 2016). Yet, since QM and its variants use separate distributions for simulations and observations in baseline and prediction periods, relying on multiple distributions parametrized with often scarce observations of extremes can impede reliable bias correction. Thus, there is a need for a correction approach that integrates the strengths of these two general approaches by harnessing information beyond extremes while focusing on extremes.

Recently, extreme value distributions leveraging ordinary-event records have gained popularity in hydrological frequency analysis due to their superior predictive capability (Boumis et al., 2024; Falkensteiner et al., 2023; Gründemann et al., 2023; Marra et al., 2018; Miniussi et al., 2020; Vidrio-Sahagún et al., 2023). In particular, the simplified Metastatistical extreme value (SMEV) distribution (Marra et al., 2019) integrates ordinary independent events along with the extremes and can incorporate nonstationary structures. The nonstationary SMEV distribution offers comparable performance to other nonstationary ordinary-event-based distributions, is easy to implement, and generally outperforms the GEV distribution (Vidrio-Sahagún et al., 2023; Vidrio-Sahagún & He, 2022a). Both frequency analysis and bias correction of extremes share the use of extreme value distributions. Yet, the SMEV distribution has not been explored for correcting extreme event projections. Therefore, incorporating the SMEV distribution for correcting extremes is appealing, as it would harness ordinary-event data, which are often more abundant in observations and better simulated by climate models than the extremes.

Acknowledging the current shortcomings in addressing nonstationarity and using multiple distributions with limitations of extreme data, we propose a novel approach for bias-correcting precipitation extremes. To the best of our knowledge, this is the first approach that (a) incorporates explicit and continuous nonstationary treatment into the QM method and (b) leverages information on extremes embedded in ordinary events using the SMEV distribution. We systematically assess and compare the proposed approach with several widely used methods (including QM, DQM, QDM, and CDF-t) using the GEV distribution as benchmarks through a simulation study. In addition, we demonstrate its practical application using observational data sets and their simulated counterparts from high-resolution regional and coarse-resolution global climate models. The methods were evaluated using corrected distributions, quantiles, and time series, considering similarity with actual distributions, accuracy, and uncertainty.

Proposed Approach

To correct biases in nonstationary extreme precipitation simulated by climate models, we propose the nonstationary QM (NS-QM) and its simplified version for consistent nonstationarity patterns (CNS-QM), both of which address nonstationarity continuously and explicitly. NS-QM can deal with different nonstationarity patterns in the baseline and prediction periods. CNS-QM is applicable when nonstationarity patterns in both periods are consistent, which could be the case in some climate model projections. Note that NS-QM can also deal with consistent nonstationarity patterns, albeit with unnecessary complexity. Figure 1 outlines the nonstationarity treatment of the proposed methods, which offer greater realism than existing stationary (QM) and quasi-stationary (DQM, QDM, and CDFt) methods. As shown in this figure, quasi-stationary methods oversimplify nonstationarity by assuming step changes between periods and stationarity within each period (for more details on these methods, refer to the Appendix A). In contrast, the proposed methods address nonstationarity explicitly considering its continuous and deterministic nature. On the other hand, the proposed and existing variants of the QM method employ distinct probability distributions for simulations and observations in baseline and prediction periods (e.g., see Figure 2). Using multiple distributions with limited data might cause large errors and uncertainty. Therefore, we propose integrating SMEV distributions into the correction process to use more information beyond the extremes and enhance prediction accuracy and certainty.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

A key advantage of continuous nonstationary correction in impact analysis is its ability to account for gradual and sustained changes in climate variables across the planning period. For instance, when heavy precipitation increases, continuous nonstationary correction prevents early overestimation and late underestimation in the planning horizon by avoiding the oversimplification of averaging changes across the period. Additionally, it enables incremental adaptive measures, such as raising riverbanks or implementing decentralized stormwater management practices as needed progressively over time, rather than implementing measures for distant future conditions at once. This can facilitate cost-effective adjustments to mitigate future hydroclimate risks, enhancing adaptation efforts amid various resource constraints (e.g., financial constraints). Ultimately, continuous nonstationary correction can provide stakeholders with more information to develop resource-efficient adaptation strategies.

Nonstationary Bias Correction of Extremes

The conventional QM method (Panofsky & Brier, 1968) relates the cumulative distribution functions $$H(\mathit{\cdot })$$ of observations and simulations by a transfer function. The calibrated transfer function is then applied to the raw simulations in the bias correction process. The QM method assumes stationarity of the $$H(\mathit{\cdot })$$ of both observations and simulations across baseline and prediction periods, and the corrected simulations of annual maxima $$Y$$ are estimated as: 1 $\begin{align*}{\widehat{y}}_{t}^{s,b}={H}_{{Y}^{o,b}}^{-1}\left({H}_{{Y}^{s,b}}\left({y}_{t}^{s,b}\right)\right)\\ {\widehat{y}}_{t}^{s,p}={H}_{{Y}^{o,b}}^{-1}\left({H}_{{Y}^{s,b}}\left({y}_{t}^{s,p}\right)\right)\end{align*}$ where the superscripts of $$Y$$ provide the full specification of $$H(\cdot )$$; $$o$$ and $$s$$ denote observations and simulations, respectively, while $$b$$ and $$p$$ indicate baseline and prediction periods, respectively; $${y}_{t}^{s,(\cdot )}$$ denotes the simulated series; and $${\widehat{y}}_{t}^{s,(\cdot )}$$ is the corrected $${y}_{t}^{s,(\cdot )}$$.

However, assuming stationarity is unrealistic if the statistics of extremes continuously change over time in response to external drivers, such as global warming. For such circumstances, we propose a new formulation of the QM method, namely NS-QM, that explicitly deals with nonstationary extremes. In NS-QM, observations and simulations are considered realizations from nonstationary stochastic processes denoted as $${Y}_{t}^{o}$$ and $${Y}_{t}^{s}$$, respectively, whose statistics continuously change over time while governed by a deterministic function. As a result, $${Y}_{t}^{o}$$ and $${Y}_{t}^{s}$$ are characterized by the nonstationary distributions $${H}_{{Y}_{t}^{o}}\left(\cdot ;{\boldsymbol{v}}_{\boldsymbol{t}}\right)$$ and $${H}_{{Y}_{t}^{s}}\left(\cdot ;{\boldsymbol{v}}_{\boldsymbol{t}}\right)$$, respectively, both governed by the covariate vector $${\boldsymbol{v}}_{\boldsymbol{t}}$$. Here, $${v}_{t}$$ denotes a generic covariate, as NS-QM is not restricted to any specific covariate. In practice, $${v}_{t}$$ could be time or a physical variable, provided that nonstationarity can be successfully attributed to it and a deterministic law of time can be identified. To account for potential changes in the nonstationarity patterns between subperiods, NS-QM incorporates the CDFt principle (Michelangeli et al., 2009) through the transformation $$T:[0,1]\to [0,1]$$ that links $${H}_{{Y}_{t}^{o}}\left(\cdot ;{\boldsymbol{v}}_{\boldsymbol{t}}\right)$$ and $${H}_{{Y}_{t}^{s}}\left(\cdot ;{\boldsymbol{v}}_{\boldsymbol{t}}\right)$$. Such a transformation can be expressed as $${H}_{{Y}_{t}^{o,b}}\left({y}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)=T\left({H}_{{Y}_{t}^{s.b}}\left({y}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)\right)$$. Besides, the time-dependent variable $${u}_{t}={H}_{{Y}_{t}^{s.b}}\left({y}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)$$ is defined, such that $${u}_{t}\mathit{\in }[0,1]\,\mathit{\forall }\,t$$, and thus $${y}_{t}={H}_{{Y}_{t}^{s.b}}^{-1}\left({u}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)$$. Hence, the transformation can be expressed as $$T\left({u}_{t}\right)={H}_{{Y}_{t}^{o,b}}\left({H}_{{Y}_{t}^{s.b}}^{-1}\left({u}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right);{\boldsymbol{v}}_{\boldsymbol{t}}\right)$$. Therefore, assuming that $$T(\cdot )$$ remains valid in the prediction period (i.e., $${H}_{{Y}_{t}^{o,p}}\left({y}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)=T\left({H}_{{Y}_{t}^{s,p}}\left({y}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)\right)$$), the CDF of $${\widehat{y}}_{t}^{s,p}$$ is estimated by: 2 $\begin{align*}{H}_{{\widehat{Y}}_{t}^{s,p}}\left({y}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)\cong {H}_{{Y}_{t}^{o,p}}\left({y}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)=T\left({H}_{{Y}_{t}^{s,p}}\left({y}_{t};{\boldsymbol{v}}_{\boldsymbol{t}}\right)\right)={H}_{{Y}_{t}^{o,b}}\left({H}_{{Y}_{t}^{s,b}}^{-1}\left[{H}_{{Y}_{t}^{s,p}}\left({y}_{t}^{s,p};{\boldsymbol{v}}_{\boldsymbol{t}}\right);{\boldsymbol{v}}_{\boldsymbol{t}}\right];{\boldsymbol{v}}_{\boldsymbol{t}}\right)\end{align*}$

Note that in Equation 2, the nonstationary distributions $${H}_{{Y}_{t}^{s,p}}(\cdot )$$, $${H}_{{Y}_{t}^{s,b}}(\mathit{\cdot })$$, and $${H}_{{Y}_{t}^{o,b}}(\cdot )$$, correspond to different periods but are used at concurrent time points in the prediction period through extrapolation of $${H}_{{Y}_{t}^{s,b}}(\cdot )$$ and $${H}_{{Y}_{t}^{o,b}}(\cdot )$$. In this way, the three distributions are virtually present across the same period, allowing the application of $$T(\cdot )$$ at every point in time. For bias-correcting $${y}_{t}^{s,p}$$, numerical quantile-quantile mapping is performed between $${H}_{{\widehat{Y}}_{t}^{s,p}}(\cdot )$$ and $${H}_{{Y}_{t}^{s,p}}(\cdot )$$. This consists of finding the quantiles $${\widehat{y}}_{t}^{s,p}$$ whose probabilities $${H}_{{\widehat{Y}}_{t}^{s,p}}(\cdot )$$ are equal to those of $${y}_{t}^{s,p}$$ in $${H}_{{Y}_{t}^{s,p}}(\cdot )$$. Such quantile-quantile mapping can be expressed as: 3 $\begin{align*}{\widehat{y}}_{t}^{s,p}={H}_{{\widehat{Y}}_{t}^{s,p}}^{-1}\left({H}_{{Y}_{t}^{s,p}}\left[{y}_{t}^{s,p};{\boldsymbol{v}}_{\boldsymbol{t}}\right];{\boldsymbol{v}}_{\boldsymbol{t}}\right)\end{align*}$

Figures 2a and 2b graphically illustrate the derivation of $${H}_{{\widehat{Y}}_{t}^{s,p}}(\cdot )\cong {H}_{{Y}_{t}^{o,p}}(\cdot )$$ in the NS-QM method. Furthermore, it is worth noting that when correcting simulations in the baseline period, the transformation $$T(\cdot )$$ is unnecessary and Equation 2 reduces to $${\widehat{y}}_{t}^{s,b}={H}_{{Y}_{t}^{o,b}}^{-1}\left({H}_{{Y}_{t}^{s,b}}\left[{y}_{t}^{s,b};{\boldsymbol{v}}_{\boldsymbol{t}}\right];{\boldsymbol{v}}_{\boldsymbol{t}}\right)$$.

When nonstationarity patterns do not change from the baseline to the prediction period, namely, showing consistent nonstationarity patterns across these periods, the $$T(\cdot )$$ is no longer required. As a result, NS-QM can be simplified to consider a consistent nonstationarity pattern, referred to as the CNS-QM method, which is given by: 4 $\begin{align*}{\widehat{y}}_{t}^{s}={H}_{{Y}_{t}^{o}}^{-1}\left({H}_{{Y}_{t}^{s}}\left[{y}_{t}^{s};{\boldsymbol{v}}_{\boldsymbol{t}}\right]\right)\end{align*}$ where $${H}_{{Y}_{t}^{o}}$$ and $${H}_{{Y}_{t}^{s}}$$ can be estimated using the entire data sets of observations and simulations without distinguishing between baseline and prediction periods, especially for $${H}_{{Y}_{t}^{s}}$$. Compared to NS-QM, the CNS-MQ method offers the advantage of using fewer distributions, albeit at the expense of assuming invariant nonstationarity patterns across baseline and prediction periods. The CNS-QM method is graphically illustrated in Figures 2c and 2d.

It is worth noting that bias correction is meaningful when model outputs are informative and realistic to some degree despite systematic biases. A key requirement for the proposed nonstationary correction methods is that climate models capture nonstationarity qualitatively, that is, by showing consistent trend directions in both observations and simulations, even if trend magnitudes are biased.

Bias Correction of Extremes Leveraging Ordinary Event Data

The GEV distribution is a popular choice for $$H(\cdot )$$ when bias-correcting extremes. In the NS-QM and CNS-QM methods, we incorporate parsimonious nonstationary structures into the GEV distribution to capture the nonstationarity, expressing the location and/or scale distribution parameters as linear functions of the selected $${v}_{t}$$ as follows: 5 $\begin{align*}{GEV}_{k,l,0}\sim \mathrm{G}EV\left({\xi }_{t},{\alpha }_{t},\kappa \right);\,\left\{\begin{array}{@{}c@{}}{\xi }_{t}={\xi }_{0}+\mathbb{I}(k){\xi }_{1}{v}_{t}\\ {\alpha }_{t}={\alpha }_{0}+\mathbb{I}(l){\alpha }_{1}{v}_{t}\\ \kappa =\text{constant}\end{array}\right.\end{align*}$ where $$\xi $$, $$\alpha $$, and $$\kappa $$ are the location, scale, and shape parameters of the GEV distribution; $$\mathbb{I}(\cdot )$$ is the indicator function, which takes the value of 1 if its argument is “1” and 0 otherwise; $$k$$ and $$l$$ indicate whether $$\xi $$ and $$\alpha $$ are expressed as a function of $${v}_{t}$$, respectively; and $${\xi }_{0}$$ and $${\xi }_{1}$$, and $${\alpha }_{0}$$ and $${\alpha }_{1}$$ are the regression coefficients of $$\xi $$ and $$\alpha $$, respectively. The $$\kappa $$ was treated as constant, as it is unrealistic to allow it to vary due to sample size constraints in practical applications. Thereby, we consider the two model structures, $${GEV}_{1,0,0}$$ and $${GEV}_{1,1,0}$$, in this paper. For instance, in $${GEV}_{1,1,0}$$, $$\xi $$ and $$\alpha $$ both vary with $${v}_{t},$$ as $${\xi }_{t}={\xi }_{0}+{\xi }_{1}{v}_{t}$$, and $${\alpha }_{t}={\alpha }_{0}+{\alpha }_{1}{v}_{t}$$.

To incorporate the information from ordinary events into the NS-QM and CNS-QM methods, we substitute the GEV distribution with an ordinary-event-based distribution of extremes. Specifically, we adopt the SMEV distribution (Marra et al., 2019) for simplicity, as it can perform equivalently to the more complex nonstationary MEV distribution (Vidrio-Sahagún et al., 2023). The SMEV distribution $${\Psi }$$ of block maxima relies on the parameter vector $$\boldsymbol{\theta }$$ of the distribution $$F(\cdot )$$ describing the magnitudes of ordinary events $$X$$ and the number of events per block ($$n$$). In this paper, we used the Weibull distribution $$\mathcal{W}(\cdot )$$ for $$F(\cdot )$$ due to its theoretical and empirical support in describing rainfall tails (De Michele & Avanzi, 2018; Marra et al., 2023; Serinaldi & Kilsby, 2014; Wilson & Toumi, 2005). The stationary and nonstationary $${\Psi }$$ distributions are given, respectively, by: 6 $\begin{align*}{\Psi }(x;\boldsymbol{\theta },n)={[\mathcal{W}(x;\boldsymbol{\theta })]}^{n}={\left[1-\mathrm{exp}\left({-\left(\frac{x}{\lambda }\right)}^{\beta }\right)\right]}^{n}\end{align*}$ 7 $\begin{align*}{\Psi }\left(x;{\boldsymbol{\theta }}_{t},{n}_{t}\right)={\left[\mathcal{W}\left(x;{\boldsymbol{\theta }}_{t}\right)\right]}^{{n}_{t}}={\left[1-\mathrm{exp}\left({-\left(\frac{x}{{\lambda }_{t}}\right)}^{{\beta }_{t}}\right)\right]}^{{n}_{t}}\end{align*}$ where $$\boldsymbol{\theta }$$ is composed of $$\lambda $$ and $$\beta $$, which are the scale and shape parameters of the $$\mathcal{W}(\cdot )$$ distribution, respectively. In NS-QM and CNS-QM, $${H}_{{Y}_{t}^{(\cdot )}}(\cdot )$$ in Equations 2–4 is replaced by $${{\Psi }}_{{X}_{t}^{(\cdot )}}$$, denoting the SMEV distribution. We adopt parsimonious nonstationary structures where the distribution parameters are expressed as linear functions of $${v}_{t}$$ as given below: 8 $\begin{align*}{SMEV}_{q,r,s}\sim \mathrm{S}MEV\left({\lambda }_{t},{\beta }_{t},{n}_{t}\right);\,\left\{\begin{array}{@{}c@{}}{\lambda }_{t}={\lambda }_{0}+\mathbb{I}(q){\lambda }_{1}{v}_{t}\\ \begin{array}{c}{\beta }_{t}={\beta }_{0}+\mathbb{I}(r){\beta }_{1}{v}_{t}\\ {n}_{t}={n}_{0}+\mathbb{I}(s){n}_{1}{v}_{t}\end{array}\end{array}\right.\end{align*}$ where $$q$$, $$r$$, and $$s$$ are model-specific subindices used to denote the nonstationary structure and are set to “0” or “1” to indicate whether $$\lambda $$, $$\beta $$, and $$n$$ are expressed as a function of $${v}_{t}$$ or time-invariant, respectively; and $${\lambda }_{0}$$ and $${\lambda }_{1}$$, $${\beta }_{0}$$ and $${\beta }_{1}$$, as well as $${n}_{0}$$ and $${n}_{1}$$ are the regression coefficients of $$\lambda $$, $$\beta $$, and $$n$$, respectively. Following this notation, the four nonstationary SMEV structures are $$S{MEV}_{1,0,0}$$, $${SMEV}_{1,1,0}$$, $${SMEV}_{1,0,1}$$ and $${SMEV}_{1,1,1}$$. For instance, in $${SMEV}_{1,1,0}$$, $$\lambda $$ and $$\beta $$ vary with $${v}_{t}$$, as $${\lambda }_{t}={\lambda }_{0}+{\lambda }_{1}{v}_{t}$$ and $${\beta }_{t}={\beta }_{0}+{\beta }_{1}{v}_{t}$$, while $$n$$ remains constant.

Parameter Estimation and Uncertainty Quantification

The parameters of the GEV distribution and the $$\mathcal{W}(X)$$ distribution in the SMEV distributions are estimated using the maximum likelihood method due to its easiness of incorporating the nonstationary structure into the probability distribution. Additionally, for the SMEV distribution, $$n$$ is estimated as the sample average of $${n}_{j}$$ when the SMEV-model-specific subindex $$s=0$$, and $${n}_{0}$$ and $${n}_{1}$$ for $${n}_{t}$$ are estimated using the Theil-Sen trend estimator when $$s=1$$ (Equation 8).

In practice, the low-magnitude ordinary events may negatively impact the estimation of the right tail of $$\mathcal{W}(X)$$ in the SMEV distribution (Marra et al., 2019; Miniussi & Marra, 2021). Thus, we employ the left-censored version of the maximum likelihood method (Cohn, 2005; Helsel, 2011; Vidrio-Sahagún & He, 2022a) for estimating the parameters of $$\mathcal{W}(X)$$ in practical applications (Section 3.3). The threshold for censoring the data is case-specific and is commonly determined by optimizing some relevant metric (Marra et al., 2019, 2020; Vidrio-Sahagún & He, 2022a). The threshold is thus specified by maximizing the fitting efficiency (Akaike Information Criterion [AIC]) while ensuring the suitability of the $$\mathcal{W}(X)$$ distribution for the left-censored data. Furthermore, the estimation of $$n$$ or $${n}_{t}$$ in the SMEV distribution also accounts for the drizzle effect of climate and atmospheric models, which refers to the overestimation of the occurrence and duration of precipitation events with very low intensity (Chen et al., 2021; Stephens et al., 2010; Trenberth, 2011). Hence, we adopt a trace precipitation threshold of 1 mm/day to define days with precipitation events in practical applications, as in several previous studies (Gomez-Garcia et al., 2019; Polade et al., 2014; Vincent et al., 2018; Zhang et al., 2011).

We quantify the estimation uncertainty of the methods using the parametric bootstrap (Efron et al., 1992) and derive the confidence intervals of the bias-corrected estimates at a $$\rho $$ significance level as follows:

• First, we obtain parametric estimates $${\widehat{H}}_{{Y}_{t}^{s,b}}$$, $${\widehat{H}}_{{Y}_{t}^{s,p}}$$, and $${\widehat{H}}_{{Y}_{t}^{o,b}}$$ (when using the GEV distribution) or $${\widehat{{\Psi }}}_{{X}_{t}^{s,b}}$$, $${\widehat{{\Psi }}}_{{X}_{t}^{s,p}}$$, and $${\widehat{{\Psi }}}_{{X}_{t}^{o,b}}$$ (when adopting the SMEV distribution) using the available samples of simulations and observations.

• Next, we generate $${N}_{samp}^{b}$$ synthetic samples of size $$n$$ from $${\widehat{H}}_{{Y}_{t}^{s,b}}$$, $${\widehat{H}}_{{Y}_{t}^{s,p}}$$, and $${\widehat{H}}_{{Y}_{t}^{o,b}}$$ (when using the GEV) or of size $$n\times M$$ from the ordinary-event distribution of $${\widehat{{\Psi }}}_{{X}_{t}^{s,b}}$$, $${\widehat{{\Psi }}}_{{X}_{t}^{s,p}}$$, and $${\widehat{{\Psi }}}_{{X}_{t}^{o,b}}$$ (when adopting the SMEV).

• We obtain parametric estimates $${\widehat{H}}_{{Y}_{t}^{s,b}}^{i}$$, $${\widehat{H}}_{{Y}_{t}^{s,p}}^{i}$$, and $${\widehat{H}}_{{Y}_{t}^{o,b}}^{i}$$ (when using the GEV) or $${\widehat{{\Psi }}}_{{X}_{t}^{s,b}}^{i}$$, $${\widehat{{\Psi }}}_{{X}_{t}^{s,p}}^{i}$$, and $${\widehat{{\Psi }}}_{{X}_{t}^{o,b}}^{i}$$ (when using the SMEV) for each synthetic sample $$i=1,2,\mathit{\dots },{N}_{samp}^{b}$$.

• Then, we estimate the $${N}_{samp}^{b}$$ sets of bias-adjusted quantiles $${\widehat{y}}_{t}^{s,{(\cdot )}^{i}}$$ based on all the parametric estimates $${\widehat{H}}_{{Y}_{t}^{(\cdot )}}^{i}$$ or $${\widehat{{\Psi }}}_{{X}_{t}^{(\cdot )}}^{i}$$.

• Finally, we derive the $$(1-\rho $$) $$100\mathit{\%}$$ confidence intervals using the $$\rho /2$$ and $$1-\rho /2$$ percentiles from all the sets of $${\widehat{y}}_{t}^{s,{(\cdot )}^{i}}$$.

Note that in the CNS-QM method, the $${\widehat{H}}_{{Y}_{t}^{s,b}}$$ and $${\widehat{H}}_{{Y}_{t}^{s,p}}$$ or $${\widehat{{\Psi }}}_{{X}_{t}^{s,b}}$$ and $${\widehat{{\Psi }}}_{{X}_{t}^{s,p}}$$ are combined into $${\widehat{H}}_{{Y}_{t}^{s}}$$ or $${\widehat{{\Psi }}}_{{X}_{t}^{s}}$$, respectively. We use a significance level of $$\rho =0.05$$ for the estimation of confidence intervals and all the statistical analyses employed in this paper.

Implementation and Evaluation of Proposed Methods

We assess and compare the proposed NS-QM and CNS-QM with benchmark methods, including the conventional QM, DQM, QDM, and CDFt in both a simulation study and several real-world applications.

Simulation Study

We assess the performance of the proposed methods through comparison with the QM, DQM, QDM, and CDFt methods as benchmarks in a simulation study, where the statistical characteristics of the underlying stochastic processes are known a priori. We generate synthetic daily series representing observations and simulations by randomly sampling $$N$$ times per block from two $$\mathcal{W}(\cdot )$$ distributions with time-dependent $${\lambda }_{t}$$ and $${\beta }_{t}$$ with different parametrizations. Thus, distinct distributions are used for observations and simulations, which reflect the biases in simulations and so align well with reality. Here, $$N$$ is kept constant, as its variations typically play a secondary role in the context of nonstationary extremes (Vidrio-Sahagún et al., 2023). Two nonstationary scenarios are considered. In one scenario, the nonstationarity is consistent across the baseline and prediction periods, that is, the trend slopes in observations and simulations are constant. In the other scenario, the trend slopes change by intensifying/increasing in the prediction period with respect to the baseline period. Both scenarios consider that simulations underestimate the degree of nonstationarity (having weaker trends), and are defined as follows: 9 $$\begin{array}{lll}\text{Consistent}\,\text{nonstationarity}\end{array}:\left\{\begin{array}{@{}l@{}}{X}_{t}\sim \mathcal{W}\left({\lambda }_{t},{\beta }_{t}\right)\quad \text{where}\,\left\{\begin{array}{@{}c@{}}{\lambda }_{t}={\lambda }_{0}+{\lambda }_{1}{v}_{t}\\ {\beta }_{t}={\beta }_{0}+{\beta }_{1}{v}_{t}\end{array}\right.\,\\ N\sim \mathcal{N}\left({\mu }_{{N}_{0}},{\sigma }_{{N}_{0}}\right)\quad \text{where}\,\left\{\begin{array}{@{}l@{}}{\mu }_{{N}_{0}}=\text{constant}\\ {\sigma }_{{N}_{0}}=\text{constant}\end{array}\right.\end{array}\right.$$ 10 $$\begin{array}{c}\text{Changing}\text{nonstationarity}\end{array}:\left\{\begin{array}{@{}l@{}}{X}_{t}\sim \mathcal{W}\left({{\Lambda }}_{t},{\mathrm{\mathrm{B}}}_{t}\right)\quad \text{where}\,\left\{\begin{array}{@{}c@{}}\begin{array}{c}{\lambda }_{t}={\lambda }_{0}+{\lambda }_{1}{v}_{t}\quad \forall t\in b\\ {\lambda }_{t}={\lambda }_{0}+{\lambda }_{2}{v}_{t}-{{\increment}}_{{\lambda }_{b:p}}\ \forall t\in p\end{array}\\ \begin{array}{c}{\beta }_{t}={\beta }_{0}+{\beta }_{1}{v}_{t}\quad \forall t\in b\\ {\beta }_{t}={\beta }_{0}+{\beta }_{2}{v}_{t}-{{\increment}}_{{\beta }_{b:p}}\ \forall t\in p\end{array}\end{array}\right.\\ \begin{array}{cc}N\sim \mathcal{N}\left({\mu }_{{N}_{0}},{\sigma }_{{N}_{0}}\right)& \text{where}\,\left\{\begin{array}{@{}c@{}}{\mu }_{{N}_{0}}=\text{constant}\\ {\sigma }_{{N}_{0}}=\text{constant}\end{array}\right.\end{array}\,\end{array}\right.$$ where $${\mathit{{\increment}}}_{{{\Lambda }}_{b:p}}$$ and $${\mathit{{\increment}}}_{{\mathrm{\mathrm{B}}}_{b:p}}$$ are adjustment terms for filling the differences in $${\lambda }_{t}$$ and $${\beta }_{t}$$ estimated using coefficients $${\lambda }_{2}$$ and $${\beta }_{2}$$ instead of $${\lambda }_{1}$$ and $${\beta }_{1}$$ at the end of the baseline to ensure distribution continuity in the transition from the baseline to the prediction period; $$\mathcal{N}(\cdot )$$ denotes the Normal distribution; and $${v}_{t}$$ is a hypothetical covariate that monotonically and linearly increases over time, representing nonstationarity deterministically and mimicking the climate change signal after removing the natural variability envelope. Thus, in the simulation study, $${v}_{t}$$ does not correspond to any specific real-world variable or entity. We used this hypothetical $${v}_{t}$$ to simplify the assessment of methods' performance without the noise and complexity of nonstationary stochastic climate variables (e.g., air temperature), while keeping practical generalizability. Table 1 shows the hyperparameters used in both nonstationary scenarios, which make the values of $$\lambda $$, $$\beta $$, and $$N$$ $$\mathit{\pm }$$ three standard deviations fall within their typical ranges of $$7\mathit{\le }{\lambda }_{t}\mathit{\le }13$$, $$0.7\mathit{\le }{\beta }_{t}\mathit{\le }1$$, and $$40\mathit{\le }N\mathit{\le }70$$ (De Michele & Avanzi, 2018; Marra et al., 2018, 2020; Miniussi & Marani, 2020) to resemble the statistical characteristics of the rainfall data set employed in the practical application. Figure 3 illustrates the ordinary-event distributions of simulations and observations and the extreme-event distributions estimated by Monte Carlo simulation (more details below). The baseline and prediction periods each cover 70 years to mimic data availability in general. For instance, the CMIP6 models consider the historical/baseline period as 1850–2014, and the future/prediction period up to the end of the century (2100).

Table 1 Distribution Specifications for Synthetic Data Generation in the Nonstationary Scenarios of Consistent and Changing Nonstationarity Patterns

Scenario	Hyperparameters for observations	Hyperparameters for simulations
Consistent nonstationarity	$${o}^{b}:\left\{\begin{array}{@{}cc@{}}{\lambda }_{0}=8,& {\lambda }_{1}=0.035;\\ {\beta }_{0}=1,& {\beta }_{1}=-0.002;\\ {\mu }_{{N}_{0}}=50;& \text{and}\quad {\sigma }_{{N}_{0}}=3\end{array}\right.$$	$${s}^{b}:\left\{\begin{array}{@{}c@{}}\begin{array}{cc}{\lambda }_{0}=7,& {\lambda }_{1}=0.025;\\ {\beta }_{0}=1.1,& {\beta }_{1}=-0.0015\,\\ {\mu }_{{N}_{0}}=60;& \text{and}\quad {\sigma }_{{N}_{0}}=3\end{array}\end{array}\right.$$$${s}^{p}\equiv {s}^{b}$$
Changing nonstationarity	$${o}^{b}:$$ same as above	$${s}^{b}:$$ same as above
$${s}^{p}:$$ same as above except for $${\lambda }_{2}=2{\lambda }_{1}$$ and $${\beta }_{2}=2{\beta }_{1}$$

[IMAGE OMITTED. SEE PDF]

In all synthetic data sets, both $${\lambda }_{j}$$ and $${\beta }_{j}$$ (estimated annually) display significant temporal trends according to the Mann-Kendall test (Kendall, 1975; Mann, 1945), aligning with the predefined nonstationary scenarios. Besides, the annual maximum series extracted from all data sets exhibit significant trends, confirming the presence of nonstationary signals in the extremes. Moreover, we confirm that the annual maximum series are suitably fitted by the nonstationary GEV distribution at $$\rho $$ by the two-sample Kolmogorov-Smirnov (2-sample-KS) test (Massey, 1951). Since this test requires homogeneity of the samples, we apply the standard transformation to the nonstationary annual maximum series using the corresponding reference distribution to allow the validation (see Coles, 2001; Ragno et al., 2019). We generate 1,000 pairs of data sets for each nonstationary scenario, resulting in a total of 4,000 data sets.

We assess the performance of all methods in both the baseline and prediction periods. The methods' evaluation encompasses distribution- and quantile-wise assessments, which examine the correction of the entire distribution and specific quantiles corresponding to certain cumulative probabilities $$H(Y)$$, respectively (see Section 3.3). In the quantile-wise assessment, the methods are assessed at the cumulative probabilities of 0.5, 0.9, and 0.99 (i.e., at the return periods of $$T$$ = 2, 10, and 100 years, respectively). Whereas in the distribution-wise assessment, the methods are evaluated based on 100 quantiles with equally spaced exceedance probability at each time slice. We obtain the actual quantiles from observations in each nonstationary scenario using Monte Carlo simulation by (a) generating $$1\times {10}^{7}$$ blocks, each with $${n}_{t}$$ ordinary events sampled from the $$\mathcal{W}(\cdot )$$ distribution corresponding to observations with $${\theta }_{t}$$ at a time $${t}_{k}$$, where $${\theta }_{t}$$ is pre-specified and $${n}_{t}$$ is randomly sampled from its latent distribution (Equations 9 and 10 and Table 1); (b) extracting the annual maxima from each block at $${t}_{k}$$, and accordingly yielding $$1\times {10}^{7}$$ realizations from the actual nonstationary extreme-event distribution at the time slice; and (c) obtaining the quantiles of interest from this extensive collection of block maxima using the unbiased Weibull plotting position formula. We repeated these steps for all time slices in the baseline and prediction periods.

Table 2 presents the abbreviations for all the bias correction methods employed in the simulation study. It is important to note that when the SMEV distribution is used in NS-QM and CNS-QM, only the actual (a priori known) nonstationary $$\mathcal{W}(\cdot )$$ distribution, that is, $$SME{V}_{1,1,0}$$ is adopted, as N is stationary. Recall that $$GE{V}_{1,0,0}$$ and $$GE{V}_{1,1,0}$$ are statistically significantly appropriate for describing the annual maximum series extracted from the synthetic daily series. Hence, both $$GE{V}_{1,0,0}$$ and $$GE{V}_{1,1,0}$$ were adopted.

Table 2 Summary Description of Bias Correction Methods and Stationary and Nonstationary Probability Distributions

Class	Bias correction method	Distribution	Nonstationary structure	Abbreviation
Stationary	QM	GEV	–	QM-GEV
SMEV	–	QM-SMEV
Quasi-stationary	DQM	GEV	–	DQM-GEV
SMEV	–	DQM-SMEV
QDM	GEV	–	QDM-GEV
SMEV	–	QDM-SMEV
CDFt	GEV	–	CDFt-GEV
SMEV	–	CDFt-SMEV
Nonstationary (consistent)	CNS-QM	GEV	1,0,0 ($${\xi }_{t}$$)	CNS-QM-GEV1,0,0
1,1,0 ($${\xi }_{t}$$ and $${\alpha }_{t}$$)	CNS-QM-GEV1,1,0
SMEV	1,1,0 ($${\lambda }_{t}$$ and $${\beta }_{t}$$)	CNS-QM-SMEV
Nonstationary (changing)	NS-QM	GEV	1,0,0 ($${\xi }_{t}$$)	NS-QM-GEV1,0,0
1,1,0 ($${\xi }_{t}$$ and $${\alpha }_{t}$$)	NS-QM-GEV1,1,0
SMEV	1,1,0 ($${\lambda }_{t}$$ and $${\beta }_{t}$$)	NS-QM-SMEV

Practical Applications

We demonstrate a detailed application of the proposed and benchmark methods using daily precipitation observations from a station and the outputs of a regional climate model at the corresponding grid cell as simulations (Figure 4). The observational data set was recorded at the hydrometeorological station of Cuitzeo in Mexico (MXN00016027 in the Global Historical Climatology Network daily (GHCNd) database (Menne et al., 2012)), where snowfall does not occur. This station has year-round data for 84 years from 1923 to 2017. The simulations are the outputs of the WRF-MPI-ESM-LR model, which is the high-resolution Weather Research and Forecasting regional climate model (Skamarock et al., 2008) forced using boundary conditions by the Max-Planck-Institute Earth System Model (global climate model) with low-resolution setup (MPI-ESM-LR) (Giorgetta et al., 2013) from the CMIP5 archive. This climate model was run over North America at a spatial resolution of 0.22° (∼25 km) for the historical period of 1950–2005 and the future period of 2006–2100 using the RCP 8.5 emissions scenario. The regional climate model output data were downloaded from the North American Coordinated Regional Climate Downscaling Experiment (NA-CORDEX) archive (Mearns et al., 2017).

[IMAGE OMITTED. SEE PDF]

We conducted exploratory data analysis assisted by deductive reasoning (Serinaldi et al., 2018; Vidrio-Sahagún et al., 2024) to assess the stationarity assumption at the selected site. The detected significant temporal trend in the mean of observed annual maxima precipitation suggests nonstationarity in the extremes (Figure 4b). Similarly, the regional climate model outputs show nonstationarity in the baseline period but underestimate it compared to observations; while in the future horizon, an insignificant and weaker trend (compared to the baseline period) was detected in the climate model projections. No significant trends were detected in the variability of annual maxima precipitation when using the moving-windows Mann-Kendall test. Both observations and simulations in the baseline period show the absence of significant trends in the $$n$$, but the climate model simulations suggest a significant decreasing trend in the future. Furthermore, the annual maxima of maximum daily temperature, which is often related to extreme precipitation, presents significant trends in both observations and simulations during the historical period; while an increase in the trend in the future was projected by the climate model. The detected consistent increases in annual maximum precipitation and temperature align with the notion that a warmer atmosphere is expected to increase its moisture-holding capacity and consequently lead to more extreme precipitation (Kharin et al., 2013; Trenberth, 2011). Thus, adopting nonstationary correction methods is considered appropriate for this study case.

We also applied the methods to other sites in different regions (Table S1 in Supporting Information S1). We used observations from the GHCNd database (Menne et al., 2012) and CMIP6 (Eyring et al., 2016) climate simulations from the Max-Planck-Institute Earth System Model MPI-ESM1-2-HR (Mauritsen et al., 2019). This global climate model with a ∼100 km nominal resolution was run for the historical period of 1850–2014 and future period of 2015–2100 under the SSP585 scenario. These additional sites in Australia, China, Japan, and the USA (Texas, Mississippi, and Georgia) have more than 50 years of observations and over 80% year-round data completeness. The climate model qualitatively captured local nonstationarity signatures in extremes at these sites, that is, trend directions, meeting the key requirement for meaningful nonstationary correction. These sites also exhibit diverse sample sizes, statistical characteristics, bias levels in extreme precipitation magnitudes and trend magnitudes (underestimation in climate simulations), and nonstationarity patterns in the ordinary event distributions (Tables S2–S3, and Figure S1 in Supporting Information S1).

Given that the distribution of extreme precipitation is virtually unknown, especially for the prediction period, we demonstrate the practical application of the methods in two ways. On one hand, we perform bias correction on the time series of climate model simulations, as is the common practice. The resulting corrected time series serves as the post-processed product used as input for impact models (e.g., hydrological models) to investigate climate change impacts. On the other hand, to delve into the distinctions among the methods, we correct three representative quantiles of the raw simulations associated with cumulative probabilities $$H(Y)=0.5,0.9,$$ and $$0.99$$, which are treated as either constant or time-varying, depending on the stationarity/nonstationarity treatment of the underlying method.

In the NS-QM and CNS-QM methods, we select the nonstationary structure of the distribution based on the identified nonstationary signals. This approach has been shown more effective in capturing nonstationary stochastic processes than the traditional selection based on performance metrics in the decomposition-based nonstationary frequency analysis (Vidrio-Sahagún & He, 2022b). In the nonstationary distributions, we used time as the covariate (i.e., $${v}_{t}=t$$), reflecting the nonstationary signals (trends) in the climate model projections of precipitation extremes. Attributing nonstationarity to physical covariates (Slater et al., 2021), though valuable, can be challenging and is beyond this paper's scope.

Performance Assessment Metrics

The distribution-wise assessment evaluates the similarity between the bias-corrected and actual distributions in two ways: (a) employing the 2-sample-KS test (Massey, 1951) to evaluate whether $${H}_{{\widehat{Y}}_{t}^{s,p}}(\cdot )$$ (when using the GEV) or $${{\Psi }}_{{\widehat{X}}_{t}^{s,p}}(\cdot )$$ (when using the SMEV) and $${H}_{{Y}_{t}^{o,p}}(\cdot )$$ follow the same distribution at the $$\rho $$ significance level at each time slice, and reporting the rejection rate ($$RR\mathit{\in }[0,1]$$); and (b) employing the overall Skill Score ($$SS\mathit{\in }(-\infty ,1]$$), which provides a global measure of accuracy and accounts for correlation and conditional and unconditional biases. The SS and RR metrics are given by: 11 $\begin{align*}SS\left({\widehat{y}}_{t,H(\cdot )},{y}_{t,H(\cdot )}\right)={\rho }_{{\widehat{y}}_{t,H(\cdot )},{y}_{t,H(\cdot )}}^{2}-{\left[{\rho }_{{\widehat{y}}_{t,H(\cdot )},{y}_{t,H(\cdot )}}-\frac{{\sigma }_{{\widehat{y}}_{t,H(\cdot )}}}{{\sigma }_{{y}_{t,H(\cdot )}}}\right]}^{2}-{\left[\frac{{\mu }_{{\widehat{y}}_{t,H(\cdot )}}-{\mu }_{{y}_{t,H(\cdot )}}}{{\sigma }_{{y}_{t,H(\cdot )}}}\right]}^{2}\end{align*}$ 12 $\begin{align*}RR=\frac{1}{M}\sum \nolimits_{t=1}^{M}{c}_{t},\text{where}\,{c}_{k}=\left\{\begin{array}{@{}c@{}}1:\text{if}\,{H}_{0}\,\text{at}\,t\,\text{is}\,\text{rejected}\,\\ 0:\text{otherwise}\end{array}\right.\end{align*}$ where $${y}_{t,H(\cdot )}$$ and $${\widehat{y}}_{t,H(\cdot )}$$ are the actual and bias-corrected quantiles associated with the effective (nonstationary) cumulative probabilities $$H(\cdot )$$ at time $$t$$, respectively; $${\rho }_{\widehat{y},y}$$ is the correlation between $${y}_{t,H(\cdot )}$$ and $${\widehat{y}}_{t,H(\cdot )}$$; $${\sigma }_{{y}_{t,H(\cdot )}}$$ and $${\sigma }_{{\widehat{y}}_{t,H(\cdot )}}$$ are the conditional (nonstationary) standard deviations of $${y}_{t,H(\cdot )}$$ and $${\widehat{y}}_{t,H(\cdot )}$$, respectively; $${\mu }_{{\widehat{y}}_{t,H(\cdot )}}$$ and $${\mu }_{{y}_{t,H(\cdot )}}$$ are the conditional mean of $${y}_{t,H(\cdot )}$$ and $${\widehat{y}}_{t,H(\cdot )}$$, respectively; and the null hypothesis $${H}_{0}$$ of the 2-sample-KS test is that both samples come from the same distribution. In the simulation study, the SS was calculated considering quantiles with equally spaced cumulative probabilities across the distribution at each time slice.

In the quantile-wise assessment, we evaluate the accuracy of the estimates using the relative root mean square error ($$RRMSE\mathit{\in }[0,\infty )$$) and relative bias ($$RBias\mathit{\in }[-\infty ,\infty )$$), and quantified the uncertainty of the estimates using the relative average bandwidth (RAW) and the relative coverage width index (RCWI) of the uncertainty bands at certain $$H(\cdot )$$ s. These metrics are given by: 13 $\begin{align*}RRMSE=\sqrt{\frac{1}{M}\sum \nolimits_{t=1}^{M}{\left(\frac{{\widehat{y}}_{t,H(\cdot )}-{y}_{t,H(\cdot )}}{{y}_{t,H(\cdot )}}\right)}^{2}}\end{align*}$ 14 $\begin{align*}RBias=\frac{1}{M}\sum \nolimits_{t=1}^{M}\frac{{\widehat{y}}_{t,H(\cdot )}-{y}_{t,H(\cdot )}}{{y}_{t,H(\cdot )}}\end{align*}$ 15 $\begin{align*}RAW=\frac{1}{M}\sum \nolimits_{t=1}^{M}\left(\frac{{C}_{t,H(\cdot )}^{U}-{C}_{t,H(\cdot )}^{L}}{{\widehat{y}}_{t,H(\cdot )}}\right)\end{align*}$ 16 $\begin{align*}CF=\frac{1}{M}\sum \nolimits_{t=1}^{M}{C}_{t,H(\cdot )},\text{where}\,{C}_{t,H(\cdot )}=\left\{\begin{array}{@{}c@{}}1\quad \forall \,t\ s.t.\ {C}_{t,H(\cdot )}^{L}\le {y}_{t,H(\cdot )}\le {C}_{t,H(\cdot )}^{U}\\ 0\quad \text{elsewhere}\end{array}\right.\end{align*}$ 17 $\begin{align*}RCWI=RAW\cdot \mathrm{exp}{([1-\rho ]-CF)}^{2}\end{align*}$ where $${C}_{t,H(\cdot )}^{U}$$ and $${C}_{t.H(\cdot )}^{L}$$ are the upper and lower uncertainty bounds at $$t$$ and $$H(\cdot )$$; $$CF$$ is the coverage frequency; and M is the sample size (years).

Results and Discussion

Simulation Study Results

Quantile-Wise Performance Assessment

Figure 5 shows the RRMSE and RBias of bias-corrected quantiles at $$H(Y)=0.5,0.9,$$ and $$0.99$$ during both the baseline and prediction periods and those of uncorrected quantiles under the scenarios of consistent and changing nonstationarity. In general, all the methods reduced the errors and biases of raw simulations. In the baseline period of the scenario of consistent nonstationarity, QM-SMEV, NS-QM-SMEV, and CNS-QM-SMEV exhibited similar medians of RRMSE and outperformed their GEV-based counterparts (QM-GEV, NS-QM-GEV_1,0,0, NS-QM-GEV_1,1,0, CNS-QM-GEV_1,0,0 and CNS-QM-GEV_1,1,0). All methods except CNS-QM-GEV_1,0,0 offered unbiased estimates, as their RBias medians were close to zero. The widths of interquartile ranges of SMEV-based methods were narrower than those of GEV-based methods, more notably for higher $$H(Y)$$ s. In the prediction period of this scenario, the methods offered more distinguishable performance. The NS-QM-GEV method showed advantages over stationary and quasi-stationary methods with the same distribution (i.e., QM-GEV, DQM-GEV, QDM-GEV, and CDFt-GEV), as it yielded a smaller RBias median but a comparable RRMSE median, and wider interquartile ranges of both metrics. The CNS-QM-GEV method yielded overall comparable or slightly smaller RBias and RRMSE medians and similar or narrower interquartile ranges compared to the stationary and quasi-stationary GEV-based methods, especially when using $$GE{V}_{1,1,0}$$. In addition, CNS-QM-GEV_1,1,0 produced a lower RRMSE median and narrower interquartile ranges of RRMSE and RBias than NS-QM-GEV in general, although with a larger median of RBias. Among all the methods, NS-QM-SMEV and CNS-QM-SMEV always yielded the best overall performance, as they yielded the lowest medians of RRMSE, and their RBias medians were close to zero across all $$H(Y)$$ s. All these demonstrate that NS-QM-SMEV and CNS-QM-SMEV were overall superior to all other methods in both the baseline and prediction periods of the scenario of consistent nonstationarity. However, the stationary and quasi-stationary SMEV-based methods always produced slightly narrower interquartile ranges in RRMSE and RBias.

[IMAGE OMITTED. SEE PDF]

In the baseline period of the scenario of changing nonstationarity, CNS-QM-GEV produced larger RRMSE and RBias (negative) medians among the GEV-based methods, especially when using $$GE{V}_{1,0,0}$$, indicating underestimation. The QM-GEV and NS-QM-GEV methods generally performed equivalently with comparable RRMSE medians and unbiasedly with small RBias medians. The interquartile ranges of all GEV-based methods were comparable in general (except CNS-QM-GEV_1,0,0, especially at high $$H(Y)$$). In addition, all SMEV-based methods were roughly unbiased and yielded similar interquartile ranges of RBias, albeit CNS-QM-SMEV produced a larger error (i.e., higher RRMSE median and wider interquartile range) than the other SMEV-based methods across all $$H(Y)$$ s. For both GEV- and SMEV-based methods, the outperformance of NS-QM over CNS-QM is not surprising, as CNS-QM combines data from both the baseline and prediction periods for parametrizations assuming consistent nonstationarity, which degrades the nonstationarity description in the baseline period of this scenario. The SMEV-based methods consistently outperformed the GEV-based, especially by their lower RRMSE medians, narrower interquartile ranges of RRMSE, and narrower interquartile ranges of RBias. On the other hand, in the prediction period of this scenario, CNS-QM-GEV was more clearly inferior to NS-QM-GEV, as it produced larger RRMSE and RBias due to its assumption of consistent nonstationarity across the periods. NS-QM-GEV was not obviously better than stationary/quasi-stationary GEV-based methods, as they had similar or higher RRMSE, albeit being roughly unbiased. This result might be ascribed to using more parameters and distributions in NS-QM-GEV, leading to parametrizations subject to more pronounced errors. NS-QM-SMEV again demonstrated superiority, as it yielded the lowest RRMSE median and the RBias median closest to zero among all the methods. All the results from both scenarios demonstrate that NS-QM-SMEV outperformed all other methods and thus argue that both correct capturing of the nonstationarity and incorporating ordinary-event records improve the accuracy of bias correction.

Figure 6 presents the uncertainty evaluated in terms of RAW and RCWI in bias-corrected quantiles at $$H(Y)=0.5,0.9,$$ and $$0.99$$. Overall, the methods yielded similar results in both scenarios of nonstationarity. In the baseline periods, all the SMEV-based methods exhibited lower uncertainty than the GEV-based methods, as evidenced by their lower RAW and RCWI medians as well as their narrower interquartile ranges. The overall lower uncertainty of SMEV-based methods could be attributed to improved parametrization owing to using more extensive data sets (i.e., of ordinary events). It is worth noting that NS-QM-SMEV and CNS-QM-SMEV showed slightly higher uncertainty than QM-SMEV. NS-QM-GEV_1,1,0 and CNS-QM-GEV_1,1,0 presented the highest RAW and RCWI medians among all methods. In the prediction period, SMEV-based methods were also less uncertain than GEV-based methods. NS-QM-SMEV and CNS-QM-SMEV produced slightly higher RAW medians and overall comparable RCWI medians compared to the stationary and quasi-stationary SMEV-based methods. Among all methods, NS-QM-GEV showed the highest uncertainty, followed by CDFt-GEV. The increased uncertainty in NS-QM-GEV could be attributed to its use of more distributions (three), more complex model structures (with four or five parameters), and/or sole use of annual maxima data compared to other methods, making its out-of-sample estimates more uncertain. CDFt-GEV also uses three distributions directly (unlike DQM and QDM, see Appendix), each with three parameters, and relies solely on annual maxima data; this results in higher uncertainty of CDFt-GEV than several other methods that use fewer distributions or more data. However, the superiority of SMEV-based methods over GEV-based methods in RCWI was not as prominent as in RAW due to their wider interquartile range of RCWI. Given that RCWI considers the coverage of the actual quantiles by the uncertainty band, caution is advised in interpreting narrow uncertainty bands of SMEV-based methods, as they might be overly narrowed. Therefore, NS-QM-SMEV and CNS-QM-SMEV enhanced bias-correction of extremes by substantially lowering uncertainty compared to their nonstationary GEV-based counterparts and were competitive with stationary/quasi-stationary SMEV-based methods.

[IMAGE OMITTED. SEE PDF]

Distribution-Wise Performance Assessment

Figure 7 presents the distribution-wise assessment of bias-corrected quantiles in terms of RR and SS for all methods under both scenarios, along with the RR and SS of raw simulations. Overall, all methods improved the RR and SS of raw simulations. In the baseline period of the scenario of consistent nonstationarity, the NS-QM-GEV and CNS-QM-GEV methods presented higher RR medians and upper RR quartiles than all other methods, especially CNS-QM-GEV_1,0,0. In terms of SS, CNS-QM-GEV_1,0,0 was inferior to all the other GEV-based methods, which preformed equivalently. All these indicate the low performance of NS-QM-GEV and CNS-QM-GEV to match the actual distribution of the extremes. In contrast, all the SMEV-based methods performed approximately optimally, as they had RR medians close to zero, SS medians close to one, and negligible interquartile ranges of both RR and SS. In the prediction period of this scenario, NS-QM-GEV and CNS-QM-GEV were not always preferable compared to their stationary and quasi-stationary GEV-based counterparts, as they produced SS and RR that were equivalent or worse. In contrast, the NS-QM-SMEV and CNS-QM-SMEV methods maintained the lowest RR median and highest SS median among all methods and were superior to all the other SMEV- and GEV-based methods in capturing the entire actual distribution of extremes.

[IMAGE OMITTED. SEE PDF]

In the scenario of changing nonstationarity in the baseline period, the methods gave similar results as in the scenario of consistent nonstationarity, except that besides CNS-QM-GEV_1,0,0, CNS-QM-GEV_1,1,0 was also not favorable in both RR and SS. The SMEV-based methods (except CNS-QM-SMEV) outperformed all the GEV-based methods. In the prediction period of this scenario, NS-QM-SMEV offered the best performance, as it yielded the lowest RR median and highest SS median with a narrower or roughly equally wide interquartile range than all other methods. Overall, the CNS-QM method, irrespective of the distribution (either GEV or SMEV), did not show an advantage over its counterparts using the same distribution. This argues the importance of duly capturing nonstationarity when bias-correcting simulations. The superior performance of the NS-QM-SMEV to reproduce the entire distribution of extremes supports that correctly modeling nonstationarity and incorporating ordinary-event data can improve bias correction under nonstationarity.

The Quantile-Quantile Matching Mechanism

All these methods should work by translating every simulation data point $${y}_{t}^{s}$$ from the distribution of simulations to the distribution of observations while maintaining an identical cumulative probability $$H(Y)$$. For example, if $${y}_{t}^{s}$$ is associated with the cumulative probability $${H}_{{Y}_{t}^{s}}(\cdot )$$, its bias-corrected estimate should be the quantile $${H}_{{Y}_{\mathrm{t}}^{o}}^{-1}(\cdot )$$ associated with the same probability. Thus, we explored the faithfulness of such quantile-quantile matching under nonstationarity for all the methods. Figure 8 presents the cumulative probabilities expressed as return periods ($$T=1/\left[1-{H}_{{Y}_{t}}(\cdot )\right]$$) used in the bias correction in the prediction periods of both scenarios and compares them with the actual $$T$$ with which they should have been corrected. An ideal bias correction would use $$T$$s in the procedure identical to the actual $$T$$s, corresponding to the 1:1 diagonal. As shown in Figure 8a, all the stationary and quasi-stationary GEV-based methods tended to correct the simulations typically using a larger $$T$$ than the actual $$T$$, especially for the most extreme events. The most critical case is the QM-GEV method, which often overestimated the $$T$$ in the procedure; for instance, the $$T=$$ 100-year quantile was often corrected as the $$T=$$ 10,000-year quantile. Figure 1 provides a visual illustration of the origin of this quantified mismatch. NS-QM-GEV and CNS-QM-GEV (Figure 8b) improved the quantile-quantile matching compared to QM-GEV but not compared to DQM-GEV, QDM-GEV, and CDFt-GEV methods, as they still produced biases toward higher $$T$$ s and wide variability envelops around the 1:1 diagonal. Stationary and quasi-stationary SMEV-based methods (Figure 8c) produced similar results to their GEV-based counterparts, but their variability envelope was reduced, indicating more robustness. In contrast, NS-QM-SMEV and CNS-QM-SMEV were roughly unbiased, as their variability envelopes were distributed around the 1:1 diagonal approximately evenly, indicating their superiority to all the other methods. Besides, these two methods presented narrower variability, showcasing their robust behavior. These results justify the better performance of NS-QM-SMEV and CNS-QM-SMEV under nonstationarity, which results from minimizing the quantile-quantile mismatch and consequently increasing the reliability in bias-correcting future extremes.

[IMAGE OMITTED. SEE PDF]

Effects of Sample Size in Bias Correction Under Nonstationarity

Apart from the bias correction methods themselves, the sample size of data sets is another relevant factor affecting their performance. As shown previously (Figures 5–8), CDFt-GEV, NS-QM-GEV_1,1,0 and CDFt-SMEV overall outperformed or performed similarly to other quasi-stationary GEV-, nonstationary GEV-, and quasi-stationary SMEV-based methods, respectively. In this section, we utilize NS-QM-SMEV and CNS-QM-SMEV alongside these three methods and their stationary counterparts to illustrate the impact of sample size on bias correction in quantile- and distribution-wise assessments.

Figure 9 shows the impact of sample size (ranging from 50 to 100 years in each period) on the RRMSE and RBias medians of the selected methods in the quantile-wise assessment. In the baseline periods of both scenarios, all methods exhibited comparable performance, generally reducing RRMSE with increasing sample size (Figures 9a–9c) and keeping their RBias close to zero (Figures 9d–9f). The exceptions were CNS-QM-SMEV in RRMSE of the scenario of changing nonstationarity and QM-GEV in RBias of both scenarios at high $$H(Y)$$. The overall improvement in accuracy with larger sample sizes suggests an increase in goodness-of-fit, even with statistical misrepresentations of the methods (e.g., excluding nonstationary treatment or neglecting the changing nonstationarity between periods). However, the performance of the methods markedly diverged in the prediction period, revealing the mismatch between their statistical descriptions and the actual nonstationary processes. In both scenarios, stationary and quasi-stationary methods, regardless of the distribution, displayed increasing or approximately constant RRMSE and negative RBias as the sample size augmented. In contrast, both NS-QM and CNS-QM, irrespective of the distribution used, consistently reduced or maintained RRMSE and RBias with increasing sample size when aligned with the nonstationary scenario (i.e., NS-QM under both scenarios and CNS-QM under consistent nonstationarity). Under consistent nonstationarity, NS-QM-SMEV and CNS-QM-SMEV consistently offered smaller RRMSE and RBias than all other methods for sample sizes of around 70 years, and their superiority was more prominent as the sample size increased. Similarly, under changing nonstationarity, NS-QM-SMEV outperformed the other methods from sample sizes of 70 years onward. This finding supports that capturing nonstationarity correctly becomes more critical for bias correction over longer periods.

[IMAGE OMITTED. SEE PDF]

Figure 10 further displays the impact of sample size on the SS and RR medians of the selected bias correction methods in the distribution-wise assessment. Generally, the SMEV-based methods performed better than the GEV-based methods in terms of both SS and RR across most sample sizes in both periods under both nonstationary scenarios. In the baseline period, SMEV-based methods were insensitive to the sample size, as SS and RR were roughly constant over the range of sample sizes, while the GEV-based methods improved their SS and RR with increasing sample sizes. However, in the prediction period, similar to the quantile-wise assessment, the performance of the methods diverged more noticeably, especially for larger sample sizes. In both scenarios, the performance of QM and CDFt, independently of the distribution, degraded with increasing sample size, as shown by their increasing RR; the performance of QM also deteriorated in terms of SS. Besides, in the scenario of changing nonstationarity, CNS-QM-SMEV also degraded in SS and RR with longer data sets. The degrading performance of QM, CDFt, and CNS-QM with increasing sample sizes shows that incorrect nonstationarity characterization progressively distorts the statistical characteristics of the entire distribution in longer prediction periods. In contrast, NS-QM and CNS-QM consistently increased SS, decreased RR, or kept optimal performance with increasing sample sizes in the scenario of consistent nonstationarity. The outperformance of NS-QM-SMEV and CNS-QM-SMEV started around a sample size of 70 years and became more evident toward larger sample sizes. Similarly, NS-QM improved or kept optimal SS and RR with increasing sample sizes in the scenario of changing nonstationarity. NS-QM-SMEV was superior to all other methods from 70 years onward in this scenario. These results reaffirm the superiority of NS-QM-SMEV and CNS-QM-SMEV, especially when the sample size is approximately 70 years or more. Thus, the data availability could constraint the implementation of nonstationary bias correction in practice.

[IMAGE OMITTED. SEE PDF]

Practical Application Results

This section demonstrates the bias correction methods for adjusting climate model projections of precipitation extremes through two applications. Application 1 provides a detailed example at a site in Mexico (MXN00016027) using the outputs of a high-resolution regional climate model driven by a coarse-scale global climate model. Application 2 summarizes results from six additional sites and shows the methods' applicability in other regions using the outputs of a coarse-scale global climate model.

Application 1: A Detailed Real-World Example

The left-censoring thresholds for Cuitzeo (MXN00016027) were set to 50% and 70% for observations and simulations, respectively, based on the AIC optimization and proper fit of the $$\mathcal{W}(\cdot )$$ distribution of ordinary events. In the baseline period, annual estimates of $$\lambda $$ of the $$\mathcal{W}(\cdot )$$ distribution for both observations and simulations show significant increasing trends (results not shown), aligning with the trends in rainfall annual maxima (Figure 4b). In the prediction period, regional climate model projections suggest a weakening trend in $$\lambda $$, similar to the projected annual maxima in this period (Figure 4b). Yet, when assuming a constant nonstationarity pattern (as in the CNS-QM method) and assessing all climate model simulations from 1950 to 2100 collectively, the upward trend in $$\lambda $$ remains statistically significant. Given the presence of a significant trend only in the mean of annual maxima, we opted for the $$GE{V}_{1,0,0}$$ in the nonstationary GEV-based method. Similarly, due to the trends presented in $$\lambda $$ and $$n$$ (Figure 4c), we adopted the $$SME{V}_{1,0,1}$$ (see Equation 8) in the nonstationary SMEV-based methods. Apart from the nonstationary GEV- and SMEV-based methods, one method of each class (i.e., stationary and quasi-stationary GEV- and SMEV-based methods) was also included for comparison purposes. These benchmark methods are QM-GEV, CDFt-GEV, QM-SMEV, and CDFt-SMEV, which were shown to generally outperform or perform comparably to their counterparts in the simulation study.

Bias-Correction of Time Series Projections

It is worth noting that, at the study site, the highest daily rainfall on record was 167 mm/d in 2015, which almost doubled the second-highest daily rainfall of 90 mm/d in 2004 (Figure 4b). Besides, the regional climate model slightly overestimated extreme precipitation in the baseline period, as it yielded an annual maxima rainfall median of 46 mm/d, while the observations' median was 39 mm/d. Figure 11 presents the time series $${y}_{t}^{s,p}$$ and $${\widehat{y}}_{t}^{s,p}$$ along with the uncertainty yielded by all the selected methods in the prediction period (2006–2010). The methods overall mitigated the overestimation of the climate model by decreasing the magnitudes of $${y}_{t}^{s,p}$$ in general. Yet, GEV-based methods yielded notably high estimates for the most extreme events. For instance, both QM-GEV and NS-QM-GEV produced a few (i.e., six) extreme events exceeding 167 mm/d, with uncertainty up to $$2.3\times {10}^{4}$$ and $$4.5\times {10}^{6}$$ mm/d, respectively. CDFt-GEV produced lower estimates for the most extreme events but still yielded four extreme events exceeding 167 mm/d and with uncertainty up to $$1.5\times {10}^{6}$$ mm/d. The very high estimates of these methods could result from their proneness to quantile-quantile mismatch (Figure 8). In contrast, CNS-QM-GEV avoided the very high estimates, with a maximum estimate of 129 mm/d and uncertainty up to 274 mm/d. As aforementioned, CDFt-GEV and NS-QM-GEV use more distributions (three) than their GEV-based counterparts, and NS-QM-GEV uses more complex distributions (five distribution parameters). Conversely, CNS-QM-GEV produced more robust estimates, which could be ascribed to using fewer distributions (two) and combining the baseline and prediction periods to enhance fitting. On the other hand, all the SMEV-based methods avoided yielding exceptionally high estimates and produced notably narrower confidence intervals. Both QM-SMEV and CDFt-SMEV produced estimates close to $${y}_{t}^{s,p}$$, but with roughly constant adjustments irrespective of time. In contrast, both CNS-QM-SMEV and NS-QM-SMEV displayed time-varying bias adjustments and increased the nonstationarity degree of $${\widehat{y}}_{t}^{s,p}$$ compared to that of $${y}_{t}^{s,p}$$. For instance, $${y}_{t}^{s,p}$$ had an insignificant trend of 0.6 mm/d/decade (Figure 4b), while $${\widehat{y}}_{t}^{s,p}$$ by CNS-QM-SMEV displayed a significant trend of 2.1 mm/d/decade. NS-QM-SMEV produced a milder nonstationarity degree than CNS-QM-SMEV (an insignificant trend of 0.9 mm/d^/decade). This feature of the NS-QM and CNS-QM methods is essential for correcting potential underestimation of nonstationarity in the extremes simulated by climate or forecast models. These results highlight that NS-QM-SMEV and CNS-QM-SMEV can correct biases more realistically, including bias correction in nonstationarity itself, and reduce estimation uncertainty and errors.

[IMAGE OMITTED. SEE PDF]

Bias Correction of Simulated Quantiles

Figure 12 presents the bias-corrected quantiles for $$H(Y)=0.5,0.9$$ and $$0.99$$ for $${y}_{t}^{s,b}$$ and $${y}_{t}^{s,p}$$ along with their uncertainties and observations. These corrected quantile estimates are of paramount importance in infrastructure design, hazard assessment, and risk mitigation/adaptation. The QM-GEV yielded constant estimates over time due to the use of time-invariant distributions, and more pronounced increases in uncertainty as $$H(Y)$$ increases. CDFt-GEV considers relative changes between the baseline and prediction periods and thus produced sudden changes in quantiles in 2005. Especially, the sudden change in quantile for $$H(Y)=0.99$$ appears unrealistic. This raises questions about the realism and utility of the stationary and quasi-stationary methods, especially when estimating rare events for the short-term future horizon. In contrast, CNS-QM-GEV provided more realistic estimates showing a continuous evolution over time and having narrower uncertainty bands than QM-GEV and CDFt-GEV. NS-QM-GEV showed the temporal evolution of quantiles, but like CDFt-GEV (which also uses nested distributions parametrized independently in each period), it yielded a large shift in quantiles at the transition of the periods too, particularly for $$H(Y)=0.99$$, and large uncertainties. In contrast, it is not surprising that all the SMEV-based methods showed very low uncertainty compared to the GEV-based methods. QM-SMEV substantially reduced uncertainty compared to QM-GEV and maintained constant estimates with a lower magnitude. Similar to CDFt-GEV and NS-QM-GEV, CDFt-SMEV and NS-QM-SMEV also yielded sudden shifts in quantiles at the transition of the periods, but they were small (e.g., changes of 2 and 3 mm/d for $$H(Y)=0.99$$ for CDFt-SMEV and NS-QM-SMEV, respectively). In the study case, the small shifts are attributed to the relative changes in annual maxima balanced with the changes in the statistical characteristics of ordinary events. Like CNS-QM-GEV, the quantiles produced by CNS-QM-SMEV continuously evolved over time, reflecting nonstationarity impacts more realistically and avoiding shifts at the transition of the periods. Compared to CNS-QM-SMEV, NS-QM-SMEV yielded a less pronounced trend in quantiles in the prediction period, aligning with the weaker trend simulated by the climate model (Figure 4b). This practical application illustrates the feasibility and utility of NS-QM-SMEV and CNS-QM-SMEV for bias correction, as they produced quantiles that coherently vary over time and are less uncertain compared to other methods.

[IMAGE OMITTED. SEE PDF]

Application 2: Additional Examples From Diverse Regions

The results of the study cases in different regions using CMIP6 outputs (Figures S2–S4 in Supporting Information S1 for time series correction and Figures S5–S7 in Supporting Information S1 for quantile correction in Supplementary materials) are consistent with those at MXN00016027. Specifically, in the time series bias correction, all methods generally reduced the climate model's bias (underestimation). GEV-based methods tended to highly overestimate the most extreme events, while CNS-QM-GEV generally mitigated this issue. SMEV-based methods consistently provided more reasonable estimates with narrower confidence intervals. CNS- and NS-QM-SMEV offered time-varying adjustments, often increasing the nonstationarity degree in corrected estimates compared to raw projections. In quantile correction, QM-GEV produced constant estimates, while CDFt-GEV exhibited abrupt changes that were often unrealistic for short-term predictions. CNS-QM-GEV yielded more realistic, gradually evolving estimates. NS-QM-GEV captured temporal changes but, like CDFt-GEV, often resulted in large shifts at period transitions, even when nonstationarity seems approximately consistent (at ASN00017014, JA000047898, and USC00226084). SMEV-based methods showed lower uncertainty than GEV-based methods. While CDFt-SMEV and NS-QM-SMEV exhibited shifts at period transitions, these were smaller than those in CDFt-GEV and NS-QM-GEV. CNS-QM-SMEV produced continuously evolving quantiles that realistically reflect nonstationarity and avoided abrupt transitions but with less flexibility to capture changes in nonstationarity. Overall, these results reaffirm that NS- and CNS-QM-SMEV can correct biases more effectively while reducing uncertainty and gross errors, like substantial overestimation and large discontinuities. However, we observed numerical issues in point and confidence interval estimates at two sites (CHM00057237 and JA000047898), which had the shortest samples, when using methods based on the CDFt principle (CDFt and nonstationary NS-QM approaches).

Limitations and Future Research

In both the simulation study and practical applications, the proposed NS- and CNS-QM-SMEV demonstrated advantages in bias correction for nonstationary extreme events. However, there are opportunities for improvement and further research. For instance, methods based on the CDFt principle (both existing quasi-stationary CDFt and proposed nonstationary NS-QM approaches) can encounter numerical issues in point and confidence interval estimates, as observed at two sites in Application 2. As noted by Lanzante et al. (2019, 2020), these issues in the CDFt principle may arise from multiple conversions between probabilities and quantiles across the four involved distributions, particularly near their tails, and the short sample sizes. Further numerical improvements are needed to address this in practical applications.

Moreover, the NS- and CNS-QM methods assume a constant (stationary) transfer function between observations and simulations. However, addressing nonstationary biases might require time-varying relationships between observations and simulations (Van De Velde et al., 2022). While nonstationary transfer functions could improve bias correction, they introduce additional challenges. For example, determining their specific structural form and suitable covariates adds complexity and may increase the risk of overfitting. Besides, since these functions can only be calibrated with historical data, their future validity is uncertain. Future research should synthesize the advantages and disadvantages of using nonstationary transfer functions. Furthermore, the developed nonstationary methods used distributional transfer functions, but other transfer functions, such as empirical transformations (e.g., based on linear/power regression or smoothing splines), could be considered.

The diagnosis and characterization of nonstationarity, for both the distributions and bias correction method (nonstationary formulation choice), are crucial aspects of implementing NS- and CNS-QM confidently in real-world applications. This can be conducted using deductive reasoning alongside exploratory data analysis (as demonstrated in this paper for the distributions) or through more detailed nonstationarity attribution approaches (Koutsoyiannis & Montanari, 2015; Merz et al., 2012; Serinaldi et al., 2018; Slater et al., 2021). However, there is a trade-off between refining nonstationarity characterization and increasing estimation uncertainty, especially with limited data availability. This paper focused on parsimonious nonstationary structures, but more complex structures (e.g., non-linear functions of covariates or adopting multiple covariates) can be explored, though they would also raise estimation uncertainty (Ouarda et al., 2019; Serinaldi & Kilsby, 2015). Developing suitable approaches and guidelines to identify proper nonstationary structures for the distributions and nonstationary schemes for the bias correction methods would be highly beneficial. Future research could also explore the potential of incorporating more complex nonstationary structures in the SMEV-based methods while maintaining practical applicability.

In addition, the proposed methods used deterministic covariates but can be extended to incorporate nonstationary stochastic physical covariates, which may require detailed nonstationarity attribution (Slater et al., 2021). Yet, nonstationary stochastic physical covariates often involve covariate-dependent characterizations (el Adlouni et al., 2007; Ouarda et al., 2018) that result in irregular estimate fluctuations over time, complicating practical applications (e.g., Hesarkazzazi et al., 2021; O’Brien & Burn, 2018; Vu & Mishra, 2019). Unlike the SMEV distribution, the nonstationary MEV distribution (Vidrio-Sahagún et al., 2023) could offer a potential solution by addressing both covariate stochasticity and nonstationarity. Future research on improving the integration of physical covariates for nonstationary correction is desired.

It is worth highlighting that extending the introduced correction methods to other hydroclimatic variables is feasible with two key considerations: the suitability of nonstationary treatment and the applicability of the Metastatistical approach. The first consideration is that the variable must exhibit nonstationary behaviour to justify using nonstationary correction methods. The second involves selecting an appropriate distribution for ordinary events and ensuring their statistical independence. For variables with strong autocorrelation, like air temperature and streamflow, additional steps would be required to extract independent ordinary events, such as applying decorrelation procedures or establishing a minimum time spacing between ordinary events.

Furthermore, it is worth noting the different yet complementary goals of the proposed methods and climate model ensembles. Global-regional climate models' ensembles leverage multiple outputs from different models to improve climate projections by seeking a consensus representation of the climate system's probability space while addressing uncertainties (Benestad et al., 2017). However, climate model ensembles may still contain systematic biases, especially for extremes (Ban et al., 2021; Pichelli et al., 2021). Bias correction aims to remove these biases in either individual or ensemble climate projections for local-scale applications. While the proposed methods apply to projections from any modeling source, investigating their benefits for ensembles, which might have lower biases than individual ensemble members, is recommended.

Conclusions

We developed a novel approach to advance the bias correction of precipitation extremes under nonstationarity from two aspects. On one hand, we introduced the NS-QM and CNS-QM methods to account for nonstationarity explicitly and continuously. NS-QM accounts for changing nonstationarity patterns between the baseline and prediction periods, while CNS-QM assumes consistent nonstationarity patterns across periods. On the other hand, we introduced nonstationary SMEV distributions into the bias correction to leverage ordinary-event data, which are often better simulated by climate models. The combination of both innovations led to the NS-QM-SMEV and CNS-QM-SMEV methods. The proposed methods were systematically assessed independently of specific climate model structures/parametrizations in a simulation study using synthetic data sets, taking the conventional QM, DQM, QDM, and CDFt methods and the GEV distribution as benchmarks. The simulation results demonstrate that NS-QM-SMEV and CNS-QM-SMEV are superior to all other methods, as their estimates are more accurate (roughly unbiased and with lower error), less uncertain, and they can capture the entire distribution more effectively. However, these methods require adequate identification of nonstationarity patterns and minimum sample sizes of around 70 years for their reliable parametrization; their superiority is improved with the sample size increase. Although conventional stationary and quasi-stationary methods may show competitive results at some quantiles or distribution characteristics, particularly for short data sets, caution is warranted as they tend to distort quantile-quantile matching, deviating from their theoretical expected functioning. This issue is substantially minimized in NS-QM-SMEV and CNS-QM-SMEV. Moreover, in the practical applications, we also demonstrated that the proposed SMEV-based nonstationary methods generally avoid producing extremely high estimates, large uncertainty, and prominent discontinuity in the estimates. Therefore, the NS-QM-SMEV and CNS-QM-SMEV methods are overall superior and practically viable for bias-correcting nonstationary hydro-climatic extremes.

Appendix

Appendix - A

We adopted the QM, DQM, QDM and CDFt methods as benchmarks. The QM method (Equation 1) may capture nonstationarity partially by assigning varying probabilities to time-varying quantiles based on their evolving positions relative to the stationary distribution. The quasi-stationary DQM, QDM, and CDFt methods indirectly address nonstationarity through relative changes between baseline and prediction periods in the form of a sudden shift in the stationary distributions. Recall that DQM, QDM, and CDFt provide estimates for the prediction period, and they reduce to QM when bias-correcting in the baseline period (see Equations A1–A3). The differences between these four benchmark methods and the NS-QM and CNS-QM methods are illustrated in Figure 1.

The DQM method shifts the estimates to match the relative change in the mean of $${y}_{t}^{s,p}$$ with respect to $${y}_{t}^{s,b}$$. Since we are bias-correcting precipitation, we adopted the scaling and rescaling (multiplicative) form of the DQM method, which is given by: A1 $\begin{align*}{\widehat{y}}_{t}^{s,p}={H}_{{Y}^{o,b}}^{-1}\left({H}_{{Y}^{s,b}}\left[\frac{{\overline{y}}^{s,b}}{{\overline{y}}^{s,p}}{y}_{t}^{s,p}\right]\right)\cdot \frac{{\overline{y}}^{s,p}}{{\overline{y}}^{s,b}}\,\end{align*}$ where $${\overline{y}}^{s,b}$$ and $${\overline{y}}^{s,p}$$ are the means of simulations in the baseline and prediction periods, respectively.

The QDM method is the multiplicative form of the equidistant quantile matching method (Li et al., 2010), which is deemed more realistic for precipitation (Wang & Chen, 2014). The QDM method aims to preserve the relative changes in all quantiles of the simulations between baseline and prediction periods and is given by: A2 $\begin{align*}{\widehat{y}}_{t}^{s,p}={H}_{{Y}^{o,b}}^{-1}\left({H}_{{Y}^{s,p}}\left[{y}_{t}^{s,p}\right]\right)\cdot {{\Delta }}_{t}^{s,p};\quad \text{where}\,{{\Delta }}_{t}^{s,p}=\frac{{y}_{t}^{s,p}}{{H}_{{Y}^{s,b}}^{-1}\left({H}_{{Y}^{s,p}}\left({y}_{t}^{s,p}\right)\right)}\end{align*}$ where the term $${{\Delta }}_{t}^{s,p}$$ accounts for the relative change in simulations' quantiles between baseline and prediction periods.

The CDFt method uses the relationship between distributions of observations and simulations in the baseline period to estimate the observations' distribution in the prediction period. The CDFt method is given by: A3 $\begin{align*}{H}_{{\widehat{Y}}^{s,p}}\left({y}_{t}\right)={H}_{{Y}^{o,b}}\left({H}_{{Y}^{s,b}}^{-1}\left[{H}_{{Y}^{s,p}}\left({y}_{t}^{s,p}\right)\right]\right)\end{align*}$ where $${H}_{{Y}^{o,p}}\left({y}_{t}\right)$$ is then applied to bias-correct $${y}_{t}^{s,p}$$ and produce $${\widehat{y}}_{t}^{s,p}$$.

Acknowledgments

This work was funded by the Flood Hazard Identification and Mapping Program of Environment and Climate Change Canada as well as the Canada Research Chair (Tier 1) awarded to Dr. Pietroniro.

Data Availability Statement

We used data from different sources. The observational data sets from the Global Historical Climatology Network daily (GHCNd) database (Menne et al., 2012) provided by the National Centers for Environmental Information (NCEI) are available at: . The MXN00016027 observational data set at Cuitzeo, spanning up to 2017, was retrieved from the source (CONAGUA) , as the data set in the GHCN database only covers up to 2005. The WRF-MPI-ESM-LR regional climate model data from the NA-CORDEX data archive is freely available at: . CMIP6 global climate model simulations provided by the Earth System Grid Federation (ESGF) are available at: . The MATLAB scripts are available at .

References

Ban, N., Caillaud, C., Coppola, E., Pichelli, E., Sobolowski, S., Adinolfi, M., et al. (2021). The first multi‐model ensemble of regional climate simulations at kilometer‐scale resolution, part I: Evaluation of precipitation. Climate Dynamics, 57(1–2), 275–302. https://doi.org/10.1007/s00382‐021‐05708‐w

Benestad, R., Sillmann, J., Thorarinsdottir, T. L., Guttorp, P., Mesquita, M. D. S., Tye, M. R., et al. (2017). New vigour involving statisticians to overcome ensemble fatigue. Nature Climate Change, 7(10), 697–703. https://doi.org/10.1038/nclimate3393

Bony, S., Stevens, B., Frierson, D. M. W., Jakob, C., Kageyama, M., Pincus, R., et al. (2015). Clouds, circulation and climate sensitivity. Nature Geoscience, 8(4), 261–268. https://doi.org/10.1038/ngeo2398

Boumis, G., Moftakhari, H. R., & Moradkhani, H. (2024). A metastatistical frequency analysis of extreme storm surge hazard along the US coastline. Coastal Engineering Journal, 66(2), 1–15. https://doi.org/10.1080/21664250.2024.2338323

Brunner, M. I., Slater, L., Tallaksen, L. M., & Clark, M. (2021). Challenges in modeling and predicting floods and droughts: A review. Wiley Interdisciplinary Reviews: Water, 8(3), 1–32. https://doi.org/10.1002/wat2.1520

Bürger, G., Sobie, S. R., Cannon, A. J., Werner, A. T., & Murdock, T. Q. (2013). Downscaling extremes: An intercomparison of multiple methods for future climate. Journal of Climate, 26(10), 3429–3449. https://doi.org/10.1175/JCLI‐D‐12‐00249.1

Cannon, A. J., Sobie, S. R., & Murdock, T. Q. (2015). Bias correction of GCM precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes? Journal of Climate, 28(17), 6938–6959. https://doi.org/10.1175/JCLI‐D‐14‐00754.1

Chen, D., Dai, A., & Hall, A. (2021). The convective‐to‐total precipitation ratio and the “drizzling” bias in climate models. Journal of Geophysical Research: Atmospheres, 126(16). https://doi.org/10.1029/2020JD034198

Christensen, J. H., Boberg, F., Christensen, O. B., & Lucas‐Picher, P. (2008). On the need for bias correction of regional climate change projections of temperature and precipitation. Geophysical Research Letters, 35(20). https://doi.org/10.1029/2008GL035694

Cohn, T. A. (2005). Estimating contaminant loads in rivers: An application of adjusted maximum likelihood to type 1 censored data. Water Resources Research, 41(7), 1–13. https://doi.org/10.1029/2004WR003833

Coles, S. (2001). An introduction to statistical modeling of extreme values. Springer‐Verlag London.

De Michele, C., & Avanzi, F. (2018). Superstatistical distribution of daily precipitation extremes: A worldwide assessment. Scientific Reports, 8(1), 1–11. https://doi.org/10.1038/s41598‐018‐31838‐z

Efron, B. (1992). Bootstrap methods: Another look at the Jackknife. In Breakthroughs in statistics (Springer S). Springer.

el Adlouni, S., Ouarda, T. B. M. J., Zhang, X., Roy, R., & Bobée, B. (2007). Generalized maximum likelihood estimators for the nonstationary generalized extreme value model. Water Resources Research, 43(3), 1–13. https://doi.org/10.1029/2005WR004545

Emmanouil, S., Langousis, A., Nikolopoulos, E. I., & Anagnostou, E. N. (2023). Exploring the future of rainfall extremes over CONUS: The effects of high emission climate change trajectories on the intensity and frequency of rare precipitation events. Earth's Future, 11(4). https://doi.org/10.1029/2022EF003039

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., & Taylor, K. E. (2016). Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization. Geoscientific Model Development, 9(5), 1937–1958. https://doi.org/10.5194/gmd‐9‐1937‐2016

Falkensteiner, M.‐A., Schellander, H., Ehrensperger, G., & Hell, T. (2023). Accounting for seasonality in the metastatistical extreme value distribution. Weather and Climate Extremes, 42, 100601. https://doi.org/10.1016/j.wace.2023.100601

Giorgetta, M. A., Jungclaus, J., Reick, C. H., Legutke, S., Bader, J., Böttinger, M., et al. (2013). Climate and carbon cycle changes from 1850 to 2100 in MPI‐ESM simulations for the Coupled Model Intercomparison Project phase 5. Journal of Advances in Modeling Earth Systems, 5(3), 572–597. https://doi.org/10.1002/jame.20038

Gomez‐Garcia, M., Matsumura, A., Ogawada, D., & Takahashi, K. (2019). Time scale decomposition of climate and correction of variability using synthetic samples of stable distributions. Water Resources Research, 55(5), 3632–3658. https://doi.org/10.1029/2018WR023053

Gründemann, G. J., Zorzetto, E., Beck, H. E., Schleiss, M., van de Giesen, N., Marani, M., & van der Ent, R. J. (2023). Extreme precipitation return levels for multiple durations on a global scale. Journal of Hydrology, 621, 129558. https://doi.org/10.1016/j.jhydrol.2023.129558

Helsel, D. R. (2011). Three approaches for censored data. In Statistics for censored environmental data using Minitab® and R (pp. 12–21). John Wiley & Sons, Inc. https://doi.org/10.1002/9781118162729.ch2

Heo, J. H., Ahn, H., Shin, J. Y., Kjeldsen, T. R., & Jeong, C. (2019). Probability distributions for a quantile mapping technique for a bias correction of precipitation data: A case study to precipitation data under climate change. Water (Switzerland), 11(7), 1475. https://doi.org/10.3390/w11071475

Hernanz, A., García‐Valero, J. A., Domínguez, M., Ramos‐Calzado, P., Pastor‐Saavedra, M. A., & Rodríguez‐Camino, E. (2022). Evaluation of statistical downscaling methods for climate change projections over Spain: Present conditions with perfect predictors. International Journal of Climatology, 42(2), 762–776. https://doi.org/10.1002/joc.7271

Hertig, E., Maraun, D., Bartholy, J., Pongracz, R., Vrac, M., Mares, I., et al. (2019). Comparison of statistical downscaling methods with respect to extreme events over Europe: Validation results from the perfect predictor experiment of the COST Action VALUE. International Journal of Climatology, 39(9), 3846–3867. https://doi.org/10.1002/joc.5469

Hesarkazzazi, S., Arabzadeh, R., Hajibabaei, M., Rauch, W., Kjeldsen, T. R., Prosdocimi, I., et al. (2021). Stationary vs non‐stationary modelling of flood frequency distribution across northwest England. Hydrological Sciences Journal, 00(00), 1–16. https://doi.org/10.1080/02626667.2021.1884685

Hui, Y., Xu, Y., Chen, J., Xu, C. Y., & Chen, H. (2020). Impacts of bias nonstationarity of climate model outputs on hydrological simulations. Hydrology Research, 51(5), 925–941. https://doi.org/10.2166/nh.2020.254

Kallache, M., Vrac, M., Naveau, P., & Michelangeli, P. A. (2011). Nonstationary probabilistic downscaling of extreme precipitation. Journal of Geophysical Research, 116(5), D05113. https://doi.org/10.1029/2010JD014892

Kendall, M. G. (1975). Rank correlation methods, book series, Charles Griffin. Oxford University Press.

Kharin, V. V., Zwiers, F. W., Zhang, X., & Wehner, M. (2013). Changes in temperature and precipitation extremes in the CMIP5 ensemble. Climatic Change, 119(2), 345–357. https://doi.org/10.1007/s10584‐013‐0705‐8

Kim, S., Joo, K., Kim, H., Shin, J. Y., & Heo, J. H. (2021). Regional quantile delta mapping method using regional frequency analysis for regional climate model precipitation. Journal of Hydrology, 596, 125685. https://doi.org/10.1016/j.jhydrol.2020.125685

Koutsoyiannis, D., & Montanari, A. (2015). Negligent killing of scientific concepts: The stationarity case. Hydrological Sciences Journal, 60(7–8), 1174–1183. https://doi.org/10.1080/02626667.2014.959959

Lafferty, D. C., & Sriver, R. L. (2023). Downscaling and bias‐correction contribute considerable uncertainty to local climate projections in CMIP6. Npj Climate and Atmospheric Science, 6(1), 158. https://doi.org/10.1038/s41612‐023‐00486‐0

Lange, S. (2019). Trend‐preserving bias adjustment and statistical downscaling with ISIMIP3BASD (v1.0). Geoscientific Model Development, 12(7), 3055–3070. https://doi.org/10.5194/gmd‐12‐3055‐2019

Lanzante, J. R., Adams‐Smith, D., Dixon, K. W., Nath, M., & Whitlock, C. E. (2020). Evaluation of some distributional downscaling methods as applied to daily maximum temperature with emphasis on extremes. International Journal of Climatology, 40(3), 1571–1585. https://doi.org/10.1002/joc.6288

Lanzante, J. R., Nath, M. J., Whitlock, C. E., Dixon, K. W., & Adams‐Smith, D. (2019). Evaluation and improvement of tail behaviour in the cumulative distribution function transform downscaling method. International Journal of Climatology, 39(4), 2449–2460. https://doi.org/10.1002/joc.5964

Li, H., Sheffield, J., & Wood, E. F. (2010). Bias correction of monthly precipitation and temperature fields from Intergovernmental Panel on Climate Change AR4 models using equidistant quantile matching. Journal of Geophysical Research, 115(10). https://doi.org/10.1029/2009JD012882

Mann, H. B. (1945). Nonparametric tests against trend. Econometrica: Journal of the Econometric Society, 13(3), 245–259. https://doi.org/10.2307/1907187

Maraun, D. (2013). Bias correction, quantile mapping, and downscaling: Revisiting the inflation issue. Journal of Climate, 26(6), 2137–2143. https://doi.org/10.1175/JCLI‐D‐12‐00821.1

Maraun, D. (2016). Bias correcting climate change simulations ‐ A critical review. Current Climate Change Reports, 2(4), 211–220. https://doi.org/10.1007/s40641‐016‐0050‐x

Maraun, D., & Widmann, M. (2017). Statistical downscaling and bias correction for climate research. Statistical downscaling and bias correction for climate research. Cambridge University Press. https://doi.org/10.1017/9781107588783

Marra, F., Amponsah, W., & Papalexiou, S. M. (2023). Non‐asymptotic Weibull tails explain the statistics of extreme daily precipitation. Advances in Water Resources, 173, 104388. https://doi.org/10.1016/j.advwatres.2023.104388

Marra, F., Borga, M., & Morin, E. (2020). A unified framework for extreme subdaily precipitation frequency analyses based on ordinary events. Geophysical Research Letters, 47(18), 1–8. https://doi.org/10.1029/2020GL090209

Marra, F., Nikolopoulos, E. I., Anagnostou, E. N., & Morin, E. (2018). Metastatistical Extreme Value analysis of hourly rainfall from short records: Estimation of high quantiles and impact of measurement errors. Advances in Water Resources, 117(May), 27–39. https://doi.org/10.1016/j.advwatres.2018.05.001

Marra, F., Zoccatelli, D., Armon, M., & Morin, E. (2019). A simplified MEV formulation to model extremes emerging from multiple nonstationary underlying processes. Advances in Water Resources, 127(April), 280–290. https://doi.org/10.1016/j.advwatres.2019.04.002

Massey Jr, F. J. (1951). The Kolmogorov‐Smirnov test for goodness of fit. Journal of the American Statistical Association, 46(253), 68–78. https://doi.org/10.2307/2280095

Maurer, E. P., & Pierce, D. W. (2014). Bias correction can modify climate model simulated precipitation changes without adverse effect on the ensemble mean. Hydrology and Earth System Sciences, 18(3), 915–925. https://doi.org/10.5194/hess‐18‐915‐2014

Mauritsen, T., Bader, J., Becker, T., Behrens, J., Bittner, M., Brokopf, R., et al. (2019). Developments in the MPI‐M Earth system model version 1.2 (MPI‐ESM1.2) and its response to increasing CO2. Journal of Advances in Modeling Earth Systems, 11(4), 998–1038. https://doi.org/10.1029/2018MS001400

Mearns, L., McGinnis, S., Korytina, D., Arritt, R., Biner, S., Bukovsky, M., et al. (2017). The NA‐CORDEX dataset, version 1.0. NCAR Climate Data Gateway. https://doi.org/10.5065/D6SJ1JCH

Menne, M. J., Durre, I., Vose, R. S., Gleason, B. E., & Houston, T. G. (2012). An overview of the global historical climatology network‐daily database. Journal of Atmospheric and Oceanic Technology, 29(7), 897–910. https://doi.org/10.1175/JTECH‐D‐11‐00103.1

Merkenschlager, C., Hertig, E., & Jacobeit, J. (2017). Non‐stationaries in the relationships of heavy precipitation events in the Mediterranean area and the large‐scale circulation in the second half of the 20th century. Global and Planetary Change, 151, 108–121. https://doi.org/10.1016/j.gloplacha.2016.10.009

Merz, B., Vorogushyn, S., Uhlemann, S., Delgado, J., & Hundecha, Y. (2012). HESS opinions: “More efforts and scientific rigour are needed to attribute trends in flood time series”. Hydrology and Earth System Sciences, 16(5), 1379–1387. https://doi.org/10.5194/hess‐16‐1379‐2012

Michalek, A. T., Villarini, G., & Kim, T. (2024). Understanding the impact of precipitation bias‐correction and statistical downscaling methods on projected changes in flood extremes. Earth's Future, 12(3). https://doi.org/10.1029/2023EF004179

Michelangeli, P. A., Vrac, M., & Loukos, H. (2009). Probabilistic downscaling approaches: Application to wind cumulative distribution functions. Geophysical Research Letters, 36(11). https://doi.org/10.1029/2009GL038401

Milly, A. P. C. D., Betancourt, J., Falkenmark, M., Hirsch, R. M., Zbigniew, W., Lettenmaier, D. P., et al. (2008). Stationarity is dead: Whither water management? Science, 319(5863), 573–574. https://doi.org/10.1126/science.1151915

Milly, P. C. D., Betancourt, J., Falkenmark, M., Hirsch, R. M., Kundzewicz, Z. W., Lettenmaier, D. P., et al. (2015). On critiques of “stationarity is dead: Whither water management?”. Water Resources Research, 51(9), 7785–7789. https://doi.org/10.1002/2015WR017408

Miniussi, A., & Marani, M. (2020). Estimation of daily rainfall extremes through the metastatistical extreme value distribution: Uncertainty minimization and implications for trend detection. Water Resources Research, 56(7). https://doi.org/10.1029/2019WR026535

Miniussi, A., Marani, M., & Villarini, G. (2020). Metastatistical extreme value distribution applied to floods across the continental United States. Advances in Water Resources, 136(December 2019), 103498. https://doi.org/10.1016/j.advwatres.2019.103498

Miniussi, A., & Marra, F. (2021). Estimation of extreme daily precipitation return levels at‐site and in ungauged locations using the simplified MEV approach. Journal of Hydrology, 603(PB), 126946. https://doi.org/10.1016/j.jhydrol.2021.126946

Mitovski, T., Cole, J. N. S., McFarlane, N. A., Von Salzen, K., & Zhang, G. J. (2019). Convective response to large‐scale forcing in the tropical western Pacific simulated by spCAM5 and CanAM4.3. Geoscientific Model Development, 12(5), 2107–2117. https://doi.org/10.5194/gmd‐12‐2107‐2019

Moustakis, Y., Papalexiou, S. M., Onof, C. J., & Paschalis, A. (2021). Seasonality, intensity, and duration of rainfall extremes change in a warmer climate. Earth's Future, 9(3). https://doi.org/10.1029/2020EF001824

O’Brien, N. L., & Burn, D. H. (2018). A nonstationary peaks‐over‐threshold approach for modelling daily precipitation with covariate‐dependent thresholds. Canadian Water Resources Journal, 1784(May), 1–24. https://doi.org/10.1080/07011784.2018.1455538

Ouarda, T. B. M. J., Charron, C., & St‐Hilaire, A. (2019). Uncertainty of stationary and nonstationary models for rainfall frequency analysis. International Journal of Climatology, 40(September), 1–20. https://doi.org/10.1002/joc.6339

Ouarda, T. B. M. J., Yousef, L. A., & Charron, C. (2018). Non‐stationary intensity‐duration‐frequency curves integrating information concerning teleconnections and climate change. International Journal of Climatology, 39(August 2018), 2306–2323. https://doi.org/10.1002/joc.5953

Panofsky, H. A., & Brier, G. W. (1968). Some applications of statistics to meteorology (p. 224). The Pennsylvania State University.

Pfahl, S., O’Gorman, P. A., & Fischer, E. M. (2017). Understanding the regional pattern of projected future changes in extreme precipitation. Nature Climate Change, 7(6), 423–427. https://doi.org/10.1038/nclimate3287

Pichelli, E., Coppola, E., Sobolowski, S., Ban, N., Giorgi, F., Stocchi, P., et al. (2021). The first multi‐model ensemble of regional climate simulations at kilometer‐scale resolution part 2: Historical and future simulations of precipitation. Climate Dynamics, 56(11–12), 3581–3602. https://doi.org/10.1007/s00382‐021‐05657‐4

Polade, S. D., Pierce, D. W., Cayan, D. R., Gershunov, A., & Dettinger, M. D. (2014). The key role of dry days in changing regional climate and precipitation regimes. Scientific Reports, 4(1), 4364. https://doi.org/10.1038/srep04364

Ragno, E., AghaKouchak, A., Cheng, L., & Sadegh, M. (2019). A generalized framework for process‐informed nonstationary extreme value analysis. Advances in Water Resources, 130(November 2018), 270–282. https://doi.org/10.1016/j.advwatres.2019.06.007

Rajczak, J., & Schär, C. (2017). Projections of future precipitation extremes over Europe: A multimodel assessment of climate simulations. Journal of Geophysical Research: Atmospheres, 122(20), 10773–10800. https://doi.org/10.1002/2017JD027176

Rummukainen, M. (1997). Methods for statistical downscaling of GCM simulations. SMHI.

Serinaldi, F., & Kilsby, C. G. (2014). Rainfall extremes: Toward reconciliation after the battle of distributions. Water Resources Research, 50(1), 336–352. https://doi.org/10.1002/2013WR014211

Serinaldi, F., & Kilsby, C. G. (2015). Stationarity is undead: Uncertainty dominates the distribution of extremes. Advances in Water Resources, 77, 17–36. https://doi.org/10.1016/j.advwatres.2014.12.013

Serinaldi, F., Kilsby, C. G., & Lombardo, F. (2018). Untenable nonstationarity: An assessment of the fitness for purpose of trend tests in hydrology. Advances in Water Resources, 111(June 2017), 132–155. https://doi.org/10.1016/j.advwatres.2017.10.015

Sherwood, S. C., Bony, S., & Dufresne, J. L. (2014). Spread in model climate sensitivity traced to atmospheric convective mixing. Nature, 505(7481), 37–42. https://doi.org/10.1038/nature12829

Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Barker, D. M., Duda, M. G., et al. (2008). A description of the advanced research WRF version 3. NCAR Technical Note, 475, 113.

Slater, L. J., Anderson, B., Buechel, M., Dadson, S., Han, S., Harrigan, S., et al. (2021). Nonstationary weather and water extremes: A review of methods for their detection, attribution, and management. Hydrology and Earth System Sciences, 25(7), 3897–3935. https://doi.org/10.5194/hess‐25‐3897‐2021

Srivastav, R. K., Schardong, A., & Simonovic, S. P. (2014). Equidistance quantile matching method for updating IDFCurves under climate change. Water Resources Management, 28(9), 2539–2562. https://doi.org/10.1007/s11269‐014‐0626‐y

Stephens, G. L., L’Ecuyer, T., Forbes, R., Gettelmen, A., Golaz, J., Bodas‐Salcedo, A., et al. (2010). Dreary state of precipitation in global models. Journal of Geophysical Research, 115(D24). https://doi.org/10.1029/2010jd014532

Teutschbein, C., & Seibert, J. (2012). Bias correction of regional climate model simulations for hydrological climate‐change impact studies: Review and evaluation of different methods. Journal of Hydrology, 456–457, 12–29. https://doi.org/10.1016/j.jhydrol.2012.05.052

Trenberth, K. E. (2011). Changes in precipitation with climate change. Climate Research, 47(1–2), 123–138. https://doi.org/10.3354/cr00953

Um, M. J., Kim, H., & Heo, J. H. (2016). Hybrid approach in statistical bias correction of projected precipitation for the frequency analysis of extreme events. Advances in Water Resources, 94, 278–290. https://doi.org/10.1016/j.advwatres.2016.05.021

Van De Velde, J., Demuzere, M., De Baets, B., & Verhoest, N. E. C. (2022). Impact of bias nonstationarity on the performance of uni‐ and multivariate bias‐adjusting methods: A case study on data from Uccle, Belgium. Hydrology and Earth System Sciences, 26(9), 2319–2344. https://doi.org/10.5194/hess‐26‐2319‐2022

Vidrio‐Sahagún, C. T., & He, J. (2022a). Hydrological frequency analysis under nonstationarity using the Metastatistical approach and its simplified version. Advances in Water Resources, 166(June), 104244. https://doi.org/10.1016/j.advwatres.2022.104244

Vidrio‐Sahagún, C. T., & He, J. (2022b). The decomposition‐based nonstationary flood frequency analysis. Journal of Hydrology, 612, 128186. https://doi.org/10.1016/j.jhydrol.2022.128186

Vidrio‐Sahagún, C. T., He, J., & Pietroniro, A. (2023). Nonstationary hydrological frequency analysis using the Metastatistical extreme value distribution. Advances in Water Resources, 176, 104460. https://doi.org/10.1016/j.advwatres.2023.104460

Vidrio‐Sahagún, C. T., Ruschkowski, J., He, J., & Pietroniro, A. (2024). A practice‐oriented framework for stationary and nonstationary flood frequency analysis. Environmental Modelling & Software, 173, 105940. https://doi.org/10.1016/j.envsoft.2024.105940

Vincent, L. A., Zhang, X., Mekis, É., Wan, H., & Bush, E. J. (2018). Changes in Canada’s climate: Trends in indices based on daily temperature and precipitation data. Atmosphere‐Ocean, 56(5), 332–349. https://doi.org/10.1080/07055900.2018.1514579

Vu, T. M., & Mishra, A. K. (2019). Nonstationary frequency analysis of the recent extreme precipitation events in the United States. Journal of Hydrology, 575(March), 999–1010. https://doi.org/10.1016/j.jhydrol.2019.05.090

Wang, L., & Chen, W. (2014). Equiratio cumulative distribution function matching as an improvement to the equidistant approach in bias correction of precipitation. Atmospheric Science Letters, 15(1), 1–6. https://doi.org/10.1002/asl2.454

Wilson, P. S., & Toumi, R. (2005). A fundamental probability distribution for heavy rainfall. Geophysical Research Letters, 32(14), 1–4. https://doi.org/10.1029/2005GL022465

Wootten, A., Terando, A., Reich, B. J., Boyles, R. P., & Semazzi, F. (2017). Characterizing sources of uncertainty from global climate models and downscaling techniques. Journal of Applied Meteorology and Climatology, 56(12), 3245–3262. https://doi.org/10.1175/JAMC‐D‐17‐0087.1

Yuan, F., Zhao, C., Jiang, Y., Ren, L., Shan, H., Zhang, L., et al. (2017). Evaluation on uncertainty sources in projecting hydrological changes over the Xijiang River basin in South China. Journal of Hydrology, 554, 434–450. https://doi.org/10.1016/j.jhydrol.2017.08.034

Zhang, X., Alexander, L., Hegerl, G. C., Jones, P., Tank, A. K., Peterson, T. C., et al. (2011). Indices for monitoring changes in extremes based on daily temperature and precipitation data. In Wiley interdisciplinary reviews: Climate change. Wiley‐Blackwell. https://doi.org/10.1002/wcc.147