# On the attribution of weather events to climate change using a fit to
extreme value distributions

###### Abstract

Increases in extreme weather events are a potentially important consequence of
anthropogenic climate change (ACC), yet, are difficult to attribute to ACC
because the record length is often similar to, or shorter than, extreme-event
return periods. This study is motivated by the “World Weather Attribution
Project” (WWAP) and their approach of fitting extreme value distribution
functions to local observations. The approach calculates the dependence of
distribution parameters on the global mean surface temperature (GMST) and uses
this dependence to attribute extreme events to ACC. Applying the WWAP method
to a large ensemble of climate simulations run without anthropogenic forcing,
we still find a strong dependence of distribution parameters on GMST. This
dependence results from internal climate variability, such as ENSO, affecting
both extreme events and GMST. Therefore, dependence on GMST does not
necessarily imply an effect of ACC on extremes.
We next re-examine three WWAP attribution cases. We consider whether an
extreme value, normal, or log-normal distribution better represents the data;
if a GMST-dependence of distribution parameters is justified using a
likelihood ratio test; and if a meaningful attribution can be made given
errors in GMST dependence. The effects of natural variability on both GMST and
extremes make it impossible to attribute the 2020 Siberian Heatwave and
Australia’s 2020–2021 bushfires to ACC. The small number of data points for
the 2019–2021 drought in Madagascar precludes a meaningful attribution
analysis. Overall, natural variability and the uncertain relationship between
GMST and extremes make attribution using the WWAP approach challenging.

##  1 Introduction

The Intergovernmental Panel on Climate Change Grubb et al. (2022) finds that
global mean surface temperature (GMST) has risen 1.26 ∘C above preindustrial
levels as of the year 2020. While the changes are manifest clearly in the mean
climate state, and there are reasons to think extreme events may change with
the mean climate, it is a far more complex task to attribute individual
extreme weather events to anthropogenic climate change. Extreme events,
whether heat waves, droughts, floods, etc., are rare by definition, so the
observational record may simply be too short to detect anthropogenic signals.
This issue is especially pervasive in lower-income countries where data are
sparse, and the historical record is short Otto et al. (2020). The situation
is particularly problematic for extreme precipitation attribution, which is
highly sensitive to the specific timescales and percentiles used to define
extreme precipitation events Pendergrass (2018).

The field of detection and attribution of extreme events has progressed
significantly over the past decade. The purpose of attribution generally is to
assess how human-induced climate change has influenced the probability or
magnitude of an extreme event. This was done using the concept of the fraction
of attributable risk (FAR, Allen, 2003) for the 2003 European heatwave Stott
et al. (2004). This approach leveraged coupled climate models to compare the
probability of the 2003 heatwave under anthropogenic influences relative to
the probability of the heatwave without them. This was followed by further
attribution studies that took advantage of simulations from coupled climate
models to estimate human influences on extreme events, using constraints
derived from optimal fingerprinting techniques Christidis et al. (2015) and
using climate models to validate observed PDFs of precipitation Perkins et al.
(2007).

There are also approaches that specifically use observations. The “Analogue-
Based Approach” takes the circulation patterns that coincide with an extreme
event and finds days with similar circulation patterns in the historical
record Lorenz (1969); Yiou et al. (2007); Vautard and Yiou (2009). Another
methodology fits statistical distributions to the observational record in what
Stott et al. (2016) refer to as “Empirical-Based Approaches.” This approach
assesses how the return period of particular events has changed over the
observed record, with the purpose of detecting potential trends in the
frequency of extreme events. These observational approaches add to our
understanding of extreme events and have been used to study trends in extreme
events such as flooding in Thailand in 2011 Van Oldenborgh et al. (2012).
Stott et al. (2016) note the importance of also including climate models in
the analysis. Climate models allow a better estimate of the magnitude and
direction of the human-induced contribution to an extreme event than can be
done using observations because models can be run without anthropogenic
emissions to aid the analysis.

Another use of models is based on single-model initial condition large
ensembles (LE), composed of a climate model run several times with slightly
different initial conditions. Averaging over multiple ensemble members can
reduce the signal of internal variability, allowing us to identify the signal
of climate change. This generally improves the detection of trends in the
magnitude or frequency of an extreme event, which are not captured well in a
single simulation or in the observational record where the signal from
internal variability may be dominant Yoshida et al. (2017); Yamaguchi et al.
(2020). These simulations have been used in the context of tropical cyclone
slowdown Yamaguchi et al. (2020), as well as drought events in South and East
Africa Lott et al. (2013); Pascale et al. (2020). The detection and
attribution literature is extensive, and for a further detailed discussion, we
refer the reader to Stott et al. (2016), Knutson et al. (2017), and
Seneviratne et al. (2021).

