The Risk of Negatively Biased and Over-confident Return Level Estimates: A Critique of the Metastatistical Approach to Extremes

1 P R E P R I N T 1 Submitted to Stochastic Environmental Research and Risk Assessment 2 3 The Risk of Negatively Biased and Over-confident Return Level 4 Estimates: A Critique of the Metastatistical Approach to Extremes 5 By Torben Schmith, Karsten Arnbjerg-Nielsen, and Bo Christiansen. 6 National Centre for Climate Research, Danish Meteorological Institute, Copenhagen, Denmark 7 Corresponding author Email: [email protected] 8 ORCID TS: https://orcid.org/0000-0003-3442-4381 , KAN: https://orcid.org/0000-0002-6221-9505 , 9 BC: https://orcid.org/0000-0003-2792-4724 10 11 Abstract 12 Classical extreme value analysis (EVA) often provides large uncertainties on estimated return 13 levels due to limited amount of data available. Marani and Ignaccolo (2015) aim to overcome this 14 by the metastatistical extreme value (MEV) approach. Here extremes are treated as large ordinary 15 events described by one common, known distribution, and therefore a much larger pool of data 16 are available for estimation. They performed Monte Carlo simulations with synthetic Weibull17 distributed rainfall series and showed that the MEV approach gives unbiased estimates of 18 extremes with a smaller uncertainty than classical EVA does. However, the MEV approach neglects 19 that many complex physical mechanisms influence rainfall. This means that the tail behavior of the 20 distribution cannot be inferred from the ordinary events. We therefore replicated their work but 21 added new Monte Carlo experiments to study the classical EVA and the MEV methodologies with a 22 slightly perturbed tail of the underlying distribution. When applying the MEV approach, i.e. fitting 23 a Weibull distribution to the perturbed Weibull series, we obtained systematically negatively 24 biased estimates with too narrow confidence intervals – MEV became over-confident. In contrast, 25 classical EVA also here produced unbiased estimates. Finally, we showed that goodness-of-fit tests 26 are not able to provide guidance on whether MEV can provide unbiased and confident return 27 levels. Further Monte Carlo simulations showed that these conclusions seem to be quite general 28 and not dependent on the specific distribution. Consequently, the MEV approach is unsuitable to 29 provide reliable return levels, and we strongly caution against using it in real-world applications. 30 Keywords 31 ‘extreme value statistics’ , metastatistical, MEV, GEV, Gumbel, bias 32 33 2 1 Introduction 34 Classical extreme value analysis (EVA) (Coles 2001) is widely used within e.g. hydrology to estimate 35 return levels and their uncertainties. The key strength of classical EVA is that it provides a 36 theoretically founded asymptotic statistical distribution of the return levels. However, the 37 parameters of this distribution are estimated from the largest values in the data only. In practical 38 applications, the available data is typically limited, as most real-world time series are only around 39 100 years long or even less. As a result, there are very few values available for estimating the 40 parameters, leading to considerable uncertainties in the estimated return levels, especially for 41 large return periods. 42 Any method which can reduce these large uncertainties is welcomed. Marani and Ignaccolo 43 (2015), henceforth MI2015, introduces the metastatistical extreme value (MEV) approach with the 44 explicit aim to avoid relying on the asymptotic distribution of extremes in classical EVA. Instead, 45 the MEV approach regards extremes as large ordinary events, which are all distributed according 46 to a common statistical distribution (the parent distribution), which is assumed to be known. This 47 seemingly solves the problem of data shortage in classical EVA by increasing the amount of data 48 available for fitting substantially. In addition, the MEV allows for interannual stochastic variations 49 of the statistical properties of the distribution, hence the term metastatistical. 50 MEV was originally applied to daily rainfall (MI2015, Zorzetto et al. 2016) assuming a Weibull 51 parent. The MEV approach has subsequently been applied to rainfall (Marra et al. 2018; 52 Schellander et al. 2019; Zorzetto and Marani 2020; Miniussi and Marani 2020; Poschlod 2021; 53 Falkensteiner et al. 2023; Zorzetto et al. 2024; Devò et al. 2025). Besides, the MEV approach has 54 also been applied to streamflow (Miniussi et al. 2020a; Mushtaq et al. 2022) , hurricane 55 occurrence (Hosseini et al. 2020), rainfall connected to tropical cyclones (Miniussi et al. 2020b), 56 and sea level (Caruso and Marani 2022; Boumis et al. 2024). 57 The MEV methodology assumes that the entire range of rainfall follows a Weibull distribution. This 58 is largely motivated by physical arguments found in Wilson and Toumi (2005) even though their 59 conclusions refer specifically to heavy rainfall. Martinez-Villalobos and Neelin (2019) also 60 employed idealized physical arguments but concluded that the whole range of daily rainfall, 61 including the tail, should be described by a gamma distribution. Papalexiou et al. (2013) fitted four 62 different distributions to the tail of thousands of rainfall stations worldwide and found that none 63 of the four distributions systematically was better fit. 64 Thus, so far, there is no agreed exact form of the tail of the distribution of rainfall, as assumed in 65 the MEV approach. In our view this is not surprising, since rainfall is a complex phenomenon, 66 involving e.g. the distinction between stratiform and convective rainfall, land surface processes, 67 and their interaction. So rather than a simple distribution, we would expect rainfall to be 68 distributed according to a mixture of distributions, of which some could be unknown (Smith et al. 69 2011; Shin et al. 2015). Also river runoff is generally influenced by several processes that will have 70 major impacts on the most extreme observations (Merz and Blöschl 2008; Merz et al. 2022). Sea 71 surges could also be caused by different processes being important for different ranges of return 72 periods (Su et al. 2024). 73 3 It is important to have a thorough understanding of the physical processes responsible for the 74 extremes. In particular, it is crucial to consider whether the largest observations are outliers 75 resulting from very rare physical circumstances. Neglecting this could lead to underestimation of 76 exceedance rates for rare events (Merz et al. 2022) as well as their uncertainties (Merz et al. 77 2015). 78 MI2015 overlook the above considerations and therefore the need for studying how sensitive the 79 estimated return levels are to the details of the tail of the distribution. Instead, MI2015 in their 80 Monte Carlo simulation study rely exclusively on the Weibull distribution for both generating the 81 synthetic series and subsequent estimating return levels. Motivated by this we want to investigate 82 if the detailed tail behavior of the distribution used to generate the synthetic series impacts the 83 results of the MEV approach and the classical EVA. We will replicate the analysis framework of 84 MI2015 and compare the MEV and the classical EVA methods on Weibull-distributed synthetic 85 series, and in addition will we will do similar analysis on synthetic series distributed according to a 86 Weibull distribution with a slightly perturbed tail. Furthermore, we will investigate whether 87 commonly used goodness-of-fit tests can distinguish between the Weibull and the perturbed 88 Weibull synthetic series to see if such tests could help determine if the MEV approach is 89 appropriate for a given time series. 90 Recently, another critique of the MEV approach (Serinaldi et al. 2025) appeared. They showed that 91 if accounting for autocorrelation, the MEV becomes strongly biased. Our critique should be 92 regarded complementary to their. 93 2 Methods 94 Section 2.1 describes the block maximum method of classical EVA and Section 2.2 the MEV 95 methodology. Section 2.3 describes the Weibull and the perturbed Weibull mixture distributions 96 used to generate synthetic time series. Section 2.4 describes the use of goodness-of-fit tests. 97 Finally, Section 2.5 describes our protocol for performing the Monte Carlo simulations. 98 2.1 Classical EVA, the block maximum method 99 The block maximum method of classical extreme value time series analysis considers blocks of 𝑛 100 independent and identically distributed random variables {𝑥𝑖}𝑖=1,𝑛 distributed according to the 101 parent cumulative distribution function (CDF) 𝐹. The aim is to find the distribution of their 102 maximum: 103 𝑦𝑛=max{𝑥1,…,𝑥𝑛}. (1) 104 The above presumptions lead to the CDF of 𝑦𝑛 being 105 𝐻𝑛(𝑦)=[𝐹(𝑦)]𝑛. (2) 106 In principle one could estimate 𝐹 from observations and use Eq. 2 to get 𝐻𝑛. However, as Coles 107 (2001) points out, very small discrepancies in the estimate of 𝐹 can lead to substantial 108 discrepancies for 𝐹𝑛. 109 4 In classical EVA it is therefore acknowledged that the details of 𝐹, and in particular its tail, are 110 largely unknown. Instead, the Fisher–Tippett–Gnedenko theorem (Fisher and Tippett 1928; 111 Gnedenko 1943) is used to obtain a limit distribution of 𝑦𝑛 as as 𝑛→∞. This limit distribution is 112 the generalized extreme value (GEV) distribution: 113 𝐺(𝑦)=exp[− (1+𝜉𝑦−𝜇 𝜎)−1 𝜉]. (3) 114 Thus, if 𝑛 is large enough, then 𝐻𝑛(𝑦)≈𝐺(𝑦) to any desired accuracy regardless of 𝐹. For further 115 details, see Coles (2001). 116 The parameters 𝜇 and 𝜎 are the location and scale parameter, respectively, while the shape 117 parameter 𝜉 distinguishes between different types of tail behavior of 𝐺. Thus, if 𝜉=0 then 𝐺 is 118 the Gumbel or EV1 distribution with exponential tail behavior. If 𝜉>0 then 𝐺 is the Fréchet or 119 EV2 distribution with a heavy tail. Finally, if 𝜉<0 then 𝐺 is the reversed Weibull or EV3 120 distribution, which has a bounded upper tail. 121 In practical applications with time series we split the original series of length 𝑁 years into 122 consecutive blocks of one year length to get a series of annual maxima, {𝑦𝑘 }𝑘=1,𝑁. Therefore 𝑛 123 from now on denotes the number of data values per year. The three parameters of 𝐺, or in the 124 Gumbel case with 𝜉=0, only the remaining two parameters, are estimated using standard 125 techniques. We will refer to GEV and Gumbel fitting to annual maxima as the GEV and Gumbel 126 method, respectively. 127 For design and risk assessment purposes, it is common to estimate the return level, 𝑦𝑇, 128 corresponding to a return period 𝑇, which is defined in terms of the exceedance probability as 129 1−𝐺(𝑦𝑇)=1 𝑇. (4) 130 This means that 𝑦𝑛, and therefore also the original time series 𝑥𝑖, exceeds 𝑦𝑇 on average one time 131 during a time span of 𝑇 years. We can isolate 𝑦𝑇 from Eq. 4 to get 132 𝑦𝑇=𝐺−1(1−1 𝑇). (5) 133 2.2 The metastatistical extreme value (MEV) methodology 134 In the MEV methodology (MI2015), all data values, including the extremes, are described by a 135 specified parent distribution, 𝐹 - the base distribution. The corresponding extreme value CDF of 136 annual maxima, 𝐻𝑛, is obtained using Eq. (2), which includes fitting the base distribution to data. 137 The MEV also allows for stochastically varying 𝑛,and parameters 𝜗1…𝜗𝑘 of 𝐹. Therefore, the MEV 138 CDF of annual maxima is defined as the ensemble average of 𝐻𝑛, 139 𝜁(𝑦)=∑∫𝑑𝜗1…𝑑𝜗𝑘 𝑛𝑔(𝑛,𝜗1…𝜗𝑘)𝐻𝑛(𝑦;𝑛,𝜗1…𝜗𝑘) 141 =∑∫𝑑𝜗1…𝑑𝜗𝑘𝑛 𝑔(𝑛,𝜗1…𝜗𝑘)[𝐹( 𝜗1…𝜗𝑘)]𝑛, (6) 140 5 which is the general formulation of the MEV extreme value distribution. Here 𝑔(𝑛,𝜗1…𝜗𝑘) is the 142 joint density distribution of the parameters 𝑛 and 𝜗1…𝜗𝑘. If we assume ergodicity, we can 143 approximate the integration over phase space by a time average over the number of years, 𝑁, 144 yielding 145 𝜁(𝑦)≈1 𝑁∑𝐻𝑛𝑗(𝑦;𝑛𝑗,𝜗1𝑗…𝜗𝑘𝑗) 𝑁 𝑗=1 =1 𝑁∑𝐹(𝑦; 𝜗1𝑗…𝜗𝑘𝑗)𝑛𝑗 𝑁 𝑗=1 , (7) 146 where ϑ1j…ϑkj and nj vary among years. The approximate formulation in Eq. 7 eliminates the 147 need for knowing 𝑔, and Eq. 7 can provide percentiles of ζ and therefore the associated return 148 levels from Eq. 5. 149 2.3 The Weibull and the perturbed Weibull parents 150 MI2015 use a two-parameter Weibull distribution with CDF 151 𝐹𝑊(𝑥;𝑤,𝐶)=1−exp[−(𝑥 𝐶)𝑤] (8) 152 as base distribution in the MEV methodology, where 𝑤 and 𝐶 are the Weibull shape and scale 153 parameter, respectively. The parameters can either be constant or vary from year to year. Do not 154 confuse these parameters with the different parameters with the same names used in the GEV 155 and Gumbel distributions. 156 We aim to investigate whether the MEV methodology and classical GEV deliver satisfactory results 157 also when data are drown from a parent distribution with a slightly perturbed tail, while still using 158 the original base distribution for fitting. We therefore construct a perturbed parent distribution, 159 𝐹𝑊𝑝, where the original (base) Weibull distribution (Eq. 8) is mixed with another Weibull 160 distribution with a stretched tail. Mathematically, this is expressed as: 161 𝐹𝑊𝑝(𝑥)= (1−𝛼)𝐹𝑊(𝑥;𝑤,𝐶)+𝛼𝐹𝑊(𝑥;𝑤,𝛽𝐶), (9) 162 where we require the mixing ratio 𝛼≪1 and the stretching factor 𝛽>1.Both MI2015 and the 163 present study use parameters derived from the long series of daily rainfall from Padova, Italy 164 (Camuffo 1984). Thus, in the simplest setup, we use constant parameters 𝐶=7.3 mm and 𝑤= 165 0.82 as in MI2015. Furthermore, we use 𝛼=0.02 and 𝛽=3. 166 6 167 Figure 1 PDF of Weibull (red line) with 𝐶=7.3 𝑚𝑚, 𝑤=0.82 and corresponding perturbed 168 Weibull (blue line) with 𝛼=0.02 and 𝛽=3. The dashed line (right y-axis) is the ratio between the 169 survival function (see text for explanation) of Weibull and perturbed Weibull distributions. 170 Figure 1 shows the PDF of both the Weibull and the perturbed Weibull as colored lines. This 171 demonstrates that the two PDFs overall have only hardly noticeable differences. As we will see in 172 Section 3.3, the difference between 𝐹𝑊 and 𝐹𝑊𝑝 is so small that commonly used goodness-of-fit 173 tests cannot separate synthetic rainfall series generated from the two distributions. However, 174 focusing on the tails of the two distributions by considering the survival function, defined as 175 Ψ=1−𝐹, 176 a different picture emerges. Thus, the ratio Ψ𝑊𝑝 Ψ𝑊 ⁄(dotted line in Figure 1) is close to unity for 177 small values of synthetic daily rainfall, but from around 50 mm it increases exponentially (note 178 logarithmic scale on the right y-axis), illustrating that F𝑊𝑝 has significantly more mass in its tail 179 compared to F𝑊. 180 We combine Eqs. (8) and (9) with (2) to get an expression for the theoretical CDFs of annual 181 maxima, assuming 𝑛=100 wet-days each year. From this, we calculate the true return levels for 182 annual maxima, shown in Figure 2. 183 7 184 Figure 2 True annual maximum return levels as a function of return period (logarithmic axis) from 185 Weibull parent with 𝐶=7.3 𝑚𝑚, 𝑤=0.82 and 𝑛=100 and from the corresponding perturbed 186 Weibull parent with 𝛼=0.02 and 𝛽=3, using Eq. 2. 187 The true return levels shown in in Figure 2 behave qualitatively different. The return levels from 188 the Weibull parent form a straight line, which means that they are close to their asymptotic 189 distribution which is a Gumbel distribution (e.g. Embrechts et al. 1997). The true return levels of 190 the perturbed Weibull do not align along a straight line, meaning that they are farther from 191 convergence (also a Gumbel distribution). 192 The true return levels of the Weibull and perturbed Weibull parents are close to each other for 193 return periods around 10 year and lower return period but the difference increases progressively 194 with a longer return period. For 100 years return period the perturbed Weibull parent return level 195 is between 30% and 40% above that of the Weibull parent, increasing to nearly 80% at 1000 years. 196 This difference in tail behavior is comparable in magnitude to the difference in tail behavior found 197 by (Papalexiou et al. 2013) when fitting different distributions to a daily rainfall series. This 198 suggests that our construction of the Weibull and perturbed Weibull distributions reflects the 199 uncertainty in tail behavior seen in real data. 200 201 8 Thus, even though the Weibull and perturbed Weibull parents are almost identical for the bulk of 202 the data, the associated extreme value distributions of annual maxima are distinctly different. 203 2.4 Goodness-of-fit test 204 Goodness-of-fit tests have been applied by Zorzetto et al. (2016) and Marra et al. (2018) to justify 205 that observed rainfall series are distributed according to the Weibull distribution, and that the 206 MEV approach therefore is reliable. We investigated whether such tests could detect that Weibull 207 and perturbed Weibull series are distributed according to different parent distributions. If the 208 tests fail to do so, they cannot be used for identifying if the MEV approach is applicable. 209 We used the two-sample version of the tests, testing the (false) null hypothesis that two random 210 series, one distributed as a Weibull and the other as a perturbed Weibull distribution, come from 211 the same distribution. The outcome of the test is a 𝑝-value that is compared with a significance 212 level 𝑝0, and if 𝑝<𝑝0 the null hypothesis is (correctly) rejected. If many pairs of series are tested, 213 we can calculate a rejection rate. Since the two series by construction are distributed according to 214 different distributions, we would in the ideal case (neglecting type II error) have a rejection rate of 215 unity, whereas we would have a rejection rate of 𝑝0 (type I error) if they were distributed 216 according to the same distribution. 217 We applied three goodness-of-fit test: the Kolmogorov-Smirnov (Kolmogorov 1933; Smirnov 218 1948), the Cramér-von Mises (Cramér 1928; von Mises 1931) and the Anderson-Darling (Anderson 219 and Darling 1954) tests to the pairs of series in the Monte Carlo protocol, generated from Weibull 220 and perturbed Weibull parent, respectively (see Section 2.5). These tests base their test statistics 221 on the differences between the empirical CDFs of the two series, but the exact form of the test 222 statistics differs between the tests, and is reflected in the properties of the tests (Razali and Wah 223 2011). 224 The Kolmogorov-Smirnov test uses the maximum absolute difference between the empirical CDFs 225 as test statistic, which makes it sensitive to differences in the central part of the distributions. Both 226 the Cramér-von Mises and the Anderson-Darling tests use weighted differences between the 227 empirical CDFs, integrated over the entire distribution as the test statistic. The weights differ 228 between the two, so that the Cramér-von Mises test is more sensitive the overall difference of the 229 distributions, while the Anderson-Darling test is more sensitive to differences in the tails of the 230 distributions. 231 2.5 Monte Carlo protocol 232 We build upon the Monte Carlo simulation protocol from MI2015 to compare the MEV approach 233 and the classical EVA (GEV and Gumbel distribution fits to annual maxima) with synthetic series 234 generated from a Weibull parent, but extend it with synthetic rainfall series generated from a 235 perturbed Weibull parent. The protocol is: 236 1. Generate 𝑁=50 years long synthetic series randomly drawn from the Weibull parent (as 237 MI2015) and from the corresponding perturbed Weibull parent, respectively. 238 9 2. Fit Weibull distributions according to Eq. 7 (MEV approach) using probability weighted 239 moments estimation to both Weibull and perturbed Weibull parent synthetic series and 240 determine return levels corresponding to 10, 100, and 100 years return period. 241 3. Derive annual maxima series from the Weibull and perturbed Weibull parent series. 242 4. Fit GEV and Gumbel distributions using maximum likelihood estimation (classical EVA) to 243 both annual maxima series and determine return levels. 244 5. Apply two-sample goodness-of-fit tests (Kolmogorov-Smirnov, Cramér-von Mises and 245 Anderson-Darling) to the Weibull/perturbed Weibull pair of series. 246 We apply, as MI2015, the above protocol for three different setups using different combinations 247 of constant or annually varying parameters. The three setups are summarized in Table 1 and 248 described in detail below. 249 Table 1: Summary of setups. In all setups, the time series are 𝑁=50 years long. 250 Setup MEV model Parent two-parameter Weibull distribution, scale (𝐶) and shape (𝑤) # wet-days per year (𝑛) A Eq. 10 Fixed at 𝐶=7.3 mm and 𝑤=0.82 Fixed at 100 B Eq. 11 Annually varying, random values obtained from Padova rainfall series Fixed at 100 C Eq. 12 Fixed at 𝐶=7.3 mm and 𝑤=0.82 Random, uniformly distributed between 21 and 50 251 In setup A both the number of wet-days per year and the parameters in both parent distributions 252 are constant through time. So the synthetic series are generated with 𝑛=100 wet-days each year 253 and constant Weibull parameters of 𝐶=7.3 mm and 𝑤=0.82, as in MI2015. In this case of 254 constant 𝐶 and 𝑤, the MEV approach (Eq. 7) simplifies to 255 𝜁𝐴(𝑦)=𝐻𝑛(𝑦)=𝐹𝑊(𝑦;𝑤,𝐶)𝑛=[1−exp(−(𝑦 𝐶)𝑤)]𝑛, (10) 256 and we estimate the values of 𝑤 and 𝐶 for the two series by fitting Weibull distributions to both 257 the Weibull and the perturbed Weibull synthetic series. 258 MI2015 argues that annual variations of the parent Weibull parameters are important. Therefore, 259 they applied setup B, where they keep the number of wet-days constant at 𝑛=100 but allow the 260 Weibull parameters to vary between years as random values obtained from the Padova rainfall 261 series (Camuffo 1984). For setup B, the MEV approach (Eq. 7) becomes 262 𝜁𝐵(𝑦)≈1 𝑁∑𝐹𝑊(𝑦;𝑤𝑗,𝐶𝑗)𝑛=1 𝑁∑[1−exp(−(𝑦 𝐶𝑗)𝑤𝑗)]𝑛 𝑁 𝑗=1 𝑁 𝑗=1 , (11) 263 16 4 Discussion 378 4.1 Validity for real-world data 379 Monte Carlo studies using synthetic time series can give valuable insights into the statistical 380 methods by enabling controlled experiments. That said, such experiments must be carefully 381 designed to properly mimic the real-world examples we want to study. Furthermore, caution is 382 essential when transferring conclusions from these studies to real data, as synthetic series may fail 383 to capture the full complexities inherent in such data. Both MI2015 and the present study assume 384 that data points for different days are independent and identically distributed (stationary) but 385 these assumptions are usually not fulfilled for real-world data. 386 Classical EVA typically uses a block length of one year and selects one data value (the largest) from 387 each block. Thus, since the decorrelation time of the time series usually is in the order of a few 388 days, the series of annual maxima can be regarded as independent. In that case the theory 389 described in Section 2.1 is still valid if the number of data points are replaced by an effective 390 number of data points. 391 Non-stationarity is another important consideration. Most prominently, most hydrological time 392 series exhibit a strong annual cycle. This means that the probability of having an extreme value 393 varies through the season, which is a major reason for using a block length of one year. That said, 394 the presence of an annual cycle requires that the results from classical EVA (see Section 2.1) be 395 modified, as discussed in Buishand (1989). Koutsoyiannis (2004) considers the illustrative example 396 of precipitation for each month are all gamma-distributed but with different parameters, which 397 leads to an overall distribution with a heavy tail and thus convergence to GEV with non-zero shape 398 parameter. 399 Interannual variations and trends are other types of non-stationarity. Classical EVA can 400 incorporate trends and influences of large-scale atmospheric circulation patterns by introducing 401 suitable covariates into the analysis (Coles 2001). 402 There is thus some justification to apply classical EVA to time series with these relaxed 403 assumptions. On the other hand, (Serinaldi et al. 2025) found the MEV method to be biased for 404 autocorrelated data. More discussion on the potential limitations of the MEV methodology would 405 be welcomed before further use of the method in practical applications. 406 4.2 Simplified MEV (SMEV) 407 A variant of MEV, simplified MEV or SMEV has been introduced by Marra et al. (2019). SMEV also 408 assumes a known parent distribution, but re-introduces the focus on the tail by applying left409 censoring defined by a threshold. The SMEV thus have similarities to the peak-over-threshold from 410 EVA . The threshold is objectively determind by a test procedure (Marra et al. 2023). 411 It may be that the SMEV with its focus on the tail does not have the shortcomings of MEV pointed 412 to in this study. On the other hand, the left-censoring limits the number of data available for 413 estimation and thus re-introduces a problem, which the the original MEV aimed to solve. 414 17 5 Summary and outlook 415 We have confirmed the conclusion of MI2015 that the MEV method provides unbiased estimates 416 for Weibull parent synthetic series for all three return periods. The distributions of the estimates 417 are more narrow compared to the classical EVA (GEV and Gumbel methods). These results led 418 MI2015 to conclude that MEV was superior to classical EVA. 419 However, this conclusion changes when studying the behavior of the perturbed Weibull parent 420 series. Now the MEV methodology severely underestimates the return levels, with the true value 421 being larger than the 95th percentile of the distribution. Consequently, the MEV method generally 422 carries the risk of being over-confident. The GEV method, on the other hand, provides unbiased 423 estimates with the true value inside the 5-95 percentile interval of the distribution. The estimates 424 of the Gumbel method are unbiased for 10 year return period, but becomes increasingly 425 negatively biased for larger return periods. This is a consequence of incomplete convergence (see 426 Figure 2) which means that a GEV fit is superior to a Gumbel fit, which leads to underestimation of 427 the higher return levels (Koutsoyiannis and Baloutsos 2000). 428 Setups B and C were formulated to account for variation of the Weibull parameters and/or 429 number of wet-days between years. Our analysis shows that in general this flexibility does not in 430 general make the MEV method more capable of catching the extremes for the Weibull mixture 431 parent series, so conclusions obtained for setup A hold also for setups B and C. This is despite the 432 assumption of stationarity, which guarantees the GEV method to be valid, is violated in setup B. 433 Another important result from the simulation study is that three commonly used goodness-of-fit 434 tests are unable to distinguish effectively between time series drawn from the Weibull and 435 Weibull mixture parent distributions. Therefore, this does not provide a way to detect series 436 where the MEV method can be used and where it should be avoided. 437 We repeated our Monte Carlo eksperiments (setup A) with different combinations of base and 438 perturbation distributions. These confirmed the conclusion that the MEV methodology is over439 confident. We therefore find et reasonable to assume that the findings of the present study 440 generalize to a large class of distributions. This implies that the risk producing biased and over441 confident return level estimates is a general caveat of the MEV methodology. 442 Acknowledgements 443 The authors would like to acknowledge the support of the Danish Government through the 444 National Center for Climate Research (NCKF) and the Danish Climate Atlas. 445 Author contributions: CRediT 446 Torben Schmith: Conceptualization; analysis; writing manuscript. Karsten Arnbjerg-Nielsen, Bo 447 Christiansen: Conceptualization; contributing to manuscript. 448 Funding sources 449 18 Funding was provided by the Danish State through the National Center for Climate Research 450 (NCKF) and the Danish Climate Atlas. 451 452 References: 453 Anderson TW, Darling DA (1954) A Test of Goodness of Fit. J Am Stat Assoc 49:765–769. 454 https://doi.org/10.1080/01621459.1954.10501232 455 Boumis G, Moftakhari HR, Moradkhani H (2024) A metastatistical frequency analysis of extreme 456 storm surge hazard along the US coastline. Coast Eng J 66:380–394. 457 https://doi.org/10.1080/21664250.2024.2338323 458 Buishand TA (1989) Statistics of extremes in climatology. Stat Neerlandica 43:1–30. 459 https://doi.org/10.1111/j.1467-9574.1989.tb01244.x 460 Camuffo D (1984) Analysis of the series of precipitation at Padova, Italy. Clim Change 6:57–77. 461 https://doi.org/10.1007/BF00141668 462 Caruso MF, Marani M (2022) Extreme-coastal-water-level estimation and projection: a comparison 463 of statistical methods. Nat Hazards Earth Syst Sci 22:1109–1128. 464 https://doi.org/10.5194/nhess-22-1109-2022 465 Coles S (2001) An introduction to statistical modeling of extreme values. Springer, London, UK 466 Cramér H (1928) On the Kolmogorov-Smirnov test for goodness of fit. Scand Actuar J 1928:1–9 467 Devò P, Caruso MF, Borga M, Marani M (2025) Estimates of Rare Rainfall Extremes in Ungauged 468 Areas. Geophys Res Lett 52:e2024GL113576. https://doi.org/10.1029/2024GL113576 469 Embrechts P, Mikosch T, Klüppelberg C (1997) Modelling extremal events: for insurance and 470 finance. Springer-Verlag, Berlin, Heidelberg 471 Falkensteiner M-A, Schellander H, Ehrensperger G, Hell T (2023) Accounting for seasonality in the 472 metastatistical extreme value distribution. Weather Clim Extrem 42:100601. 473 https://doi.org/10.1016/j.wace.2023.100601 474 Fisher RA, Tippett LHC (1928) Limiting forms of the frequency distribution of the largest or smallest 475 member of a sample. Math Proc Camb Philos Soc 24:180–190. 476 https://doi.org/10.1017/S0305004100015681 477 Gnedenko BV (1943) Sur La Distribution Limite Du Terme Maximum D’Une Serie Aleatoire. Ann 478 Math 44:423 479 Hosseini SR, Scaioni M, Marani M (2020) Extreme Atlantic Hurricane Probability of Occurrence 480 Through the Metastatistical Extreme Value Distribution. Geophys Res Lett 481 47:2019GL086138. https://doi.org/10.1029/2019GL086138 482 19 Kolmogorov A (1933) Sulla determinazione empirica di una leggi di distribuzione, Giorn. 1st Ital. 483 Attuari 4:91 484 Koutsoyiannis D (2004) Statistics of extremes and estimation of extreme rainfall: I. Theoretical 485 investigation / Statistiques de valeurs extrêmes et estimation de précipitations extrêmes: I. 486 Recherche théorique. Hydrol Sci J 49:3. https://doi.org/10.1623/hysj.49.4.575.54430 487 Koutsoyiannis D, Baloutsos G (2000) Analysis of a long record of annual maximum rainfall in 488 Athens, Greece, and design rainfall inferences. Nat Hazards 22:29–48 489 Marani M, Ignaccolo M (2015) A metastatistical approach to rainfall extremes. Adv Water Resour 490 79:121–126. https://doi.org/10.1016/j.advwatres.2015.03.001 491 Marra F, Amponsah W, Papalexiou SM (2023) Non-asymptotic Weibull tails explain the statistics of 492 extreme daily precipitation. Adv Water Resour 173:104388. 493 https://doi.org/10.1016/j.advwatres.2023.104388 494 Marra F, Nikolopoulos EI, Anagnostou EN, Morin E (2018) Metastatistical Extreme Value analysis of 495 hourly rainfall from short records: Estimation of high quantiles and impact of measurement 496 errors. Adv Water Resour 117:27–39. https://doi.org/10.1016/j.advwatres.2018.05.001 497 Marra F, Zoccatelli D, Armon M, Morin E (2019) A simplified MEV formulation to model extremes 498 emerging from multiple nonstationary underlying processes. Adv Water Resour 127:280– 499 290. https://doi.org/10.1016/j.advwatres.2019.04.002 500 Martinez-Villalobos C, Neelin JD (2019) Why Do Precipitation Intensities Tend to Follow Gamma 501 Distributions? J Atmospheric Sci 76:3611–3631. https://doi.org/10.1175/JAS-D-18-0343.1 502 Merz B, Basso S, Fischer S, et al (2022) Understanding Heavy Tails of Flood Peak Distributions. 503 Water Resour Res 58:e2021WR030506. https://doi.org/10.1029/2021WR030506 504 Merz B, Vorogushyn S, Lall U, et al (2015) Charting unknown waters—On the role of surprise in 505 flood risk assessment and management. Water Resour Res 51:6399–6416. 506 https://doi.org/10.1002/2015WR017464 507 Merz R, Blöschl G (2008) Flood frequency hydrology: 2. Combining data evidence. Water Resour 508 Res 44:. https://doi.org/10.1029/2007WR006745 509 Miniussi A, Marani M (2020) Estimation of Daily Rainfall Extremes Through the Metastatistical 510 Extreme Value Distribution: Uncertainty Minimization and Implications for Trend 511 Detection. Water Resour Res 56:e2019WR026535. 512 https://doi.org/10.1029/2019WR026535 513 Miniussi A, Marani M, Villarini G (2020a) Metastatistical Extreme Value Distribution applied to 514 floods across the continental United States. Adv Water Resour 136:103498. 515 https://doi.org/10.1016/j.advwatres.2019.103498 516 20 Miniussi A, Villarini G, Marani M (2020b) Analyses Through the Metastatistical Extreme Value 517 Distribution Identify Contributions of Tropical Cyclones to Rainfall Extremes in the Eastern 518 United States. Geophys Res Lett 47:e2020GL087238. 519 https://doi.org/10.1029/2020GL087238 520 Mushtaq S, Miniussi A, Merz R, Basso S (2022) Reliable estimation of high floods: A method to 521 select the most suitable ordinary distribution in the Metastatistical extreme value 522 framework. Adv Water Resour 161:104127. 523 https://doi.org/10.1016/j.advwatres.2022.104127 524 Papalexiou SM, Koutsoyiannis D, Makropoulos C (2013) How extreme is extreme? An assessment 525 of daily rainfall distribution tails. Hydrol Earth Syst Sci 17:851–862. 526 https://doi.org/10.5194/hess-17-851-2013 527 Poschlod B (2021) Using high-resolution regional climate models to estimate return levels of daily 528 extreme precipitation over Bavaria. Nat Hazards Earth Syst Sci 21:3573–3598. 529 https://doi.org/10.5194/nhess-21-3573-2021 530 Razali NM, Wah YB (2011) Power comparisons of shapiro-wilk, kolmogorov-smirnov, lilliefors and 531 anderson-darling tests. J Stat Model Anal 2:21–33 532 Ross R (1994) Formulas to describe the bias and standard deviation of the ML-estimated Weibull 533 shape parameter. IEEE Trans Dielectr Electr Insul 1:247–253 534 Schellander H, Lieb A, Hell T (2019) Error Structure of Metastatistical and Generalized Extreme 535 Value Distributions for Modeling Extreme Rainfall in Austria. Earth Space Sci 6:1616–1632. 536 https://doi.org/10.1029/2019EA000557 537 Serinaldi F, Lombardo F, Kilsby CG (2025) Non-asymptotic distributions of water extremes: much 538 ado about what? Hydrol Earth Syst Sci 29:1159–1181. https://doi.org/10.5194/hess-29539 1159-2025 540 Shin J-Y, Lee T, Ouarda TBMJ (2015) Heterogeneous Mixture Distributions for Modeling 541 Multisource Extreme Rainfalls. J Hydrometeorol 16:2639–2657. 542 https://doi.org/10.1175/JHM-D-14-0130.1 543 Smirnov N (1948) Table for estimating the goodness of fit of empirical distributions. Ann Math Stat 544 19:279–281 545 Smith JA, Villarini G, Baeck ML (2011) Mixture Distributions and the Hydroclimatology of Extreme 546 Rainfall and Flooding in the Eastern United States. J Hydrometeorol 12:294–309. 547 https://doi.org/10.1175/2010JHM1242.1 548 Su J, Poulsen B, Nielsen JW, et al (2024) Integrating historical storm surge events into flood risk 549 security in the Copenhagen region. Weather Clim Extrem 45:100713. 550 https://doi.org/10.1016/j.wace.2024.100713 551 21 von Mises R (1931) La distribution de la moyenne d’un échantillon d’observations indépendantes 552 et identiques. Ann Inst H Poincaré 1931:271–300 553 Wilson PS, Toumi R (2005) A fundamental probability distribution for heavy rainfall. Geophys Res 554 Lett 32:. https://doi.org/10.1029/2005GL022465 555 Zorzetto E, Botter G, Marani M (2016) On the emergence of rainfall extremes from ordinary 556 events. Geophys Res Lett 43:8076–8082. https://doi.org/10.1002/2016GL069445 557 Zorzetto E, Canale A, Marani M (2024) A Bayesian non-asymptotic extreme value model for daily 558 rainfall data. J Hydrol 628:130378. https://doi.org/10.1016/j.jhydrol.2023.130378 559 Zorzetto E, Marani M (2020) Extreme value metastatistical analysis of remotely sensed rainfall in 560 ungauged areas: Spatial downscaling and error modelling. Adv Water Resour 135:103483. 561 https://doi.org/10.1016/j.advwatres.2019.103483 562 563