One particular approach that has been used in recent years as a rapid response
to extreme events attempts to fit the observed data to an extreme value
probability distribution function. The fit allows for a GMST dependence of the
distribution parameters (e.g., Ciavarella et al., 2021; Van Oldenborgh et al.,
2021; Harrington et al., 2022). A GMST dependence is interpreted as a response
of extreme events to anthropogenic global warming. This dependence is then
used to attribute extreme events and quantify the effects of anthropogenic
forcing on the probability of their occurrence. This approach has been used by
the World Weather Attribution Project (WWAP,
https://www.worldweatherattribution.org/). The WWAP also makes use of climate
model simulations in interpreting extreme events, but our focus here will be
on the more novel aspects of their approach to the empirical inference of
extreme value attribution.

Motivated by these pioneering efforts based on the fit to extreme event PDFs,
we reexamine three test cases analyzed by the WWAP: the Siberian Heatwave
analysis of Ciavarella et al. (2021, hereafter CC2021), the Australian
Bushfire analysis of Van Oldenborgh et al. (2021, hereafter OK2021) and the
Madagascar Drought analysis of Harrington et al. (2022, hereafter HW2022).
Each of these studies leverages either a Generalized Extreme Value (GEV) or
Generalized Pareto Distribution (GPD) fit to assess potential contributions to
extreme events from anthropogenic climate change. While these studies do
compute probabilities from climate model simulations, the primary findings are
based on observations, which we also focus on in our analysis, adding Large
Ensemble analyses into the mix. In each of the three cases, we attempt to
evaluate the merits and weaknesses associated with the proposed approaches and
use the CESM Large Ensemble Kay et al. (2015) to put the use of the observed
instrumental record, which is short, into perspective.

Our study addresses four questions: (1) Is a GMST dependence of the PDF
parameters necessarily an indication of the signal of anthropogenic climate
change? (2) Is a GEV or GPD necessary to fit the available observations, or is
a more standard distribution (normal/log-normal) sufficient? (3) Is the
addition of GMST-dependent distribution parameters, which form the basis for
the attribution analysis, justified statistically, and is the error level in
this dependence sufficiently small to allow attribution? (4) What do Single-
Model Initial-Condition Large Ensembles tell us about the amount of data
needed to determine the GMST dependence with confidence?

The rest of the paper is structured as follows. Section 2 describes the
statistical models, observation sets, and methodologies. Section 3a shows that
a GMST-dependence of distribution parameters does not necessarily represent
the effects of anthropogenic climate change using an analysis of a CESM
control simulation (i.e., without anthropogenic emissions). We then examine
the three attribution cases in the remainder of section 3. Section 4 provides
a summary and discussion of the results.

##  2 Methods

The first two attribution cases examined here involve the recent extreme high-
temperature events over Siberia (following CC2021) and Australia (following
OK2021) and use a fit to a Generalized Extreme Value (GEV) distribution. The
third case analyzes a two-year drought event over Madagascar using a
Generalized Pareto Distribution (GPD) model following HW2022. In this section,
we briefly describe these distributions (sections 2.1 and 2.2), and discuss
uncertainty estimates and empirical cumulative distribution functions (CDFs)
used in our analysis (2.3). Adding GMST-dependence of distribution parameters
means having more parameters to fit the observations, which necessarily leads
to a better fit. But that does not necessarily mean the additional parameter
is justified by the observations. We apply the deviance statistic/likelihood
ratio test to examine the justification for adding such GMST-dependence
parameters (2.4). The data used are described in section (2.6).

###  2.1 The Generalized Extreme Value (GEV) distribution

Extreme value distributions are a family of probability distribution functions
(PDFs) that can be used to represent the statistical behavior of block maxima
or minima of a record, e.g., of an annual time series of the maximum daily
temperature at each year (the “block” being one year in this case). It can be
shown that the distribution of such maxima follows one of three classes of
PDFs: Gumbel, Fréchet, or Weibull Coles et al. (2001), depending on the shape
of the tail of the distribution of the events for which block maxima are
calculated. Rather than explicitly specifying one of these three classes in
fitting extreme data, one can combine these distributions into one functional
form, represented by the Generalized Extreme Value (GEV) distribution. Large
quantities of data are typically needed to accurately fit a GEV and calculate
its parameters (e.g., Trevino et al., 2020).

The PDF of a variable $x$ whose statistics are governed by the GEV is given
by,

$GEV(x,\mu,\sigma,\xi)=\frac{1}{\sigma}t(x,\mu,\sigma,\xi)^{\xi+1}e^{-t(x,\mu,\sigma,\xi)},$
(1)

where

$t(x,\mu,\sigma,\xi)=\left(1+\xi\frac{x-\mu}{\sigma}\right)^{-1/\xi},$ (2)

under the assumption that $\xi\neq 0$.

In some of the attribution studies we follow (specifically, CC2021 and
OK2021), the location parameter $\mu$ is assumed to vary linearly with the
global mean surface temperature $T$, as,

$\mu=\mu_{0}+\alpha T.$ (3)

A statistically robust dependence on the GMST is considered in these studies
an indication that the events are affected by anthropogenic climate change,
and we examine this issue below. The determination of the parameter $\alpha$
is then the key to the attribution of extreme events. Notably, the
determination of the $\alpha$ component may be limited due to the short
observational record Zwiers et al. (2011); De Paola et al. (2018). To estimate
the mean $\mu_{0}$, scale parameter $\sigma$, the dependence on the global
mean surface temperature $\alpha$, and the shape parameter $\xi$, we follow
Coles et al. (2001) and the WWAP and maximize the log-likelihood of the data
points $x_{t}$ using a GEV distribution,

$\displaystyle l(\mu_{0},\sigma,\xi,\alpha)$
$\displaystyle=\sum_{t=1}^{m}\log\left(GEV(x_{t},\mu_{0},\sigma,\xi,\alpha)\right).$
(4)

We apply this approach, following CC2021 and OK2021, to observations and to
data from the CESM1 Large Ensemble Kay et al. (2015). The CESM1 Large Ensemble
contains 42 simulations of the climate over the period 1920–2100 under
historical and RCP 8.5 external forcing, of which we use the 35 that were run
on the NCAR supercomputer. Simulations differ only in small perturbations to
the initial conditions. We also perform this analysis with a 1,800-year CESM1
preindustrial control run.

###  2.2 Generalized Pareto Distribution (GPD)

In some applications, one is interested only in points in a time series that
exceed a specified threshold, whose statistics are represented by the
Generalized Pareto Distribution Coles et al. (2001). In order to study the
extreme drought events in Madagascar using a record of precipitation, HW2022
fit a Generalized Pareto Distribution (GPD) to the bottom 20th percentile of
precipitation. Low tails corresponding to precipitation drought years
belonging to the bottom 20th percentile were converted to high tails by
multiplying the precipitation data by $-1$. Although this converts the lowest
precipitation values to maxima, the transformed data have a hard maximum at
zero. We note that the validity of using a GPD, in this case, is questionable
because no such hard maximum is allowed by—or represented in the functional
form of—the GPD.

The PDF of the GPD is given by

$GPD(x,\mu,\sigma,\xi)=\frac{1}{\sigma}\left(1+\xi\frac{x-\mu}{\sigma}\right)^{-(\frac{1}{\xi}+1)}.$
(5)

HW2022 assumed that the location parameter $\mu$ and scale parameter $\sigma$
both depend exponentially on the 4-year running-averaged GMST, $\overline{T}$

$\displaystyle\mu$
$\displaystyle=\mu_{0}\exp({\alpha\overline{T}}/{\mu_{0}}),$
$\displaystyle\sigma$
$\displaystyle=\sigma_{0}\exp({\alpha\overline{T}}/{\mu_{0}}).$ (6)

While this dependence may have been intended to reflect the exponential
dependence on temperature in the Clausius-Clapeyron relation, we find below
the above dependence to be effectively linear given the values of the
calculated parameters $\alpha$ and $\mu$, when reproducing the observational
analysis results of the WWAP. We obtain the optimized GPD parameters
($\mu_{0},\sigma,\xi,\alpha$) by again maximizing the log-likelihood (Eq.
4.10, Coles et al., 2001),

$l(\mu,\sigma,\xi)=\sum_{t=1}^{m}\log(GPD(x_{t},\mu,\sigma,\xi)),$ (7)

using the trust-region constrained algorithm in the Scipy minimize Python
package (‘trust-constr’ method) applied to the negative of the log-likelihood.

The GPD requires the bound $\sigma>0$ and consistency relations of the form
$x\geq\mu$ (for $\xi\geq 0$), and $\mu\leq x$ and $x\leq\mu-\sigma/\xi$ (for
$\xi<0$). As the optimization iteratively attempts to search for the optimal
parameter values, these constraints are occasionally violated, leading to
complex numbers or NaNs (not a number) in the likelihood. Due to the
possibility of encountering such complex numbers or NaNs, we find the GPD
optimization is notably more prone to failure and sensitive to the initial
guesses of $\alpha,\mu_{0},\xi,\sigma$ than the GEV optimization. In order to
avoid such search failures, we replace complex/NaN terms in the log-likelihood
sum with large penalty terms whenever they occur due to constraints being
violated. In addition, if the optimization fails, we re-run the optimization
again with new randomly-selected initial guesses until convergence to a
solution that maximizes log-likelihood and that satisfies the constraints
occurs.

###  2.3 Uncertainty Estimates and empirical CDFs

We follow the WWAP studies and use non-parametric bootstrapping to estimate
uncertainty ranges for the optimized parameters. In particular, we use the
uncertainty in $\alpha$ as one way to evaluate the validity of the
attribution. We also compare the appropriately fitted distributions (GEV and
normal distributions for extreme temperatures or GPD and log-normal for
droughts) to the empirical cumulative distribution functions (CDFs) of the
observations to assess the quality of the fit. The empirical CDF is calculated
by sorting extreme values (either temperature or precipitation) from smallest
to largest, calculating the cumulative sum, and normalizing to a maximum CDF
value of one. We calculate error bars in the empirical CDFs based on bootstrap
resampling with replacement (such that multiple of the same data point may be
drawn) of the observations. The empirical CDF is calculated for each sampling,
and we then calculate the 90th and 95th percentiles from the resulting
distribution at each value of the variable whose CDF is estimated. In order to
estimate uncertainty in the parameters estimated from our GEV or GPD maximum
likelihood fit, we similarly perform a 5,000-bootstrap resampling of the GMST
and extreme temperature. For the Madagascar precipitation case, we perform
5,000 bootstrap resamples when analyzing the observations, but 1,000 bootstrap
resamples when analyzing the CESM-Large Ensemble due to the slow convergence
and expensive computations using the GPD log-likelihood maximization.

###  2.4 Deviance Statistics (Likelihood Ratio Test)

In general, when fitting a model to data, if the number of fit parameters is
increased, one expects a better fit because of the greater flexibility of the
fitted model. Specifically, in the case considered here, adding a parameter
$\alpha$ to include the GMST dependence of the location parameter will
increase the maximum log-likelihood of the data. The salient question is
whether the improvement in the fit exceeds that expected simply on account of
the larger number of model parameters. To evaluate the improvement of the fit,
Coles et al. (2001, theorem 2.7) recommend the use of a deviance statistic,
also known as a likelihood ratio test. This test uses a measure, $D$, equal to
twice the difference of the sum of log-likelihoods of the model with $k$
additional parameters minus that of the simpler model. Because more parameters
lead to a better fit and thus higher log-likelihood, $D$ is non-negative. $D$
is then compared to a $\chi_{k}^{2}$ distribution to assess if the $k$
additional parameters significantly improve the likelihood. For example, with
$k=1$, the 95th percentile of $\chi_{1}^{2}$ is equal to 3.8. If $D$ is larger
than this threshold, we conclude that adding an additional parameter
significantly improved the fit. It can then be inferred that the added
parameter—in our specific case, the GMST dependence, $\alpha$—meaningfully
explains features of the data. We also estimate the distribution and
uncertainty of $D$ through bootstrap resampling of the data 5,000 times.

###  2.5 Mean residual life plot

In the process of fitting a GPD, a minimum threshold is chosen to select the
data that represent extreme values that should be fit. The GPD fit can be
sensitive to the selection of this threshold, where too low of a threshold
could violate the asymptotic nature of the GPD, and too high of a threshold
would provide only a few relevant data points and lead to an unstable fit. In
order to determine an appropriate threshold, one leverages the GPD-based
expectation value of the data (Eq. 4.8, Coles et al., 2001). For a given
threshold, $u_{0}$, the expectation value of the data $X$ that exceed the
threshold can be written for $\xi$ ¡ 1 as,

$E(X-u_{0}|X>u_{0})=\frac{\sigma_{u_{0}}}{1-\xi}.$

It can be shown that for all thresholds $u>u_{0}$, this expectation value is
(Eq. 4.9, Coles et al., 2001),

$E(X-u|X>u)=\frac{\sigma_{u_{0}}+u\xi}{1-\xi}.$ (8)

That is, for $u>u_{0}$, $E(X-u|X>u)$ is a linear function of $u$. This
expected linearity of 8 with $u$—for values of $u$ for which the GPD fit is
appropriate—can be used to test the validity of the chosen threshold used in
the GPD analysis via a plot which we use in our results section known as the
Mean Residual Life (MRL) plot.

###  2.6 Data

We briefly review the observations used and how they were processed for each
of the three case studies reconsidered here. In addition, we describe the
model Large Ensemble and preindustrial simulations that we analyze.

#### Siberian Heatwave.

This study relies on two sets of observations: (1) daily maximum temperature
extremes for each year from Verkhoyansk as discussed in CC2021; and (2) 4-year
smoothed GMST anomalies from GISTEMP Lenssen et al. (2019). CC2021 also
performed a GEV analysis on temperature anomalies over a region in Siberia,
which we do not focus on here in order to study specifically the extreme
values detected at the Verkhoyansk station. At the time of the CC2021 study,
the daily maximum temperature extreme in 2020 at the Verkhoyansk station was
recorded as 38 ∘C. This temperature extreme currently appears as a missing
data point on the National Climatic Data Center website
(https://www.ncei.noaa.gov/cdo-web/search, station RSM00024266, VERHOJANSK,
data downloaded Aug 2023) but we still use this confirmed record in our
analysis. We also study the GMST anomalies and extreme temperatures over the
Siberian region defined in CC2021 within the CESM1 Large Ensemble Kay et al.
(2015) over the period 1926–2019. Similarly, we study analogous variables
within a 1,800-year CESM1 preindustrial simulation, yielding 18 segments of 94
years for consistency with the observations.

#### Australian Bushfire.

This case study, following OK2021, again uses two sets of observations for the
years 1920–2019: (1) annual (July–June) maxima of 7-day moving average daily-
maximum surface temperature data from the Australian Water Availability
Project (AWAP) over a region in Australia defined in OK2021; and (2) 4-year
smoothed GMST anomalies from GISTEMP. We again apply the same analysis to the
CESM-LE and to the long CESM control run, in this case using 18 segments of 98
years, for consistency with the observations.

#### Madagascar Drought.

Here we use, following HW2022, (1) A 24-month running mean of precipitation
from ERA5 Hersbach et al. (2020) over the region of Madagascar that was
defined in HW2022; and (2) 4-year smoothed GMST anomalies from GISTEMP. The
data range from 1951–2020, and the bottom 20% precipitation from that period
corresponds to the following 14 two-year periods: 1956–1957, 1957–1958,
1958–1959, 1959–1960, 1962–1963, 1990–1991, 1991–1992, 1992–1993, 2008–2009,
2009–2010, 2015–2016, 2016–2017, 2017–2018, 2019–2020. That is, the analysis
following HW2022 relies on a fit to only 14 data points. We discuss the issues
that are involved in fitting extreme distribution functions to such a small
number of data points below. We also study the equivalent variables within the
CESM LE over the period 1950–2019, consistent with the period of observations.
We again apply the same analysis to the CESM-LE and to the long CESM control
run, in this case using 18 segments of 70 years, for consistency with the
observations.

We note that in each of these three case studies, the WWAP leverages many
datasets, including ERA5 reanalysis, station observations, and the multimodel
CMIP5 ensemble. We use only a subset of the data they study and also the CESM
Large Ensemble data for self-consistency in comparing the three examples that
we analyze here.

##  3 Results

Our analysis of the attribution of extreme events based on a fit to extreme
value distributions is presented as follows. First, in subsection 33.1, we
assess the null hypothesis that a GMST dependence of the distribution
parameters reflects an influence of anthropogenic climate change. We do so by
analyzing a 1,800-year CESM control simulation (i.e., without anthropogenic
forcing). In the following three sections (33.2, 33.3, 33.4) we then consider
each of the three test cases and first assess the need for a GMST dependence
and for using a GEV or GPD fit to the observations relative to that of a
normal or log-normal distribution. We then evaluate whether there is a
statistically significant GMST dependence using the deviance statistic.
Finally, we estimate the uncertainty level of the parameter $\alpha$ used for
the attribution and evaluate the GEV or GPD fit to the CESM Large Ensemble in
order to deduce how much data would be necessary to obtain a significant GMST
dependence.

###  3.1 Extremes analysis of simulations with no anthropogenic forcing

The approach of the WWAP attribution methodology is to examine the dependence
of the distribution parameters of extreme events on GMST and, if these
parameters significantly differ from zero, to infer that anthropogenic climate
change affects these extreme events. The implicit null hypothesis is one of
zero dependence between extreme events and GMST in the absence of
anthropogenic influences on climate. Furthermore, the parameter representing
the dependency on GMST is then used to quantify how much a given extreme event
was made more likely as a consequence of climate change. Our objective here is
to point out that it is possible, however, for internal climate variability
such as ENSO to produce covariance between GMST and extremes even without
anthropogenic influences.

To examine this potential role of internal variability, we start with the
Australian Bushfire case and perform a GEV analysis using a 1,800-year CESM
control run, using 18 data segments of length 94 years, equal to the length of
the observed record and analogous to having 18 simulations of historical
climate with no anthropogenic emissions. The distribution of $\alpha$ values
for the GEV control analysis are shown by the green bars in Fig. 1b. Despite
no external forcing being applied to the control run (aside from a constant
volcanic forcing and an 11-year time-varying solar component), we find a large
and positive $\alpha$ calculated from essentially all model output segments
centered around a mean value of 5.1. The large $\alpha$ value means that for
every 1 ∘C GMST warming, extreme summer temperatures increase by 5.1 ∘C. It is
difficult to imagine a physical causal mechanism that would lead to such an
enormous amplification of the GMST signal in Australian heat waves. Instead,
internal variability, likely ENSO, increases both extreme temperatures over
Australia and, to a much smaller degree, affects GMST, creating this seeming
strong dependence between extreme events and GMST. It is indeed known that
several modes of internal variability affect extreme temperatures in
Australia, including ENSO, the Indian Ocean Dipole, and the Southern Annular
Mode Hendon et al. (2007); Ummenhofer et al. (2009); King et al. (2020). The
fact that the covariance between extremes and GMST found in the CESM Large
Ensemble analysis results from a third factor affecting both Australian
temperatures and global mean temperatures highlights a major difficulty with
analyses based on a fit of distribution parameters to GMST.

A similar analysis for Siberia (green bars in Fig. 1a) shows that the value
for $\alpha$ deduced from a record of a similar length to that of the
observations ranges between 1 to 4, although there is again no CO2 increase in
the corresponding model run. This large and variable response shows that a
non-zero $\alpha$ does not necessarily imply an effect of anthropogenic
climate change. The base hypothesis that a nonzero GMST dependence implies an
anthropogenic signal is contradicted by the above analysis. A more realistic
null hypothesis for the parameter dependence on GMST may involve a
distribution of values for $\alpha$.

For the Madagascar Drought example (Fig. 1c), we find a large spread in both
the pre-industrial control and CESM Large Ensemble, centered around zero in
both cases. The fact that the pre-industrial control $\alpha$ distribution is
centered around zero suggests that a null hypothesis of $\alpha=0$ could be
appropriate in this case.

Analysis of greenhouse-forced simulations raises a further issue (Figs. 3c,
7c, SI-4b) wherein the distribution of $\alpha$ in forced simulations is
narrower and closer to zero than in the unforced simulations (Fig. 1),
particularly in the case of the Australian Bushfire. This can be understood in
that internal variability, at least in these examples, evidently projects more
strongly onto regional extremes than GMST as compared to anthropogenic forcing
and can lead to a range of GMST dependencies. Thus, in simulations with
anthropogenic forcing, the expected value of $\alpha$ is smaller, at least in
the three examples that we consider. The implication is that the magnitude of
$\alpha$ may not be a good indicator of anthropogenic influence.

This result calls into question the validity of the empirical tests that WWAP
has used in at least the three studies that we review. It is possible, of
course, that a strong anthropogenic signal in a few decades could
significantly affect extreme events and overwhelm the effects of natural
variability in cases such as the Siberian Heatwave. This attribution would be
more difficult in cases like the Australian example, given the large $\alpha$
detected in the pre-industrial control simulation. It is, therefore, necessary
to study when and how such an emergence of the ACC signal occurs for this
attribution analysis to be used with confidence. We proceed, regardless, to
consider other uncertainties involved in the attribution analysis based on a
fit of extreme value distributions.

Figure 1: Effect of natural climate variability on $\alpha$ distribution. PDFs
$\times 100$ of $\alpha$ for (a) the Siberian Heatwave, (b) the Australian
Bushfire, and (c) the Madagascar Drought attribution cases based on the
1,800-year CESM pre-industrial control run with 94-, 98-, and 70-year
segments, respectively (green) and the 35 Large Ensemble members forced using
the anthropogenic greenhouse scenario (blue). The lengths of these segments
correspond to the extent of the observed records used in CC2021, OK2021, and
HW2022. The mean values for the pre-industrial control and Large Ensemble for
each of the cases are indicated in the boxes.

###  3.2 Siberian Heatwave

We now explore other facets of the attribution of the Siberian Heatwave in the
summer of 2020 to anthropogenic change. This analysis is presented in addition
to the concerns raised by the results of the control runs in the previous
section, indicating that the null hypothesis that $\alpha\neq 0$ implies an
anthropogenic effect on extreme Siberian heat waves is not justified.

Notwithstanding the foregoing finding that the interpretation of $\alpha$ is
complicated, empirical inference of whether an extreme event is attributable
to anthropogenic climate change remains a potentially important line of
inquiry. We, thus, explore other features of the WWAP approach in the context
of their analysis of the Siberian heatwave. We focus specifically on their
analysis of the June daily maximum temperature because this is the longest
record that is available and on the question of whether the data justify the
use of non-stationary (GMST-dependent) extreme value distributions.

While it is understood based on theoretical considerations that a GEV is the
appropriate representation of block-maxima of data Coles et al. (2001), it is
not obvious for a given problem what minimum block size justifies using a GEV,
and specifically if a year of daily maximum data suffices. This is especially
an issue given that the temperature is auto-correlated and, therefore, more
data are needed per averaging block to justify the use of a GEV. It is,
therefore, worthwhile checking how well a more standard (say normal)
distribution does in fitting the data.

A cumulative distribution estimate of the highest June daily-maximum
temperature observed each year between 1926–2020 at the Verkhoyansk station is
shown in Fig. 2a along with a 95% confidence interval of the CDF obtained by
bootstrapping the sample. In the first step of our analysis, we consider
cumulative distribution function fits to the observed distribution assuming
the data follows either a generalized extreme value or normal distribution.
The fitted distributions are specified here to be stationary by setting
$\alpha=0$ and $\mu=\mu_{0}$ in Eqns. 1–3), fits are obtained by maximizing
the likelihood of the data, and the results are shown in Fig. 2a. Both
distributions generally fit the observations well, with the normal fit
completely within the 95% range, and with the GEV fit being good except at the
highest temperature, where it slightly underestimates the probability of
obtaining the highest temperatures. There is only one data point above 35 °C,
so the empirical CDF is poorly constrained there in any case. Interestingly,
the normal distribution and GEV CDFs look similar and lie within the
confidence intervals defined for the empirical CDF, despite the normal
distribution having fewer parameters than the GEV.

The quality of fits seems to suggest that an extreme value distribution may
not be needed as a normal distribution does just as well. Similarly, a non-
stationary distribution with $\alpha\neq 0$ may not be necessary given that
the stationary distributions seem to provide a satisfactory fit to the
empirical cumulative distributions, and we further examine this next.

Figure 2: Siberian Heatwave: (a) Testing the need for a GEV fit. The red curve
indicates the empirical CDF of the highest daily maximum temperature for each
year for the Verkhoyansk station data. The purple and blue lines denote fits
by a GEV (with a constant location parameter) and by a normal distribution,
respectively. The light and dark green shadings correspond to 90th and 95th
percentile confidence intervals (section 22.3). (b) Justification for using a
GMST-dependent mean. Box plots of the deviance statistic between a GEV
distribution without a GMST-dependent mean relative to a GEV distribution with
a GMST-dependent mean using the Siberian heatwave data with 5,000 bootstrap
resamples. The edges of the boxes correspond to the first and third quartiles.
The deviance statistic for the station observations is shown in green as the
leftmost box. The deviance statistic was also computed for each of the 35
members from the CESM Large Ensemble (blue). The 95th percentile significance
level from a $\chi_{k=1}^{2}$ distribution with one degree of freedom is shown
by the horizontal red line, corresponding to the addition of one degree of
freedom by including $\alpha$ in the fitted model.

We now proceed to examine the role of GMST in the fit to a GEV distribution
(eq. 3) in this particular application, following CC2021. The value of
$\alpha$ and its uncertainty range for the Siberian Heatwave do not seem to be
explicitly discussed by CC2021. A linear regression of GMST anomalies against
daily max temperatures gives a slope of 0.84 ∘C per ∘C (Fig. 3a), where the
positive value indicates that the annual maximum daily temperatures at
Verkhoyansk scale with the GMST. The squared Pearson correlation, however, is
only $r^{2}=0.01$, and the p-value from a two-sided test is 0.3, indicating
that GMST is only a weakly significant indicator of trends in Siberian
temperature maxima.

A maximum likelihood fit of a nonstationary GEV, including $\alpha$, gives
results that are consistent with the simple linear regression between GMST and
temperature maxima. Specifically, we find that the $\alpha$ without resampling
in the observations is 1.5 ∘C per ∘C, similar to the slope of 1.1 obtained
from the linear regression. The 95% confidence interval is [$-1.2$, 4.8] (Fig.
3b), indicating that $\alpha$ is not significantly different from zero. Note
that we use a two-sided test here because we would also have evaluated a
scenario where the maximum likelihood value of $\alpha$ was less than zero.
The deviance statistics (section 22.4), $D$, indicate that for the
observations, the addition of one more model parameter, $\alpha$, representing
the GMST dependence, is justified as the value is above the threshold
indicated by the red horizontal line (Fig. 2a, leftmost bar), yet the
uncertainty estimate include a large range of non-significant values of the
test. The solution for all GEV parameters for the observations is given in
Supplementary Fig. SI-1.

Figure 3: Siberian Heatwave: uncertainty in $\alpha$. (a) Scatter plot of the
daily max temperature at the Verkhoyansk station against GMST. The red line
indicates a linear regression fit to the data, and the blue line indicates a
linear fit using the location parameter and $\alpha$ derived from the GEV
maximum likelihood estimate. (b) PDF $\times 100$ of the GEV-fit $\alpha$
distributions from the station data, using 5,000 bootstrap resamplings of the
observed time series. The $\alpha$ value from the observations (no resampling)
is shown by the vertical red line. (c) Shown is a histogram of the GEV-fit
$\alpha$ PDFs $\times 100$ derived from the CESM Large Ensemble, where the
$\alpha$ values are calculated from each ensemble member with no resampling
(blue). The average $\alpha$ value across the Large Ensemble is indicated with
the blue line. There are 94 annual data points for each ensemble member.

To provide context for the interpretation of the observational results, we
also examine the greenhouse-forced Large Ensemble output. We find that all
ensemble members within the CESM LE indicate a positive value of $\alpha$,
with a mean value of 2.60 ∘C per ∘C (Fig. 3c). Similarly, 33 of 35 ensemble
members give a deviance statistic that supports the use of a non-stationary
component (that is, a GMST-dependence) in the GEV statistic (Fig. 2b).
Unsurprisingly, if multiple ensemble members are analyzed simultaneously,
which mimics the availability of more data, the distribution of the GEV
parameters narrows, generally converging around the mean value of the
parameters across individual ensemble member fits (Fig. 4). The foremost
reason that the Large Ensemble suggests a stronger link between Siberian
extremes and GMST may be that the LE, on average, simulates an Arctic-wide
summer warming rate of 1.80 ∘C per ∘C, whereas GISTEMP observations only
indicate a rate of 1.3 ∘C per ∘C (Fig. 5a,b). One expects heat waves to be
strongly influenced by the JJA climatology, and Fig. 5c shows that different
ensemble members, corresponding to different realizations of natural
variability, show very different dependencies of the Arctic JJA climatology on
GMST. It is difficult, therefore, to estimate from the single realization
corresponding to the observed record to what degree the Arctic JJA climatology
is affected by anthropogenic climate change versus natural variability, adding
to the uncertainty in attributing Siberian heat waves.

The overall result for the GMST dependence for the Siberian Heatwave analysis,
which is the base of the attribution to anthropogenic climate change, seems
superficially consistent with the results found from CC2021. However, the
important caveats are that (1) we found in the previous subsection that
internal variability can lead to GMST dependence that invalidates the null
hypothesis (Fig. 1a); (2) a normal distribution is able to fit the observed
data well without the need for a GEV and without GMST-dependence of the
distribution parameters (Fig. 2a). While a GEV fit may, a priori, make sense
to fit to annual daily temperature extremes, due to both the relatively short
observational record and also the autocorrelated nature of the temperature
extremes, a normal distribution does not appear to be ruled out, and we find
that it does not perform meaningfully worse than the GEV. (3) the uncertainty
range for $\alpha$ includes negative values (Fig. 3b).

Figure 4: Siberian Heatwave: uncertainty in GEV parameters as a function of
the record length. PDFs (times 100) of the GEV parameters as a function of the
number of ensemble members used. The data were derived from the highest daily
maximum temperature in June, as in CC2021, using the CESM Large Ensemble.
Contours are drawn every 0.3 for (a), 2 for (b), 1.5 for (c), and 0.7 for (d),
all starting at 0.05. We denote the 95% confidence range for each number of
ensemble members by the red asterisks. Figure 5: Natural variability and the
attribution of Siberian heatwaves. (a) A scatter plot of the annual average
JJA Arctic (north of $70^{\circ}$ N) surface temperatures against the 4-year
running mean of GISTEMP annual GMST. (b) Same, for the 35 members of the CESM
Large Ensemble over the period 1926–2019. (c) PDF $\times 100$ of the slopes
of annual average JJA Arctic surface temperatures against the 4-year running
mean of annual GMST within each member of the CESM Large Ensemble with 5,000
bootstrap resamples (red). The same distribution is indicated for GISTEMP
observations in blue. The vertical red line indicates the average slope across
the CESM Large Ensemble, and the vertical blue line indicates the slope
determined from observations.

###  3.3 Australian Bushfire

We now consider the attribution analysis of the Australian Bushfires of the
summer of 2020–2021, again using both observations and the CESM Large
Ensemble, following the observational analysis of OK2021. While Fig. 1b showed
a very large value of $\alpha$ for Australian heat waves in a control run with
a fixed CO2, invalidating the basic assumption that a nonzero GMST dependence
is an indication of anthropogenic global warming, we proceed regardless. We
again find that a normal distribution and a GEV distribution (both with a
constant mean/location parameter) fit the 7-day moving mean of daily
temperature maxima in the observations without a GMST dependence (Fig. 6a). In
this case, the fit is exceptionally good, seemingly invalidating the need for
a GEV or for a GMST dependence. The deviance statistic/likelihood ratio test
for the model LE shows that only five of the 35 ensemble members justify
adding GMST dependence to the distribution parameters (Fig. 6b). This raises
the question of whether the influence of GMST, and therefore of anthropogenic
climate change in the WWAP formulations, is a necessary part of the
statistical description of extreme temperatures in Australia as represented in
the model large ensemble. The deviance statistic for the observations does
suggest a significant GMST component (leftmost bar in Fig. 6b), but given the
possibility that this component is driven by natural variability rather than
by anthropogenic change (Fig. 1a), one must conclude again that even if
$\alpha$ is justified statistically, it might be best not to use it for the
purpose of attribution analysis.

Figure 6: Australian Bushfire: (a) Testing the justification for a GEV. The
red curve indicates the empirical CDF of the highest annual 7-day running mean
of daily maximum temperature over the area defined in Van Oldenborgh et al.
(2021) for observations. The purple and blue lines denote GEV (with a constant
mean), and normal distribution fits, respectively. The light and dark green
shadings correspond to 90th and 95th percentile confidence intervals. (b)
Justification for using a GMST-dependent mean. Box plots of the deviance
statistic between a GEV without a GMST component and a GEV with a GMST
component using the Australia Bushfire example with 5,000 bootstrap resamples.
The edges of the boxes correspond to the first and third quartiles. The
deviance statistic was computed for each of the 35 members from the CESM Large
Ensemble (blue) and compared with the 95th percentile significance level from
a $\chi^{2}$ distribution with one degree of freedom (red line). A box plot
for the deviance statistic for observations is shown in green.

The regression analysis of extreme temperatures and GMST, (Fig. 7a) as well as
the range of $\alpha$ values estimated from observations (Fig. 7b) both show a
GMST dependence of 2.0 ∘C change in extreme temperatures per degree increase
in GMST. As mentioned above, a sensible null hypothesis for heat waves in a
warming climate is that the extreme values shift with the mean warming. The
GMST dependence calculated here suggests that Australian heat waves increase
at nearly twice the rate of the GMST. Such a stronger response requires a
physical mechanism which we are unable to come up with in this case. The
alternative is that this is a result of natural variability (likely ENSO)
affecting both Australian heat extremes and the GMST, consistent with Fig. 1.
Again, this suggests that using this value of $\alpha$ for attribution may
lead to unreliable results. Finally, Fig. 7c shows a large range of GMST-
dependencies for different ensemble members in the CESM-LE, again highlighting
that a large component of this dependence most likely represents natural
variability rather than anthropogenic change.

Figure 7: Australian Bushfire: distribution of $\alpha$ values. (a) Scatter
plot of the highest annual 7-day running mean of daily maximum temperatures
against GMST (following Van Oldenborgh et al., 2021). The red line indicates a
linear regression fit to the data, and the blue line indicates a linear fit
using the location parameter and $\alpha$ derived from the GEV maximum
likelihood estimate. (b) PDF of the GEV-fit $\alpha$ distributions from
observations, using bootstrap resampling. The $\alpha$ value from the
observations without resampling is shown by the red line (at 2.0) (c) Effect
of the natural variability and number of data from a model large ensemble on
the uncertainty in $\alpha$ for the Australian Bushfire example. Shown is a
histogram of the GEV-fit $\alpha$ PDFs derived from the $\alpha$ value from
each ensemble member with no resampling (blue). The average $\alpha$ value
across the Large Ensemble is indicated with the blue line. The PDFs for each
ensemble member calculated by bootstrapping are shown by the gray curves.

###  3.4 Madagascar Drought

Finally, we consider the attribution analysis for the 2019–2021 drought in
Madagascar, following the lead of HW2022. This case differs from the two
examined above in that it analyzes precipitation extremes that are positive by
definition. HW2022, therefore, used a different extreme value distribution
function, the GPD, which is meant to represent the high tail of the
distribution of extreme events over a specified threshold Coles et al. (2001).
As discussed in Section 2, HW2022 analyzed the driest 20% of the 2-year
averages of precipitation data from 1951–2020 multiplied by minus one to turn
minima (droughts) into maxima. We follow their example, with the important
caveat made in Section 22.2 that the validity of using a GPD this way for
precipitation data cannot be taken for granted given that there is a hard
maximum of zero for these data of negative precipitation rates, a hard maximum
that seems inconsistent with the definition of the GPD. We also note that
taking the lowest 20% of 2-year precipitation data points amounts to only 14
data points, likely too low a number for a reliable attribution.

We find again that an extreme value distribution is not required by the data
distribution, and a GMST dependence is not needed: the GPD fit to the bottom
20th percentile precipitation data is no better than a log-normal distribution
fitted to all of the precipitation data, both with constant distribution
parameters that do not depend on the GMST (Fig 8a). A log-normal distribution
seems to be a reasonable null hypothesis for the positive precipitation data.
In fact, the empirical CDF is practically unconstrained at its lower tail
because there are no data points below a value of 1.25 mm/day.

Figure 8: Madagascar Drought: (a) Evaluating the need for a GPD fit. The
empirical CDF is shown (red line, with 90th and 95th percentile confidence
intervals shown by red shadings), based on 24-month running means of daily
averaged annual precipitation over the Madagascar area defined in HW2022 using
ERA5 precipitation. A GPD fit to the bottom 20th percentile precipitation
(black) and a log-normal fit to the full distribution (blue). (b) The deviance
statistic between a GPD with and without a GMST component in the CESM Large
Ensemble (blue) and observations (green). (c) Scatter plot of the 2-year
running mean of ERA5 daily precipitation over the Madagascar area defined in
HW2022 against the 4-year running mean of annual GMST for the lowest 20
percent precipitation data. The dashed blue line shows the exponential GPD
log-likelihood fit for the location parameter for the observations without
resampling (eqn. 6) following HW2022. (d) MRL plot (section 22.5) for
different precipitation thresholds. The blue shading defines a region of
confidence calculated as $\pm$ one standard deviation. The red line denotes
the threshold used in the GPD fit of HW2022 (i.e., the 20th percentile). The
orange line represents the slope of Equation 8 ($\xi/(1-\xi)$) with the line
shifted to have the same value as the MRL line at its smallest value. The
green line represents the full linear equation from Equation 8
($\sigma+u\xi/(1-\xi)$).

While the deviance statistic (likelihood ratio) for the observations is
marginally above the red significance line (leftmost green bar, Fig. 8b), none
of the 35 ensemble members exceed the significance level. Overall, Figs. 8a,b
therefore suggest that either the CESM Large Ensemble does not adequately
capture the relevant features, or internal variability plays a dominant role
in extreme precipitation in Madagascar and the significant anthropogenic role
in the observed data is a spurious correlation with the particular sequence of
internal variability that took place historically.

As a reminder, the location parameter of the GPD, corresponding to the maximum
precipitation value considered as an extreme drought year in this attribution
case (or the minimum negative precipitation), was assumed by HW2022 to vary
exponentially with the GMST ($\overline{T}$), as
$\mu=\mu_{0}\exp(\alpha\overline{T}/\mu_{0})$ with $\mu_{0}<0$. If droughts
become more severe with increasing GMST, we expect the location parameter for
the distribution of the negative precipitation to get larger (less negative)
with increasing GMST, and therefore $\alpha/\mu_{0}$ to be negative and thus
$\alpha$ to be positive. A scatter plot of ERA5 precipitation against GMST
anomalies is shown in Fig. 8c. The dashed blue line indicated in the plot
represents an exponential fit using the $\alpha$ and $\mu_{0}$ coefficients
from the maximum log-likelihood GPD fit. This fit looks notably linear,
consistent with the Taylor expansion of the exponential function, as below,

$\displaystyle\mu=\mu_{0}\exp(\alpha T/\mu_{0})\approx\mu_{0}(1+\alpha
T/\mu_{0})=\mu_{0}+\alpha T,$

where the parameter $\alpha$ is used to fit the negative of the precipitation.
The value of $-{\alpha}_{\mathrm{GPD}}$ calculated from observations is small
and indistinguishable from zero (Fig. 8c).

Fig. 8d shows the MRL plot (section 22.5) for the Madagascar Drought case. As
a reminder, this plot needs to show a linear ranger within which the threshold
value for the GPD analysis is chosen (section 22.5). For small and large
values of precipitation, the MRL curve appears to curve nonlinearly,
indicating that a GPD fit is likely invalid using those thresholds. For
precipitation thresholds ranging from roughly 1.6 to 2 mm/day, the MRL curve
plateaus. The orange line shows the expected MRL slope, $\xi/(1-\xi)$; that
slope does not alight with the plateaued curve around the value of 1.9
corresponding to the 20th percentile precipitation threshold chosen in
Harrington et al. (2022). In general, it is difficult to identify the required
linear range in this plot that justifies using a GPD distribution and allows
selecting a threshold.

The small value of $\alpha$ found for observations suggests that the increase
in GMST has so far hardly affected the statistics of droughts in Madagascar.
While this is also the conclusion reached by HW2022, we found above that the
use of a GMST-dependence is not justified by the data in any case (Fig. 8a,b),
a GPD fit is not justified (Fig. 8a,d), and the number of data points involved
in the fit (14) seems too small to draw any meaningful conclusions on an
extreme value distribution as it leads to a large uncertainty interval shown
by the shading. The fit, including the offset (green line), does not
approximate the expected value curve. Supplemental Fig. SI-6a indicates that
even with 35 ensemble members’ worth of data, the uncertainty range for
$\alpha$ includes both positive and negative values. One could, therefore,
conclude that the approach used here is not likely to lead to confident
attribution in this case in the next decades and beyond.

##  4 Conclusions

The detection and attribution of extreme weather events represent an important
goal of research in order to understand the present and future impacts of
anthropogenic climate change. We were inspired in this work by a series of
attribution studies by the “World Weather Attribution Project” (WWAP) and
evaluate, in particular, the empirical methodology used by the WWAP for
purposes of attribution. Our goal is to evaluate the methodology and offer
some additional tests that can increase confidence in the meaningfulness of
the attribution results. The WWAP method involves fitting an extreme value
distribution function to an observed record, with the distribution parameters
depending on the global mean surface temperature (GMST). For example, the
location parameter of an extreme value distribution (GEV) is made a linear
function of GMST, $\overline{T}$, as $\mu=\mu_{0}+\alpha\times\overline{T}$.
The GMST dependence $\alpha$ is assumed to represent the effects of
anthropogenic climate change and is then used to calculate the change in
return time of extreme events and the degree to which the event was made more
or less likely due to anthropogenic global warming. This method was applied to
many extreme events such as heat waves, cold spells, floods, droughts, and
more (e.g., van Oldenborgh et al., 2016; van der Wiel et al., 2017; van
Oldenborgh et al., 2017; Ciavarella et al., 2021; Van Oldenborgh et al., 2021;
Harrington et al., 2022). We used observations, a CESM preindustrial
simulation, and the CESM Large Ensemble model output and re-examined three
specific recent attribution cases using this methodology as example
applications of our offered refinements.

We begin by examining the GMST dependence of extreme heatwave events in
Siberia (following CC2021) and Australia (following OK2021) and Madagascar
(following OK2021), using data from a CESM1 control simulation where CO2 is
kept at a preindustrial level. We find a significant GMST dependence in the
first two cases, indicating that such dependence cannot be used to
unequivocally identify the effects of anthropogenic climate change. The GMST
dependence, in these cases, is a result of regional natural variability such
as ENSO, which affects the extreme events being examined, and also weakly
affects the GMST. This seems to invalidate the hypothesis that a non-zero
$\alpha$ is a direct result of, and can be used to quantify, the effect of
anthropogenic climate change due to the increase in greenhouse forcing. The
value of $\alpha$ in the Australian case is especially large, seemingly
indicating that every degree increase in GMST leads to 5 degrees increase in
Australian heat waves. Instead, we infer that the large values of $\alpha$ are
a result of ENSO affecting GMST by 0.2 degrees for every 1-degree effect it
has on Australian heat waves.

We proceed, regardless, to examine three specific attribution cases in order
to examine other aspects of the proposed attribution methodology. We begin the
analysis of each attribution case by testing if an extreme value distribution
is indeed needed. We plot the empirical cumulative distribution function (CDF)
with bootstrapping error bars (Figs. 2a, 6a, 8a) and superimpose the fits of
an extreme value distribution (GEV or GPD) without allowing a GMST dependence
of the distribution parameters, as well as a standard (normal or a log-normal)
distribution. If the standard distribution falls within the error bars of the
empirical distribution, this serves as one indication that there may not be a
need for an extreme value distribution. If the normal or extreme distributions
perform well without a GMST dependence, this indicates that such dependence
may not be justified by the data.

As the next step, we further test if the addition of a GMST dependence of the
distribution parameters is justified statistically. The distribution
parameters are estimated by the WWAP via a fit to the data using log-
likelihood maximization. When adding a parameter (specifically $\alpha$), one
expects the fit to improve, but the question is whether the improvement is
beyond what is expected simply due to the addition of one more model
parameter, so that one is convinced that the additional parameter represents a
genuine effect in the data. The importance of the additional model parameter
($\alpha$) can be tested using the deviance statistic (aka log-likelihood
ratio test). This test provides critical information given that the entire
attribution effort in this approach is based on the estimated value of
$\alpha$. Finally, the GPD distribution is meant for high values above a
threshold, and the expected mean needs to vary linearly with the threshold
chosen. We use this constraint via a mean residual life (MRL) plot as a
consistency test to see if a GPD is justified and if the threshold is chosen
appropriately.

We find for the Siberian Heatwave of 2020 (CC2021) that (1) both a normal
distribution and a GEV fit the data successfully without a GMST-dependent
mean; (2) the deviance statistic (i.e., the likelihood ratio test) suggests
that making the distribution location parameter $\mu$ dependent on GMST via
$\alpha$ is justified for observations and for the large model ensemble,
although one worries about the previous finding that alpha is non zero also
without the effects of ACC.

For the high regional temperatures that preceded the Australian bushfires of
2020–2021 (OK2021), a GEV is not needed based on a test fit to an empirical
CDF. The deviance statistic, though, seems to suggest a clear GMST dependence
(Fig. 6b). This seeming dependence is also relatively strong in magnitude: a 1
∘C increase in GMST would correspond to a 2.0 ∘C increase in temperature
extremes. Such amplification of the anthropogenic warming signal over
Australia requires a physical mechanism we are not aware of. The most likely
explanation is that the signal is a result of ENSO affecting Australia as well
as the GMST, as found for the CESM control runs discussed above.

The Madagascar drought attribution (HW2022) was based on only 14 2-year
averaged data points of the lowest precipitation rates, hardly sufficient to
construct a meaningful empirical CDF or to constrain the four parameters of
the GPD distribution used for that purpose. Similarly, the MRL plot does not
demonstrate that these drought events are well represented by a GPD. While our
results—that attribution of the Madagascar drought using this methodology is
not possible—agree with those of HW2022, we emphasize our different
perspective that it is unlikely that attribution would be possible based on
such a small number of data points, as well as the other difficulties pointed
out above. These results do not rule out, however, that using additional data
sources and a different approach could lead to a successful attribution of the
Madagascar Drought.

An important issue raised by this work is the discrepancy between the modeled
and observed findings, as seen, for example, in the deviance statistic results
for the observations versus model large ensemble (Figs. 2b, 6b, 8b). This may
be due to a particularly favorable sequence of internal variability in the
observations, measurement error in the observational record, or issues with
the simulation of extreme weather in the CESM Large Ensemble. It would be
helpful to test if these discrepancies also occur with other large model
ensembles.

Overall, our results suggest uncertainty in the interpretation of the
dependence of the parameters of extreme event distributions on the global mean
surface temperature as reflecting anthropogenic climate change. The proposed
deviance statistic, MRL plot, and a careful analysis of uncertainty ranges
demonstrated here may help avoid misguided confidence in attribution results.
As noted in the introduction, WWAP does not exclusively rely upon the
empirical methodology evaluated here in order to reach conclusions, and also
evaluates other historical data, model simulations, and the dynamical context.
Nevertheless, the WWAP empirical methodology examined here figures prominently
in their analyses, and we have highlighted several concerns with its
interpretation, foremost that internal variability can give rise to the
appearance of strong covariance between GMST and local extremes. Our results
do not provide evidence against a relationship between GMST and local
extremes, but rather that internal variability needs to be taken into account
in evaluating the statistical significance of any such a relationship. It,
therefore, seems that it may be best to use the WWAP empirical methodology in
conjunction with other attribution approaches, which are summarized in the
Introduction (e.g., Stott et al., 2016; Knutson et al., 2017; Seneviratne et
al., 2021).

###### Acknowledgements.

PS is supported by the Harvard Global Institute. ET is supported by DOE grant
DE-SC0023134, and thanks the Weizmann Institute for its hospitality during
parts of this work. PH is supported by NSF Award 2123295. The CESM project is
supported primarily by the National Science Foundation (NSF). We acknowledge
the CESM Large Ensemble Community Project and supercomputing resources
provided by NSF/CISL/Yellowstone. We also acknowledge instructive
conversations with Angela Rigden. All data used in this paper are available in
public repositories. The climatological output from the CESM-LE and CESM
control runs used in this paper are publicly available using the Intake-ESM
Python package. All codes and data used here are available under the public
open science foundation depository https://osf.io/fycu2/, DOI
10.17605/OSF.IO/FYCU2.

## References

  * Allen (2003) Allen, M., 2003: Liability for climate change. Nature, 421 (6926), 891–892.
  * Christidis et al. (2015) Christidis, N., P. A. Stott, and F. W. Zwiers, 2015: Fast-track attribution assessments based on pre-computed estimates of changes in the odds of warm extremes. Climate Dynamics, 45 (5), 1547–1564.
  * Ciavarella et al. (2021) Ciavarella, A., and Coauthors, 2021: Prolonged Siberian heat of 2020 almost impossible without human influence. Climatic Change, 166 (1), 1–18.
  * Coles et al. (2001) Coles, S., J. Bawa, L. Trenner, and P. Dorazio, 2001: An introduction to statistical modeling of extreme values, Vol. 208. Springer.
  * De Paola et al. (2018) De Paola, F., M. Giugni, F. Pugliese, A. Annis, and F. Nardi, 2018: Gev parameter estimation and stationary vs. non-stationary analysis of extreme rainfall in african test cities. Hydrology, 5 (2), 10.3390/hydrology5020028, URL https://www.mdpi.com/2306-5338/5/2/28.
  * Grubb et al. (2022) Grubb, M., and Coauthors, 2022: Introduction and Framing. IPCC, 2022: Climate Change 2022: Mitigation of Climate Change. Contribution of Working Group III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, P. Shukla, J. Skea, R. Slade, A. A. Khourdajie, R. van Diemen, D. McCollum, M. Pathak, S. Some, P. Vyas, R. Fradera, M. Belkacemi, A. Hasija, G. Lisboa, S. Luz, and J. Malley, Eds., Cambridge University Press, Cambridge, UK and New York, NY, USA, chap. 1, 10.1017/9781009157926.003.
  * Harrington et al. (2022) Harrington, L. J., and Coauthors, 2022: Attribution of severe low rainfall in southern madagascar, 2019-21. World Weather Attribution, URL https://www.worldweatherattribution.org/wp-content/uploads/ScientificReport˙Madagascar.pdf.
  * Hendon et al. (2007) Hendon, H. H., D. W. Thompson, and M. C. Wheeler, 2007: Australian rainfall and surface temperature variations associated with the Southern Hemisphere annular mode. Journal of Climate, 20 (11), 2452–2467.
  * Hersbach et al. (2020) Hersbach, H., and Coauthors, 2020: The era5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 146 (730), 1999–2049, https://doi.org/10.1002/qj.3803, URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3803, https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/qj.3803.
  * Kay et al. (2015) Kay, J. E., and Coauthors, 2015: The Community Earth System Model (CESM) large ensemble project: A community resource for studying climate change in the presence of internal climate variability. Bulletin of the American Meteorological Society, 96 (8), 1333–1349.
  * King et al. (2020) King, A. D., A. J. Pitman, B. J. Henley, A. M. Ukkola, and J. R. Brown, 2020: The role of climate variability in Australian drought. Nature Climate Change, 10 (3), 177–179.
  * Knutson et al. (2017) Knutson, T., J. Kossin, C. Mears, J. Perlwitz, and M. Wehner, 2017: Detection and attribution of climate change. Climate Science Special Report: Fourth National Climate Assessment, Volume I, U.S. Global Change Research Program, Washington, DC, USA.
  * Lenssen et al. (2019) Lenssen, N., G. Schmidt, J. Hansen, M. Menne, A. Persin, R. Ruedy, and D. Zyss, 2019: Improvements in the gistemp uncertainty model. J. Geophys. Res. Atmos., 124 (12), 6307–6326, 10.1029/2018JD029522.
  * Lorenz (1969) Lorenz, E. N., 1969: Atmospheric predictability as revealed by naturally occurring analogues. Journal of Atmospheric Sciences, 26 (4), 636–646.
  * Lott et al. (2013) Lott, F. C., N. Christidis, and P. A. Stott, 2013: Can the 2011 east african drought be attributed to human-induced climate change? Geophysical Research Letters, 40 (6), 1177–1181, https://doi.org/10.1002/grl.50235, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/grl.50235.
  * Otto et al. (2020) Otto, F. E. L., and Coauthors, 2020: Challenges to understanding extreme weather changes in lower income countries. Bulletin of the American Meteorological Society, 101 (10), E1851 – E1860, 10.1175/BAMS-D-19-0317.1, URL https://journals.ametsoc.org/view/journals/bams/101/10/bamsD190317.xml.
  * Pascale et al. (2020) Pascale, S., S. B. Kapnick, T. L. Delworth, and W. F. Cooke, 2020: Increasing risk of another cape town ”day zero” drought in the 21st century. Proceedings of the National Academy of Sciences, 117 (47), 29 495–29 503.
  * Pendergrass (2018) Pendergrass, A. G., 2018: What precipitation is extreme? Science, 360 (6393), 1072–1073, 10.1126/science.aat1871, URL https://www.science.org/doi/abs/10.1126/science.aat1871, %**** ̵main.bbl ̵Line ̵125 ̵****https://www.science.org/doi/pdf/10.1126/science.aat1871.
  * Perkins et al. (2007) Perkins, S., A. Pitman, N. J. Holbrook, and J. McAneney, 2007: Evaluation of the ar4 climate models’ simulated daily maximum temperature, minimum temperature, and precipitation over australia using probability density functions. Journal of climate, 20 (17), 4356–4376.
  * Seneviratne et al. (2021) Seneviratne, S., and Coauthors, 2021: Weather and climate extreme events in a changing climate. Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA.
  * Stott et al. (2004) Stott, P. A., D. Stone, and M. R. Allen, 2004: Human contribution to the european heatwave of 2003. Nature, 432 (7017), 610–614.
  * Stott et al. (2016) Stott, P. A., and Coauthors, 2016: Attribution of extreme weather and climate-related events. WIREs Climate Change, 7 (1), 23–41, https://doi.org/10.1002/wcc.380, URL https://wires.onlinelibrary.wiley.com/doi/abs/10.1002/wcc.380, https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcc.380.
  * Trevino et al. (2020) Trevino, A. M., K. A. McKinnon, and P. J. Huybers, 2020: Extremely wrong: When limited data leads to biases in general extreme values distribution estimates. AGU Fall Meeting Abstracts, Vol. 2020, A147–0010.
  * Ummenhofer et al. (2009) Ummenhofer, C. C., M. H. England, P. C. McIntosh, G. A. Meyers, M. J. Pook, J. S. Risbey, A. S. Gupta, and A. S. Taschetto, 2009: What causes southeast Australia’s worst droughts? Geophysical Research Letters, 36 (4).
  * van der Wiel et al. (2017) van der Wiel, K., and Coauthors, 2017: Rapid attribution of the august 2016 flood-inducing extreme precipitation in south louisiana to climate change. Hydrology and Earth System Sciences, 21 (2), 897–921, 10.5194/hess-21-897-2017, URL https://hess.copernicus.org/articles/21/897/2017/.
  * van Oldenborgh et al. (2016) van Oldenborgh, G. J., F. E. L. Otto, K. Haustein, and K. AchutaRao, 2016: The heavy precipitation event of december 2015 in chennai, india. Bulletin of the American Meteorological Society, 97 (12), S87 – S91, 10.1175/BAMS-D-16-0129.1, URL https://journals.ametsoc.org/view/journals/bams/97/12/bams-d-16-0129.1.xml.
  * Van Oldenborgh et al. (2012) Van Oldenborgh, G. J., A. Van Urk, and M. Allen, 2012: The absence of a role of climate change in the 2011 thailand floods. Bull. Amer. Meteor. Soc, 93, 1047–1049.
  * van Oldenborgh et al. (2017) van Oldenborgh, G. J., and Coauthors, 2017: Attribution of extreme rainfall from hurricane harvey, august 2017. Environmental Research Letters, 12 (12), 124 009, 10.1088/1748-9326/aa9ef2, URL https://dx.doi.org/10.1088/1748-9326/aa9ef2.
  * Van Oldenborgh et al. (2021) Van Oldenborgh, G. J., and Coauthors, 2021: Attribution of the Australian bushfire risk to anthropogenic climate change. Natural Hazards and Earth System Sciences, 21 (3), 941–960.
  * Vautard and Yiou (2009) Vautard, R., and P. Yiou, 2009: Control of recent european surface climate change by atmospheric flow. Geophysical Research Letters, 36 (22).
  * Yamaguchi et al. (2020) Yamaguchi, M., J. C. Chan, I.-J. Moon, K. Yoshida, and R. Mizuta, 2020: Global warming changes tropical cyclone translation speed. Nature Communications, 11 (47).
  * Yiou et al. (2007) Yiou, P., R. Vautard, P. Naveau, and C. Cassou, 2007: Inconsistency between atmospheric dynamics and temperatures during the exceptional 2006/2007 fall/winter and recent warming in europe. Geophysical Research Letters, 34 (21).
  * Yoshida et al. (2017) Yoshida, K., M. Sugi, R. Mizuta, H. Murakami, and M. Ishii, 2017: Future changes in tropical cyclone activity in high-resolution large-ensemble simulations. Geophysical Research Letters, 44 (19), 9910–9917, https://doi.org/10.1002/2017GL075058, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL075058.
  * Zwiers et al. (2011) Zwiers, F. W., X. Zhang, and Y. Feng, 2011: Anthropogenic influence on long return period daily temperature extremes at regional scales. Journal of Climate, 24 (3), 881 – 892, 10.1175/2010JCLI3908.1, URL https://journals.ametsoc.org/view/journals/clim/24/3/2010jcli3908.1.xml.

## Supplementary material

Figure SI-1: Siberian Heatwave: distribution of observation GEV values.
Probability distribution functions (times 100) of the GEV parameters using a
5,000-times bootstrap resampling of the annual observational time series.
Figure SI-2: Australia Bushfire: distribution of GEV parameter values
calculated for the observed record. Probability distribution functions (times
100) of the GEV parameters using a 5,000-times bootstrap resampling of the
observed annual time series. The mean values for each parameter are indicated
with the red lines within each panel. The double peak structure of the $\xi$
distribution is purely a result of one data point seen as the highest 7-day
running mean of daily maximum temperatures at GMST $\approx 0$ in Fig. 7a.
Removing that data point leads to a single peak $\xi$ distribution connecting
the two peaks seen here. The distributions of the other GEV parameters are not
significantly affected by the removal of that data point. Figure SI-3:
Australian Bushfire: dependence of confidence in GEV parameters on the number
of ensemble members. Probability distribution functions (times 100) of the GEV
parameters as a function of the number of ensemble members used. The data were
derived from the highest annual 7-day running mean of daily maximum
temperatures, as in Van Oldenborgh et al. (2021), but using the CESM Large
Ensemble. Smaller values over a column indicate a greater spread in the
distribution of the parameter. Contour levels are indicated every 0.3 for (a),
2 for (b), 1.5 for (c), and 0.7 for (d), all starting at 0.1. We denote the
95% confidence range for each number of ensemble members by the red asterisks.
Figure SI-4: Madagascar Drought: distribution of $\alpha$ values. (a) PDF of
the GEV-fit $\alpha$ distributions from observations, using bootstrap
resampling. The mean $\alpha$ value for the observations is indicated with the
red line (at $0.026$). (b) Effect of the number of data on the uncertainty in
$\alpha$ for the Madagascar Drought example. Shown is a histogram of the GEV-
fit $\alpha$ PDFs derived from the $\alpha$ value from each ensemble member
with no resampling (blue). The resampling of the data within each ensemble
member is indicated with the gray curves. Figure SI-5: Madagascar Drought:
distribution of observation GPD values. Probability distribution functions
(times 100) of the GPD parameters using a 5,000-times bootstrap resampling of
the annual observational time series. The mean values for each parameter are
indicated with the red lines within each panel. Figure SI-6: Madagascar
Drought: dependence of confidence in GPD parameters on the number of ensemble
members. Probability distribution functions (times 100) of the GPD parameters
as a function of the number of ensemble members used. Contour levels are
indicated every 0.8 for (a), 0.4 for (b), 1 for (c), and 2 for (d), all
starting at 0.1. We denote the 95% confidence range for each number of
ensemble members by the red asterisks.