Convergence of extreme value statistics in a two-layer quasi-geostrophic atmospheric model

Con v ergence of extreme v alue statistics in a t w o-la y er quasi-geostrophic atmospheric mo del V era Melinda G´ alﬁ (1,2), T am´ as B´ odai (3), V alerio Lucarini (1 , 3 , 4) (1) Meteorological Institute, CEN, Universit y of Hamburg, Ham burg, German y , (2) IMPRS-ESM, Max Planck Institute for Meteorology , Ham burg, Germany , (3) Departmen t of Mathematics and Statistics, Universit y of Reading, Reading, UK, (4) Cen ter for Environmen tal Policy , Imp erial College London, London, UK July 11, 2018 Abstract W e search for the signature of univ ersal prop erties of extreme even ts, theoretically predicted for Axiom A ﬂo ws, in a chaotic and high dimensional dynamical system b y studying the con vergence of GEV (Generalized Extreme V alue) and GP (Generalized P areto) shap e parameter estimates to a theoretical v alue, expressed in terms of partial dimensions of the attractor, whic h are global prop erties. W e consider a tw o lay er quasi- geostrophic (QG) atmospheric mo del using t wo forcing lev els, and analyse extremes of diﬀeren t types of physical observ ables (lo cal, zonally-a veraged energy , and the av erage v alue of energy ov er the mid-latitudes). Regarding the predicted universalit y , w e ﬁnd closer agreement in the shap e parameter estimates only in the case of strong forcing, pro ducing a highly chaotic b eha viour, for some observ ables (the lo cal energy at ev ery latitude). Due to the limited (though very large) data size and the presence of serial correlations, it is diﬃcult to obtain robust statistics of extremes in case of the other observ ables. In the case of weak forcing, inducing a less pronounced c haotic ﬂow with regime b eha viour, w e ﬁnd worse agreemen t with the theory developed for Axiom A ﬂows, whic h is unsurprising considering the properties of the system. 1 In tro duction and motiv ation The in vestigation of extreme even ts is extremely relev ant for a range of disciplines in mathe- matical, natural, so cial sciences and engineering. Understanding the large ﬂuctuations of the system of interest is of great imp ortance from a theoretical p oint of view, but also when it comes to assess the risk asso ciated to low probabilit y and high impact even ts. In man y cases, in order to gauge preparedness and resilience prop erly , one would lik e to b e able to quantify the return times for even ts of diﬀerent in tensity , and tak e suitable measures for preven ting the expected impacts. Prominen t examples are w eather and climate extremes, whic h can ha ve a h uge impact on human so ciet y and natural ecosystems. The present uncertaint y in the future pro jections of extremes mak es their study even more urgent and crucial [1]. In practical terms, the main goal behind the study of the extreme ev ents is to understand the 1 prop erties of the highest quan tiles of the v ariable of interest. A fundamen tal drawbac k comes from the fact that extreme ev ents are rare, so that it is diﬃcult to collect satisfactory statistics from the analysis of a time series of ﬁnite length. Additionally , in the absence of a strong mathematical framework, it is virtually impossible to mak e quan titativ e statemen ts about the probabilit y of occurrence of ev ents larger than observ ed. Therefore, statistical inference based on empirical mo dels tends to suﬀer from the lack of predictive p o w er. A robust theoretical framew ork for analysing extreme even ts is provided b y extreme v alue theory (EVT). After the early contributions b y Fisher and Tipp ett [2], EVT w as introduced by Gnedenk o [3], who disco vered that under rather general conditions the extreme ev ents asso ciated to sto c hastic v ariables can b e describ ed by a family of parametric distributions. As a result, the problem of studying the extremes of a stochastic v ariable can be reduced to estimating the parameters of a kno wn probability distribution. The most imp ortant parameter of such a distribution is called the shap e parameter, which determines the qualitativ e prop erties of the distribution, and, in particular, whether it has a ﬁnite or inﬁnite upper p oin t, or, more concretely , whether extremes are b ounded b y an absolute ﬁnite maximum or not. The t w o most important prop erties of extremes, i.e., that they occur rarely and their magni- tude is un usually high or lo w, constitute the basis of tw o popular methods of EVT, the blo c k maxima (BM) and the p eak-ov er-threshold (POT) approaches. The BM approach aims at ﬁnding the limiting distributi on of maxima M n of independent iden tically distributed random (i.i.d.r.) v ariables separated into blo c ks of size n , as n → ∞ . Under rather general condi- tions and for many diﬀeren t parent distributions, the limiting distributions of blo c k maxima b elong to the Generalized Extreme V alue (GEV) distribution family; this is the p oin t of view originally prop osed b y Gnedenko [3]. The POT approach aims at understanding the statistical prop erties of the exceedances of a sto c hastic v ariable X ab ov e a given threshold u . Under the same conditions as for the BM approach, if u is v ery high, the distribution of the threshold exceedances X − u belongs to the Generalized Pareto (GP) distribution family [4, 5]. The existence of w ell-deﬁned functional forms for the distributions describing extreme ev ents pro vides predictiv e p ow er: one can in principle compute the return time of some y et unobserv ed even ts. It is remark able that these t wo p oin ts of view on extremes (whic h lead to diﬀeren t selection pro cedures and diﬀerent choice of the ev ents classiﬁed as extremes) are, in fact, equiv alen t. Problems in applying EVT to actual time series result from the fact that, typically , the ob- serv ed data feature a certain degree of serial correlations [6]. Note that the setting b ehind the construction of EVT can b e extended by relaxing the h yp othesis of indep endence of the ran- dom v ariables, as long as correlations are deca ying suﬃcien tly fast. This is of clear relev ance when trying to use EVT for studying observ ables of deterministic dynamical systems. In this case, in fact, the underlying dynamics determines the existence of correlations b et ween the v alues of observ ables at diﬀerent times, and one can easily guess that, when the dynamical system is c haotic, there is go o d hop e of deriving EVT for its observ ables [7]. Obtaining the true limiting EVT can b e extremely hard, ev en in simple dynamical systems [8]. When analysing ﬁnite time series, the con vergence of the estimated GEV or GP shap e parameters to the asymptotic v alues can b e very slo w. The sp eed of conv ergence dep ends on the t yp e of parent distribution [9], and can b e additionally slow ed down b y correlations [10, 11]. Due to the fact that the data size is alw a ys limited, there is typically a diﬀerence b etw een the asymptotic GEV or GP parameters and the estimated ones; ﬁnite-size estimates are generally biased. F or example, the GEV shap e parameter of a simple Gaussian pro cess is 0, but, for 2 an y ﬁnite time series, we w ould estimate t ypically a negative shap e parameter [2]. When performing statistical inference using the BM or POT metho d (ﬁtting the GEV or GP mo del, resp ectiv ely , to data), it is crucial to hav e an appropriate proto col of selection of “go o d” candidates for extremes [11]. On the one hand, if the chosen blo c ks (for the BM metho d) are to o short or the threshold (for the POT metho d) is to o low, the appro ximation of the limit model is lik ely to b e inappropriate, leading to false parameter estimates. Hence, the v eriﬁcation of the agreement betw een the statistical mo del and the a v ailable data is essen tial, which is often done based on go o dness-of-ﬁt tests, lik e the Kolmogorov-Smirno v [12], Anderson-Darling [13] or Pearson’s chi-squared tests [14]. On the other hand, if the blo c ks are to o large or the threshold is to o high, the n umber of extremes ma y b e insuﬃcient for a reliable estimation of the parameters, and uncertain t y b ecomes very high. As discussed later in the pap er, Coles [11] sho ws how to derive an optimal c hoice for the v alue of the block size or the threshold, in such a wa y as to verify that we are close to the asymptotic lev el as required by EVT but w e use the a v ailable data as eﬃcien tly as p ossible. Classical EVT has been extended and adapted to analyse extremes of observ ables of c haotic dynamical systems, where the sensitiv e dep endence on initial conditions is fundamen tally resp onsible for generating a de-facto sto chastic pro cess. The reader is referred to Lucarini et al. [7] for a detailed o verview of the ﬁeld of EVT for dynamical systems. It is possible to establish an EVT for observ ables of chaotic dynamical systems when one considers Axiom A dynamical systems. These are rather sp ecial c haotic dynamical systems that are uniformly h yp erb olic on their attractor (so that stable and unstable directions are w ell separated), whic h supp orts a so-called Sinai-Ruelle-Bow en (SRB) measure. Such an inv ariant measure has ph ysical relev ance because it is stable against w eak stochastic perturbations [15, 16]. One of the great merits of Axiom A dynamical systems is that they allo w for deriving rigorous and robust statistical mec hanical prop erties for purely deterministic bac kground dynamics. Despite ha ving deterministic dynamics, when lo oking at their observ ables, they behav e just lik e generators of sto c hastic processes. While Axiom A systems are rather sp ecial and indeed not generic, they ha ve great relev ance for applications if one takes in to accoun t the chaotic h yp othesis, which indicates that high-dimensional c haotic syste ms b eha v e at all practical purp oses as if they were Axiom A [17, 18]. Sev eral studies dealing with EVT for dynamical systems rev eal a link b et w een the statistical prop erties of the extremes and geometric (and p ossibly in turn global dynamical) character- istics of the system pro ducing these extremes [8, 19, 20, 21, 22]. The main ﬁndings are that when suitable observ ables are chosen for the dynamical system of interest, it is p ossible to re- late the GEV or GP parameters describing the extremes to basic prop erties of the dynamics, and esp ecially to the geometry of the attractor. In particular, depending on the c hoice of the observ able, one can asso ciate the most imp ortan t parameter of the GEV or GP distribution to the information dimension of the attractor or to the partial information dimension along the stable and unstable directions of the ﬂo w [22]. These partial dimensions are well-deﬁned ev erywhere on the chaotic attractor, possibly with a v ariation with lo cation, also for non- uniformly hyperb olic systems [23], b eside Axiom A systems. Ho wev er, Axiom A systems p ossess an ergo dic SRB measure whic h lends itself to a universalit y of the shap e parameter for all suﬃcien tly smooth observ ables; the lo cal or point-wise (partial) dimensions taking the same v alue almost everywhere [24]. In this case the uniform shap e parameter can b e related to the (partial) Kaplan-Y orke dimension(s) whic h is (are) deﬁned by the global dynamical c haracteristic num bers, the Lyapuno v exp onents. Clearly , this is an asymptotic result, and 3 one must exp ect that diﬀerences emerge on pre-asymptotic lev el when diﬀerent observ ables are studied. In other w ords, this theory do es not make predictions regarding the con vergence of shap e parameter estimates, the analysis of whic h is the main ob jectiv e of this pap er. Via the connection with fractal dimensions, it can b e said that the analysis of extremes acts as a microscop e able to assess the ﬁne scale properties of the in v arian t measures. Some preliminary numerical tests sho w that the con vergence to the asymptotic shap e pa- rameter is slow in low-dimensional cases, in systems well known not to be Axiom A [22]. B´ odai [25] examined the conv ergence to the GEV distribution in the case of extremes of site v ariables in the Lorenz 96 mo del [26], inv estigating separately a range of cases extending from weak to strong chaos. He found that when considering conﬁgurations supp orting weak c haos with a lo w-dimensional attractor, the theoretical results obtained in the con text of the Axiom A hypothesis are hard to v erify . F or low er dimensions, up to a dimension of ab out 5, shap e parameter estimates ﬂuctuate greatly rather than conv erge, while blo c k maxima data can b e shown not to conform to a GEV mo del; and for somewhat larger dimensions, up to 9 in the study , estimates could diverge from the predicted v alue while data already conform to a GEV mo del. Go od agreement with the theory was found only in the highly turbulent case p ossessing a higher-dimensional attractor, ab out 30, supp orting the basic idea b ehind the chaotic h yp othesis. Also in this case, nonetheless, very slow con vergence w as found. In previous analysis p erformed on higher dimensional, in termediate complexity mo dels with O(10 2 − 10 3 ) degrees of freedom, very slow (if an y) conv ergence to EVT distributions could b e found in the case of extremes of lo cal temp erature observ ables [27]. In another analysis of a similar mo del [28], the agreeme n t of the distribution of global energy extremes with a mem b er of the GEV family was indeed go od, yet large uncertaint y remained on the v alue of the shap e parameter, and no stringen t test w as made to mak e sure that the estimate w as stable against changes in the blo c k size considered in the BM analysis. Clearly , the sp eciﬁc c hoice of the observ able and the degree of chaoticit y of the underlying dynamics is of primary relev ance regarding the con vergence to the limiting GEV or GP distribution. In this work, we use a quasi-geostrophic (QG) atmospheric mo del of in termediate complexity featuring 1056 degrees of freedom, to analyse extremes of diﬀeren t types of observ ables: lo cal energy (deﬁned at eac h grid point), zonally-a v eraged energy , and the av erage v alue of energy o ver the mid-latitudes. Our main ob jective is to compare the estimated GEV and GP shap e parameters with a shap e parameter deriv ed, based on the theory referred to ab o v e, from the prop erties of the attractor along the stable, unstable, and neutral directions. W e refer to this as the “theoretical shap e parameter”. Thus we explore n umerically the link b et ween the purely statistical prop erties of extreme ev ents based on EVT and the dynamical properties of the system pro ducing these extreme s. W e p erform sim ulations applying tw o diﬀeren t lev els of forcing: a strong forcing, pro ducing a highly c haotic behaviour of the system, and a w eak forcing, pro ducing a less pronounced c haotic b ehaviour. The dimensionalit y of the attractor is muc h larger in the former than in the latter case. This work go es b ey ond the previously men tioned studies based on more simple dynamical systems, in a sense that with our mo del w e can study the conv ergence for observ ables b eing diﬀerent ph ysical quantities, or, repre- sen ting diﬀeren t spatial scales/characteristics of the same physical quantit y . Additionally , compared to previous studies also p erformed on in termediate complexity models, w e consider longer time series and a v ariet y of observ ables. Our mo del is simple compared to a GCM (General Circulation Mo del), but contains t wo of the main processes relev an t for mid-latitude atmospheric dynamics, i.e., baro clinic and barotropic instabilities. Hence, w e contribute to 4 bridging the gap b et ween the analysis of extremes in simple and very high dimensional dy- namical systems, as in the case of the GCMs used for atmospheric and climate sim ulations, b y using a mo del that sim ulates to a certain degree Earth-like atmospheric pro cesses and allo ws also for computing with feasible computational costs some dynamical system prop erties, lik e Ly apunov Exp onen ts or Kaplan-Y orke Dimensions. The prop erties of the mo del hav e b een extensiv ely studied by Sc hubert and Lucarini [29, 30]. Based on n umerical results (Sebastian Sc hubert, p ersonal comm unication), the mo del is ex- p ected to b e non-hyperb olic, but one can assume that the chaotic hypothesis applies to it (in case of suﬃciently high forcing lev els inducing a chaotic b ehaviour of the system), and so in analysing the conv ergence of shap e parameter estimates one can take the predicted theo- retical shap e parameter as reference. W e also assume that the symmetry of the model with resp ect to longitude, that introduces a central direction b esides the directions of expansion and contraction in phase space, do es not alter the ergo dicit y of the system at a practical lev el, and hence the true shap e parameter remains uniform. Although we use an idealised mo del, our results are transferable to time series obtained from more realistic mo del simulations or from measuremen ts. By understanding the diﬀerences among the analysed observ ables, we gain insight into the statistical prop erties of extremes of geoph ysical observ ables with diﬀeren t spatial scales. By using t w o forcings, we are able to study the con vergence to theoretical shap e parameters related to diﬀerent chaotic systems: one exhibiting fast decaying correlations and another one c haracterised b y slo wer deca ying correlations. These asp ects are relev an t in the case of geoph ysical applications where one deals also with time series on sev eral spatial scales and with diﬀerent degrees of correlations. The s tructure of this article is as follows. Sec. 2 giv es a theoretical o verview, describing the blo c k maxima and the p eak ov er threshold approaches. In Sec. 3 we present our mo del, the p erformed sim ulations and the applied metho ds. In Sec. 4 we discuss our results regarding the statistics of extremes for strong forcing and for weak forcing. W e summarise and discuss our results in Sec. 5. 2 Some Elemen ts of Extreme V alue Theory Let us consider M n = max { X 1 , ..., X n } , where X 1 , ..., X n is a sequence of i.i.d.r. v ariables with common distribution function F ( x ). The extremal t yp es theorem [2, 3] states that if there exist sequences of constants { a n > 0 } and { b n } , so that the distribution of normalised M n , i.e., Pr { ( M n − b n ) /a n ≤ z } , con verges for n → ∞ to a non-degenerate distribution function G ( z ), then G ( z ) is one of three possible t yp es of so-called extreme v alue distributions, ha ving the cumulativ e distribution function G ( z ) = ( exp n −  1 + ξ  z − µ σ  − 1 /ξ o for ξ 6 = 0 , exp  − exp  −  z − µ σ  for ξ = 0 , (1) where −∞ < µ < ∞ , σ > 0, 1 + ξ ( z − µ ) /σ > 0 for ξ 6 = 0 and −∞ < z < ∞ for ξ = 0 [11]. G ( z ) represen ts the GEV family of distributions with three parameters: the location param- eter µ , scale parameter σ , and shape parameter ξ . The shape parameter ξ describ es the tail 5 b eha viour, and determines to whic h one of the three t yp es of extreme v alue distributions G ( z ) b elongs. If ξ = 0, the tail decays exponentially , G ( z ) is a t yp e I extreme v alue distribution or Gum b el distribution. If ξ > 0, the tail deca ys p olynomially , and G ( z ) belongs to the t yp e I I or F r ´ echet distribution. If ξ < 0, the domain of the distribution has an upp er limit, and is referred to as a t yp e I I I or W eibull distribution. Under the same conditions, for which the distribution of M n con verges to the GEV dis- tribution, the exceedances y = X − u of a threshold u reaching the upp er right p oint of the distributions of X , giv en that X > u , are asymptotically distributed according to the Generalized Pareto (GP) distribution family [11] H ( y ) =    1 −  1 + ˜ ξ y ˜ σ  − 1 / ˜ ξ for ˜ ξ 6 = 0 , 1 − exp  − y ˜ σ  for ˜ ξ = 0 , (2) where 1+ ˜ ξ y / ˜ σ > 0 for ˜ ξ 6 = 0, y > 0, and ˜ σ > 0. H ( y ) has t w o parameters: the scale parameter ˜ σ and the shap e parameter ˜ ξ . The shap e parameter ˜ ξ describ es again the tail behaviour, and determines to which one of the three types of GP distributions H ( y ) b elongs. If ˜ ξ = 0, the tail of the distribution deca ys exp onen tially; if ˜ ξ > 0, the tail decays p olynomially; and if ˜ ξ < 0 the distribution is bounded [4, 5, 31]. If con vergence to the GEV and GP distributions is realised, ˜ ξ = ξ and ˜ σ = σ + ξ ( u − µ ). As a result, once w e estimate the parameters for the GEV, we can derive the corresp onding GP parameters, and vice versa [11]. F rom the v alues of the GEV or GP parameters it is possible to infer the exp ected return lev els or extreme quantiles. Return lev els z p are obtained from the GEV distribution b y inv erting equation (1): z p = ( µ − σ ξ h 1 − y − ξ p i for ξ 6 = 0 , µ − σ log y p for ξ = 0 , (3) where y p = − log (1 − G ( z p )), and 1 /y p represen ts the return p erio d. In the case of the GP distribution, the m -observ ation return level z m (i.e., the level that is exceeded on a verage ev ery m observ ations) can be deriv ed from equation (2): z m = ( u − ˜ σ ˜ ξ h 1 − ( 1 mζ u ) − ˜ ξ i for ˜ ξ 6 = 0 , u − ˜ σ log( 1 mζ u ) for ˜ ξ = 0 , (4) where ζ u represen ts the probability of an individual observ ation exceeding the threshold u [11]. By plotting the GEV (GP) return level z p ( z m ) against the return p erio d 1 /y p ( mζ u ) on a logarithmic scale, the plot is linear if ξ = 0 ( ˜ ξ = 0), is conv ex if ξ > 0 ( ˜ ξ > 0), and is conca ve if ξ < 0 ( ˜ ξ < 0). In the case of a correlated stationary sto c hastic pro cess the same GEV limit la ws apply as for i.i.d.r. v ariables if certain conditions, regarding the deca y of serial correlation, are fulﬁlled [32, 33, 7]. By stationary , we refer to a sequence of correlated v ariables whose joint probabilit y distribution is time-inv ariant. How ev er, an imp ortan t restriction is that, as an eﬀect of serial correlation, an eﬀective blo ck size can b e deﬁned whic h is smaller than the n umber of observ ations in a blo c k. This can enhance the bias in the parameter estimation, app earing as a slo wer or delay ed con v ergence of the block maxima distribution to the limiting GEV distribution [11, 10]. Another p ossible eﬀect of serial correlation is the app earance of 6 extremes at consecutiv e time steps (clusters). If an extreme v alue law do es exist in this case, then G ∗ ( z ) = G ( z ) θ , where θ is called the extremal index and 0 < θ < 1 ( G ∗ ( z ) denotes the limiting distribution of BM from the correlated sequence and G ( z ) the one from an uncorrelated sequence, having the same marginal distribution). Clusters of extremes represen t a problem especially when applying the POT approac h. A widely-adopted method to get rid of the correlated extremes is declustering, whic h basically consists of iden tifying the maxim um excess within eac h cluster, and ﬁtting the GP distribution to the cluster maxima [33, 34, 35]. As men tioned b efore, several studies on EVT for observ ables of dynamical systems relate the GEV and GP shape parameters to certain prop erties of the attractor. In the case of the so-called “distance” observ ables, one can relate the GEV and GP parameters to basic geometrical prop erties of the attractor [8, 19, 20]. The distance observ ables g (dist( x, x 0 )) are functions of the Euclidean distance b et ween one p oin t on the attractor x 0 and the orbit x . The function g ( y ) is c hosen in a w a y to ha ve a global maxim um for y = 0, so that large v alues of g corresp ond to recurrences of the orbit near x 0 . Dep ending on the c hoice of the function g ( y ), the extremes of the distance observ ables can hav e p ositive, negative, or v anishing v alue for the shap e parameter. In particular, when g ( y ) is c hosen to b e a p ow er law, the shap e parameter is non-zero, and it is prop ortional to the inv erse of the Kaplan-Y orke dimension of the attractor; see below for further details [19, 20, 22]. While recurrence prop erties are indeed imp ortan t for characterising a system, distance ob- serv ables are not well suited for studying some basic physical prop erties, such as, in the the case of ﬂuids, energy or enstroph y . Hence, Holland et al. [21] studied the extremes of smo oth functions A = A ( x ) which tak e their maxim um on the attractor in a point where the corresp onding lev el surface of A ( x ) is tangen tial to the unstable manifold of the attrac- tor, referring to them as “ph ysical” observ ables. They found a relationship betw een the GEV shap e parameter and some geometrical prop erties of the attractor dealing with the properties of the unstable and stable directions in the tangen t space. The results of Holland et al. [21] w ere re-examined b y Lucarini et al. [22], using the POT approach for physical observ ables of Axiom A systems. They considered the time-contin uous time series of physical observ ables, and found that for all non-pathological physical observ ables A the shap e parameter can b e written as: ξ δ = − 1 δ , (5) with δ deﬁned as δ = d s + ( d u + d n ) / 2 , (6) where d s , d u , and d n are the partial dimensions of the attractor restricted to the stable, unstable, and neutral (i.e., cen tral) directions. As mentioned in Sec. 1, these local or p oin t- wise (partial) dimensions tak e the same v alue almost everywhere on the attractor if one considers smooth observ ables of Axiom A systems. d u is equal to the n umber of positive Ly apunov exponents [24], d n is equal to the n umber of zero Ly apunov exponents, and d s = d K Y − d u − d n , where d K Y = n + P n k =1 λ k | λ n +1 | is the Kaplan-Y orke dimension with λ k denoting the Ly apunov exponents of the system, in a descending order, and n is suc h that P n k =1 λ k is p ositiv e and P n +1 k =1 λ k is negative. W e remark that a more general point of view, taking into consideration p ossible geometrical degeneracies, suggests that − 1 /ξ δ < d K Y < − 2 /ξ δ , and, additionally , d K Y / 2 = ( d s + d u + d n ) / 2 ≤ δ ≤ d s + ( d u + d n ) / 2 [7]. 7 According to Eq. (5) the shape parameter is alwa ys negativ e (due to the compactness of the attractor), and it is close to zero in the case of systems having large v alues of the Kaplan- Y ork e dimension. F urthermore, it sho ws a universal prop ert y of extremes, whic h do es not dep end on the chosen observ able but only on the geometry of the attractor. In what follo ws w e will fo cus on comparing Eq. (5) with statistically inferred GEV and GP shap e parameters in the case of energy extremes of the mo del of in terest in this study , which is describ ed next. 3 Mo del description and metho ds W e consider a sp ectral quasi-geostrophic (QG) 2-lay er atmospheric mo del similar to the one in tro duced b y Phillips [36]. Sp eciﬁcally , our mo del, including the simulation co de, is the same as in Sch ub ert and Lucarini [29], and is a modiﬁed version of the one presen ted b y F risius [37]. The mo del represents synoptic scale mid-latitude atmospheric dynamics based on the quasi-geostrophic approximation, whic h assumes hydrostatic balance and allo ws only small departures from the geostrophic balance. The mo del features baro clinic con version and barotropic stabilisation pro cesses, and sim ulates a turbulent jet-like zonal ﬂo w when suitable v alues are chosen for the parameters of the system. The reader is referred to Holton [38] for a detailed ph ysical and mathematical description of quasi-geostrophic approximation for mid-latitude atmospheric dynamics. 3.1 Mo del description The mo del domain is a rectangular channel with latitudinal and longitudinal co ordinates ( x, y ) ∈ [0 L x ] × [0 L y ]. y = 0 represents the Equator, and y = L y corresp onds to the north P ole. W e assume p eriodicity along the x -direction, so that L x corresp onds to the length of the parallel at 45 ◦ N. The v ertical structure of the mo del atmosphere consists of only tw o discrete la yers: this is the minimal v ertical resolution needed to represen t baro clinic pro cesses [38]. Fiv e v ertical pressure levels deﬁne the t wo la y ers with boundaries at: p 2 . 5 = 1000 hPa (surface lev el), p 2 = 750 hPa, p 1 . 5 = 500 hPa, p 1 = 250 hPa, p 0 . 5 = 0 hPa (top lev el). The geostrophic stream function ψ is deﬁned at lev els p 1 and p 2 , ψ ( p 1 ) = ψ 1 and ψ ( p 2 ) = ψ 2 , where the quasi-geostrophic vorticit y equation for the mid-latitude β -plane (7) - (8) is applied, while the vertical velocity ω is sp eciﬁed at lev el p 1 . 5 , where the thermo dynamic energy equation (9) is v alid. The mo del is describ ed by the following equations in terms of the barotropic stream function ψ M = ( ψ 1 + ψ 2 ) / 2, baro clinic stream function ψ T = ( ψ 1 − ψ 2 ) / 2, and temp erature T : ∂ ∂ t ( ∇ 2 ψ M ) = − J ( ψ M , ∇ 2 ψ M + β y ) − J ( ψ T , ∇ 2 ψ T ) − r ∇ 2 ( ψ M − ψ T ) + k h ∇ 4 ψ M , (7) ∂ ∂ t ( ∇ 2 ψ T ) = − J ( ψ T , ∇ 2 ψ M + β y ) − J ( ψ M , ∇ 2 ψ T ) + r ∇ 2 ( ψ M − ψ T ) + k h ∇ 4 ψ T + f 0 ∆ p ω , (8) 8 ∂ ∂ t ( T ) = − J ( ψ M , T ) + S p ω + r R ( T e − T ) + κ ∇ 2 T (9) In the ab ov e, we expressed the advection in terms of the Jacobian op erator deﬁned as J ( A, B ) = ∂ A ∂ x ∂ B ∂ y − ∂ A ∂ y ∂ B ∂ x . S p represen ts the static stability parameter [38]. W e deﬁne the stability parameter S = RS p ∆ p 2 f 2 0 = L 2 D , where L D is the Rossby radius of deformation. The name and v alues of model parameters are listed in T able 1. The vertical v elo cit y is set to 0 at the top level, ω 0 = 0, and is deﬁned through Ekman pumping at the surface lev el, ω 2 . 5 = ∆ p f 0 2 r ∇ 2 ψ 2 , whic h parameterises the dissipativ e pro cesses o ccurring in the boundary la yer. Subgrid-scale pro cesses are represented b y momen tum and heat diﬀusion terms. The system is driv en by a Newtonian co oling term that in volv es the restoration temp erature ﬁeld: T e = ∆ T 2 cos π y L y . (10) ∆ T denotes the forced meridional temp erature diﬀerence, and quan tiﬁes the external forcing in the model. In the p erformed simulations, no time-dep endence of ∆ T is assumed, with the aim of creating time series of a deterministic equiv alent of a stationary process. If ∆ T is suﬃcien tly large, the system reac hes a steady state featuring a turbulent atmospheric ﬂow with sensitiv e dep endence on initial conditions. The ph ysical processes responsible for limited predictabilit y are in general the baro clinic and barotropic instabilit y . The Newtonian co oling pro vides the so-called baroclinic forcing to the system, and activ ates a set of energy exc hanges and transformations summarised b y the framew ork of the Lorenz energy cycle. See discussion in Holton [38]. Using hydrostatic appro ximation [38], w e obtain for our vertical discretization T = 2 f 0 R ψ T . Th us, the three model equations (7) - (9) can b e reduced to tw o equations with t wo v ariables ψ M and ψ T . The mo del equations can be transformed [29] into a non-dimensional form using the scaling factors in T able 1. In the following, w e use non-dimensional quan tities, if not indicated otherwise. As mentioned, the channel is perio dic in the x -direction. A t the meridional b oundaries, we set the meridional v elo cit y and the zonally-integrated zonal v elo cit y to 0, v = ∂ ψ ∂ x = 0 and R 2 π /a 0 d x ∂ ψ ∂ y    y =0 ,π = 0. F or these b oundary conditions, the solution of the mo del equations is ψ ( x, y , t ) = N x ,N y X k,l =1  ψ r ( k , l , t )cos( ak x ) + ψ i ( k , l , t )sin( ak x )  sin( ly ) + N y X l =1 ψ r (0 , l, t )cos( l y ) , (11) with ψ denoting the non-dimensional stream function, index r represents the real and i the imaginary co eﬃcien ts. W e apply a sp ectral cut-oﬀ at N x = N y = 16 in b oth x and y directions. Hence, the total dimension of the model phase space is 2 N y (2 N x + 1) = 1056. 3.2 Metho ds W e substitute Eq. (11) in to the ev olution equations, p erform a Galerkin pro jection, and ev entually in tegrate n umerically the non-dimensional model equations in sp ectral space using the fourth-order Runge-Kutta scheme. W e p erform sim ulations with tw o diﬀerent forced 9 meridional temp erature diﬀerences, ∆ T = 133 K and ∆ T = 40 K. In the case of strong forcing (∆ T = 133 K), the system has a Kaplan-Y ork e (KY) dimension of 585.95 with 222 p ositiv e Lyapuno v exp onen ts, so that d u = 222, d n = 2, and d s = 361 . 95. Note that the presence of a second neutral direction is related to the existence of a rotational symmetry in the system and to the fact that we consider a spectral model. This feature is of little relev ance for the analysis b elo w. W e pro duce stationary time series of 96,576 years with a time step of 0.7 hours. In the case of w eak forcing (∆ T = 40 K), the system has a KY dimension of 39.31 with 17 p ositive Lyapuno v exp onen ts, so that d u = 17, d n = 2, and d s = 20 . 31. W e pro duce stationary time series of 485,760 y ears with a time step of 2.8 hours. The sp ectral co eﬃcien ts ψ r/i of the stream functions are recorded ev ery 5.5 hours with either forcings. The Ly apunov exp onen ts are obtained b y the same sim ulation co de as the one used by Sch ubert and Lucarini [29], based on the method of Benettin et al. [39]. T able 1: List of symbols and parameter v alues for the QG mo del, Eq. (7) – (11). V ariable Symbol Unit Scaling factor Stream function ψ m 2 s − 1 L 2 f 0 T emperature T K 2 f 2 0 L 2 /R V elocity v ms − 1 Lf 0 Energy e Jkg − 1 L 2 f 2 0 Par ameter Symbol Dimensional V alue Non-dimensional v alue Scaling F actor F orced meridional temp erature diﬀerence ∆ T 133 & 40 K 0.188 & 0.0564 2 f 2 0 L 2 /R Ekman friction r 2 . 2016 × 10 − 6 s − 1 0.022 f 0 Eddy-momentum diﬀusivity k h 10 5 m 2 s − 1 9 . 8696 × 10 − 5 L 2 f 0 Eddy-heat diﬀusivit y κ 10 5 m 2 s − 1 9 . 8696 × 10 − 5 L 2 f 0 Thermal damping r R 1 . 157 × 10 − 6 s − 1 0.011 f 0 Stability parameter S 3 . 33 × 10 11 & 2 . 52 × 10 11 m 2 0.0329 & 0.0247 L 2 Coriolis parameter f 0 10 − 4 s − 1 1 f 0 Beta ( d f /dt ) β 1 . 599 × 10 − 11 m − 1 s − 1 0.509 f 0 /L Aspect ratio a 0.6896 0.6896 - Meridional length L y 10 7 m π L Zonal length L x 2 . 9 × 10 7 m 2 π /a L Speciﬁc gas constant R 287 . 06 Jkg − 1 K − 1 2 R/ 2 V ertical pressure diﬀerence ∆ p 500 hP a 1 ∆ p Time scale t 10 4 s 1 1 /f 0 Length scale L 10 7 /π m 1 10 7 /π The sp ectral output of the mo del is transformed in to the grid p oin t space using F ast F ourier T ransform resulting n x × n y grid p oin ts with n x = n y = 36 in the x and y directions. W e refer to the grid p oin ts by indices ( i x , i y ) where i x = i y = 0 , ..., 35. W e analyse extremes of total energy observ ables deﬁned in non-dimensional form b elo w. F or our extreme v alue analysis we consider only the “mid-latitudes” of the QG mo del, whic h w e deﬁne as the region b et w een the latitudes i y = 9 and i y = 26, i.e., the latitudinally central 0.5 fraction of the whole domain. The total energy is the sum of the kinetic energy of the lo wer and upper lay ers and of the av ailable potential energy: e ( i x , i y , i t ) = e k 1 ( i x , i y , i t ) + e k 2 ( i x , i y , i t ) + e p ( i x , i y , i t ) , (12) where i t represen ts the discrete time co ordinate. The comp onen ts of the right side of (12) are deﬁned for each grid p oin t as: e k 1 ( i x , i y , i t ) = 1 2  u 1 ( i x , i y , i t ) 2 + v 1 ( i x , i y , i t ) 2  , (13) 10 e k 2 ( i x , i y , i t ) = 1 2  u 2 ( i x , i y , i t ) 2 + v 2 ( i x , i y , i t ) 2  , (14) e p ( i x , i y , i t ) = 2 λ 2 ψ T ( i x , i y , i t ) 2 , (15) with the zonal comp onen t of the horizontal v elo cit y u 1 = − ∂ ψ 1 ∂ y , u 2 = − ∂ ψ 2 ∂ y , the meridional comp onen t of the horizontal v elo cit y v 1 = ∂ ψ 1 ∂ x , v 2 = ∂ ψ 2 ∂ x , and λ 2 = 1 / (2 S ). W e obtain the zonally-a veraged energy b y taking the zonal av erage of the lo cal energy (Eq. (12)): e z ( i y , i t ) = 1 n x n x − 1 X i x =0 e ( i x , i y , i t ) , (16) and the a verage mid-latitude energy by a veraging the local energy o v er the area corresp onding to the mid-latitudes: e ml ( i t ) = 2 n x n y n x − 1 X i x =0 26 X i y =9 e ( i x , i y , i t ) . (17) The energy observ ables are analysed in their non-dimensional form. The ph ysical v alues expressed in units of J/Kg (J / m 2 ) can b e obtained b y m ultiplying the non-dimensional v alues b y the factor L 2 f 2 0 = 1 . 013 × 10 5 ( L 2 f 2 0 ∆ p/g = 5 . 164 × 10 9 ). Although we record the mo del output, as stated ab o v e, ev ery 5.5 hours, w e sav e only the maxim um v alues o ver one month in the case of strong forcing and o ver three months in the case of w eak forcing. W e estimate the GEV and GP parameters based on block maxima and threshold exceedances obtained from the mon thly , resp ectively 3-mon thly , maxima series. Suc h an operation has no eﬀect on the subsequent GEV analysis. Instead, it might mo destly impact the GP analysis, as some abov e threshold ev en ts migh t b e lost, because they could b e mask ed b y a larger even t occurring within the same 1-mon th or 3-months perio d. Nonetheless, since we consider very high thresholds and an extremely lo w fraction of ev ents, the risk of losing information is negligible. The GEV and GP parameters are inferred by maximum lik eliho o d estimation (MLE), as describ ed by Coles [11]. W e estimate the GEV and GP parameters, as well as the conﬁdence interv als, using the MA TLAB functions gevfit and gpfit . The computed conﬁdence in terv als con tain the true v alue of the parameters with a probabilit y of 95%. The auto-correlation co eﬃcien ts and histograms are obtained based on 1000 years of the “ra w” sim ulated time series. 4 Extreme v alue statistics in the QG mo del As mentioned before, the sim ulations are p erformed using tw o diﬀerent conﬁgurations, where the v alue of the parameter ∆ T , describing the baro clinic forcing, is set to 133 K and 40 K, resp ectiv ely . As w e see b elow, in the case of strong forcing we ﬁnd goo d agreemen t betw een the results of the statistical inference and the theory , for lo cal observ ables at least, ev en if the sp eed of conv ergence of the estimated shape parameters (not predicted by the theory), to the v alue that is predicted by the theory is rather diverse among the considered observ ables. In the case of weak forcing and the resulting weakly turbulen t b eha viour, the results of the statistical inference analysis are in not so go od agreement with the theory , and we ﬁnd 11 that for the diﬀeren t observ ables the shap e parameter estimates ha ve v arious non-monotonic dep endence on the blo c k size. W e will inv estigate p ossible reasons for such a b ehaviour. 4.1 Strong forcing ( ∆ T = 133 K) Before presen ting the results related to the statistics of extreme ev ents, w e outline some general statistical prop erties of the analysed observ ables. As emphasised in Sec. 1 and 2, correlations hav e an eﬀect on the con v ergence of the distribution of blo c k maxima (threshold exceedances) to the GEV (GP) distribution. The auto-correlation co eﬃcien t for a stationary signal f ( t ) can b e deﬁned as follows: ρ f ( l ) = E[( f ( t ) − µ f )( f ( t + l ) − µ f )] σ 2 f , (18) where t represen ts the time and l the time lag, µ f denotes the mean and σ f the v ariance of f ( t ) [40]. By taking the ergo dic hypothesis, we estimate the auto-correlation co eﬃcien t for the lo cal energy according to (18) and obtain for each grid p oint ρ e = ρ e ( i x , i y , l ). W e calculate the in tegrated auto-correlation time scale [41] τ e = τ e ( i x , i y ) according to τ e = P n l l =0 | ρ e | 1 . W e set n l = 604 (corresponding to about 140 da ys) as an upp er limit for the in tegration in order to av oid the noisy tail of the auto-correlation coeﬃcient. W e proceed the same w ay in the case of the zonally-a veraged and av erage mid-latitude energy to obtain τ e z ( i y ) = P n l l =0 | ρ e z ( i y , l ) | and τ e ml = P n l l =0 | ρ e ml ( l ) | . The in tegrated auto-correlation time scales, expressed in da ys, are shown in Fig. 1 (a). As exp ected, the weak est auto-correlations are recorded for the lo cal observ ables, yielding ab out 1-2 da ys. Because of the information propagation along parallels due to the prev ailing zonal winds, the zonal av erage time series are impacted b y a lo w-pass ﬁltering as result of a veraging along a latitudinal band, thus the correlations b ecome stronger. F or these zonally-a veraged observ ables, as opp osed to the lo cal ones, the in tegrated auto-correlation changes substantially with latitudes. W e observe a minim um in the middle of the channel ( ≈ 3 . 5 da ys) and an increase outw ards to the b oundaries ( ≈ 15 da ys). Through av eraging ov er the area of mid-latitudes, the zonally-a veraged time series with diﬀerent prop erties are merged together. The resulting time series has an integrated auto-correlation time scale of ab out 3 . 6 da ys. Note that if in a time series of length N with reasonably fast (e.g. similar to exp onen tial) decay of correlations the integrated auto- correlation time scale is τ , then one can deduce that the time series has approximately N × ∆ t/τ eﬀectively independent en tries, where ∆ t is the time interv al. This can hav e imp ortan t eﬀects in determining when the asymptotic b eha viour of the EVT statistics is v alid. Fig. 1 (b) – (d) illustrates the histograms and the appro ximated probabilit y densit y functions (PDFs) of our observ ables. Although, due to the β -eﬀect, the dynamical prop erties of the ﬂo w as a function of latitudes are not exactly symmetric with resp ect to the meridional middle of the c hannel, our estimations of statistical quan tities and density functions exhibit appro ximate meridional symmetries. Therefore, in the case of the lo cal and zonally-av eraged 1 If the deca y of the auto-correlation is exp onen tial, the integrated auto-correlation time scale is equal to the e-folding time. 12 observ ables only half of the c hannel’s meridional extension (at every second latitude) is shown. The strongest skewness and the longest right tails are observed in the case of the PDFs of the lo cal observ ables. After spatial av eraging, the PDFs b ecome more symmetric and almost similar to a Gaussian distribution (which would lo ok lik e a parab ola on a semi-logarithmic scale), according to the cen tral limit theorem. 9 11 13 15 17 19 21 23 25 0 5 10 15 20 Int. auto-corr. [days] a) b) 0 0.1 0.2 0.3 0.4 e 0 5 10 Density lat 9 lat 11 lat 13 lat 15 lat 17 c) 0.02 0.04 0.06 0.08 0.1 0.12 e z 0 10 20 30 Density d) 0.04 0.05 0.06 0.07 0.08 e ml 10 0 10 1 Density + 1 Figure 1: Statistical prop erties of the total energy for ∆ T = 133 K. (a) Integrated auto- correlation time scales: zonal a verage of τ e (dashed line with star mark ers), τ e z (dotted line with circle markers), τ e ml (con tinuous line); histograms of the (b) lo cal, (c) zonally-av eraged, and (d) av erage mid-latitude observ ables. In the case of (b) – (d), the contin uous lines show the approximation of the PDFs by k ernel smo othing ( ksdensity function in MA TLAB), the colours mark diﬀeren t latitudes according to the legend. In what follows, w e present the results of the extreme ev ents analysis starting with the local observ ables. W e ﬁrst discuss the con vergence of the shape parameter for GEV and GP , then the conv ergence of the GP mo diﬁed scale parameter (to b e introduced b elow), and, at the end, the conv ergence of return lev els. T aking adv antage of the fact that statistics are uniform in the zonal direction, w e concatenate the mon thly maxima series for ev ery second longitude one after the other in the x -direction, th us increasing the data length to about 1 . 7 × 10 6 (from ab out 9 . 6 × 10 4 ) years. Therefore, we can estimate the GEV and GP shap e parameters for larger blo c k sizes and higher thresholds than in the case of the zonally-av eraged or a verage mid-latitude observ ables. Although the time series at ev ery second longitude are correlated with eac h other, the correlation almost v anishes at the blo c k size of 8 y ears, b eing below 13 0.15 at every latitude. In other w ords, correlations are v ery w eak at extreme levels, which is according to Coles [11] the only imp ortan t condition for the GEV limit laws to apply in the case of a stationary pro cess. Blo c k sizes smaller than 8 years are not relev ant for our analysis, since (as presented below) m uch larger ones are needed to approach the theoretical shap e parameter. In the case of the POT approach, w e use the same argument of choosing v ery high thresholds, abov e which the correlations are extremely w eak. The theory discussed in Sec. 2 indicates that the true (asymptotic) GEV shap e parameter is given by ξ δ , as expressed by Eq. (5), whic h corresp onds to approximately -0.002, and is indicated b y the straigh t line in Fig. 2. Note that the range of the theoretical v alues deriv ed taking into consideration possible geometrical degeneracies, according to what is describ ed in Sec. 2, is to o small to be visible in this case. W e deﬁne the pr e cision P ( n ) of estimation by considering half of the width of the 95 % maxim um lik eliho o d conﬁdence interv al. It is more common to deﬁne the precision b y the standard deviation of the estimates. F or a Gaussian distribution the 95 % conﬁdence interv al is larger than the standard deviation by a factor of approximately 2. Ho wev er, this distinction do es not matter for our purp oses. Besides, w e hav e a single estimate only , and so we can obtain only the conﬁdence interv al not the standard deviation of the distribution of estimates. Additionally to the precision, we deﬁne the trueness of a single estimate by the distance b et w een ξ δ and ξ ( n ): T ( n ) = | ξ δ − ξ ( n ) | . Note that the latter is diﬀerent from the usual deﬁnition in that the reference from whic h w e measure the distance is ξ δ , not the true v alue of the distribution from which the BM data is dra wn. In fact, strictly , the BM data is not drawn from a GEV distribution that we are ﬁtting, and hence we cannot even talk ab out the true v alue of an underlying GEV distribution. W e emphasise that our in terest is the con v ergence to the asymptotic v alue, whic h is why we take a reference v alue in our deﬁnition other than customary . Accordingly , we shall refer to the ‘bias’ of the estimator, again diﬀeren t from customary , as the exp ected trueness. W e remark that, since our estimates are obtained based on one realisation instead of several realisation yielding a distribution of estimates, our trueness T approximates the bias of the estimates as long as P  T . W e emphasise that we are able to calculate T here because w e know the true ξ δ ; this is not the case in practice when facing just a measured time series. Ob viously , we aim at obtaining a joint optimisation by ha ving a bias and a precision as small as pos sible. Clearly , optimalit y requires a compromise b et w een these tw o requiremen ts. When we apply the BM method and increase the block size n , the num ber of blo c ks and of BM decreases, th us the estimation of ξ ( n ) becomes more and more uncertain, and P increases monotonically . At the same time, for increasing n we exp ect a (not necessarily monotonic) con vergence of our estimated shap e parameter to the true v alue, so that the actual bias should (on the long run) decrease with n . Clearly , instead, our appro ximation T decreases only un til a certain blo c k size, ab o ve which it b ecomes more uncertain with increasing v alues of n , b ecause less BM are av ailable. W e choose as optimal blo c k size n = n ∗ the smallest blo ck size for which the estimate of the bias is low er than the estimate of the precision n ∗ = min( n ; T ( n ) < P ( n )). On the scale of v ariation ranges of P ( n ) and T ( n ) w e hav e T ( n ∗ ) ≈ P ( n ∗ ). With this w e obtain a single n um b er that can quantify the ac cur acy of estimation. This measure of accuracy provides here a basis for comparing diﬀeren t observ ables with regard to the sp eed of conv ergence, or a basis for assessing the degree of non-uniformity of estimates of v arious observ ables of interest (in terms of the range of accuracy v alues), as a ﬁnite-data-size deviation form the uniformity predicted b y theory . W e ha ve v eriﬁed that the optimal c hoice for n ∗ is virtually unc hanged when we use an alternative deﬁnition of the accuracy suc h as T 2 + P 2 , b orro wing an idea concerning the optimality of MLE estimators (not sho wn). 14 First w e assess the uniformity for the lo cal observ ables. Figure 2 (a) shows the GEV shap e parameter estimates against exp onen tially increasing blo c k sizes of n = 2 i y ears ( i = − 2 , − 1 , ..., 13), for diﬀerent latitudes. The estimated GEV shap e parameters ξ ( n ) seem to con verge monotonically for every latitude to ξ δ . The monotonic conv ergence is p oin ted out also in panel (b) in terms of T ( n ). In this diagram we displa y P ( n ) to o, b y whic h w e can determine the optimal n ∗ and the accuracies of estimation. These accuracies, dep ending on the latitude, hav e a range of 5 × 10 − 3 − 2 × 10 − 2 . A t the same time, the v alue of n ∗ ranges from sev eral tens of years to a few h undreds of y ears depending on the latitude we are considering. This is unsurprising, because the sp eed of conv ergence to the asymptotic lev el is not univ ersal. As a consequence, when ﬁnite blo ck sizes are considered, extremes of diﬀeren t observ ables can feature rather distinct prop erties. The slo w conv ergence suggests that cus- tomary choices lik e y early maxima are not alw a ys go od enough for an accurate mo delling of extremes. Figure 2 (c), giving a diﬀeren t view of the same data seen in panel (a), illustrates the estimated GEV shap e parameter as a function of latitudes for v arious blo c k sizes. F or small blo c k sizes, we observe a slight latitudinal dep endence of the shap e parameter. This latitudinal structure ﬂattens as one increases the blo c k size, and the estimated shape param- eters get closer to the theoretical v alue. According to Fig. 2 (c), universalit y emerges as w e approac h the asymptotic lev el. T o asses s the go odness of ﬁt, we p erform a one-sample Kolmogoro v-Smirnov-test (KS-test) [12] at 5 % signiﬁcance lev el using the MA TLAB function kstest . W e remark that the KS-test is p erformed in case of each blo c k size based on the whole BM data, meaning that this amount of data decreases as w e increase the blo c k size. The shap e parameter v alues for which the KS-test p-v alue p is ab o ve 0.05 (i.e., the hypothesis that the distribution of BM is a GEV distribution cannot b e rejected) are marked by circle mark ers in Fig. 2 (a). W e deﬁne n K S as the smallest blo ck size for whic h p > 0 . 05, n K S = min( n ; p > 0 . 05), and ξ K S = ξ ( n K S ). Figure 2 p oin ts out that the KS-test suggests a go od ﬁt already at smaller blo c k sizes than the optimal block size, n K S < n ∗ , and for lo w er shap e parameter v alues than the b est estimate, ξ K S < ξ ( n ∗ ). Th us, a very imp ortant conclusion is that the p-v alue of the KS-test is not an appropriate measure for the conv ergence to the limiting distribution. More precisely , it indicates that w e ha ve indeed agreement with a mem b er of the GEV family of distributions, but we cannot say what is the error from the asymptotic v alue of the parameters. W e emphasise that n K S , just lik e n ∗ , dep ends on the time series length, and it w ould b e ev en smaller if shorter time series w ere considered. This implies that in the case of applications with less data, the results of the KS-test are ev en less reliable. The misleading prop ert y of p-v alues was also sho wn by B´ odai [25], who studied the conv ergence to the GEV distribution of extremes of site v ariables in the Lorenz 96 mo del, and found p-v alues ab o v e the signiﬁcance level in cases where the theoretical prediction did not even apply , and the shap e parameter did not con v erge. The go odness-of-ﬁt test was in the men tioned study a P earson’s c hi-squared test. Misleading p-v alues based on the KS-test were p oin ted out also b y F aranda et al. [8] in the case of the BM approac h in simple systems. A slo w con v ergence of the estimated GEV shap e parameters and a po or qualit y of diagnostic tools (return lev el and quan tile plots) in case of small blo ck sizes were also found b y V annitsem [27] in case of lo cal temp erature extremes in a three-lay er QG mo del with orography. Figure 2 (d) illustrates the GP shap e parameter estimates as a function of decreasing ex- ceedance ratio (the fraction of abov e-threshold data) r , whic h is equiv alent with an increasing threshold. T o ensure direct comparability b et ween the BM and POT approac hes of EVT, 15 sample v alues of the threshold are chosen corresp onding to the sam ple v alues of the blo c k size, in such a w ay that r = 1 n × m y , where m y is the data amoun t in a year. Thus the n umber of threshold exceedances is equal to the num b er of blo c k maxima. By comparing the GP shap e parameter (Fig. 2, d – f ) with the GEV shape parameter (Fig. 2, a – c), w e generally observ e the same characteristics. More precisely , the changes of the GP shap e parameter as a function of exp onen tially decreasing exceedance ratio are very similar to the v ariation of the GEV shap e parameter according to an exp onen tially increasing blo c k size. Both GEV and GP shap e parameters seem to conv erge to ξ δ . This is also consistent with theoretical results according to which the tw o distributions are asymptotically equiv alen t [11, 22]. How- ev er, w e exp ect that in case of ﬁnite blo c k sizes (i.e., in case of every practical application) diﬀerences might emerge in the estimates of the GEV and GP shap e parameters. Although, in the case of consisten t estimations one would exp ect that at large blo c k sizes, corresp ond- ing to low exceedance ratios, the diﬀerence b etw een them should b e small, as it is the case for our estimators. Besides the mentioned similarities, we observe some diﬀerences b et w een the estimates of the GEV and GP shap e parameters. These diﬀerences concern, for exam- ple, their latitudinal dep endence (less pronounced in case of the GP shape parameter) or the width of the conﬁdence interv als (larger in case of the GP shap e parameter, indicating larger estimation uncertaint y). The most relev an t diﬀerence is, how ev er, that the GP shap e parameter seems to conv erge faster to ξ δ . This is unsurprising as in man y applications it is usually suggested to use the POT o ver the BM metho d as the former is less data-hungry and pro vides (usually) a faster conv ergence [7]. W e p erform another test to c heck whether the GP distribution is a go o d appro ximation for the distribution of threshold exceedances based on our data, and consider the GP mo diﬁed scale parameter. The GP scale parameter depends on the c hosen threshold according to σ u = σ u 0 + ξ ( u − u 0 ) [11], where ξ = ξ δ represen ts the asymptotic shape parameter, u 0 is the lo west threshold at which the GP distribution is a reasonable mo del for exceedances, and u represen ts any other threshold u > u 0 . The scale parameter can b e reparameterized yielding the mo diﬁed scale ˆ σ = σ u 0 − ξ u 0 = σ u − ξ u , whic h should con verge to a nonzero v alue. Figure 3 illustrates the modiﬁed scale parameter estimate (calculated based on the ﬁnite-size GP parameter estimates, i.e., taking threshold dep endent GP shap e parameter estimates instead of ξ δ ) as a function of the exceedance ratio r . W e observ e estimates of ˆ σ relativ ely stable to further decreases of r (for r < r ∗ , r ∗ = max( r ; T ( r ) < P ( r ))), whic h supp orts the conclusions dra wn before that w e are indeed close to asymptotic lev els as required by EVT. Note that in this case there is no univ ersality in the v alue of the mo diﬁed scale parameter, as for stochastic v ariables one has that the upp er right endpoint of the distribution is giv en b y A max = − σ u 0 /ξ + u 0 = − ˆ σ /ξ . Such an endp oint is clearly observ able-sp eciﬁc. Ha ving practical applications in mind, the BM and POT metho ds aim at obtaining statistical estimates of either return lev els or expected return p eriods, for ev en unobserved extreme ev ents. Figure 4 (a) and (b) show GEV and GP return level plots for the lo cal observ ables based on a ﬁxed blo c k size, n = 128 years, and corresp onding (as explained) r = 5 × 10 − 4 %, resp ectiv ely , at ﬁv e diﬀeren t latitudes (ev ery second latitude from the southern meridional b oundary to the channel cen tre). W e compute the GEV return levels according to Eq. (3), the GP return lev els based on Eq. (4), and estimate the 95 % conﬁdence in terv als using the delta method describ ed b y Coles [11]. The GEV and GP return level plots lo ok very similar, except tw o minor diﬀerences. One emerges simply from the diﬀeren t equations for the GEV and GP distributions, leading to slightly diﬀerent deﬁnitions of return lev els (as described in 16 block size [yr] -0.2 -0.1 0 0.1 0.2 GEV sh.p. a) lat 9 lat 10 lat 11 lat 12 lat 13 lat 14 lat 15 lat 16 lat 17 lat 18 lat 19 lat 20 lat 21 lat 22 lat 23 lat 24 lat 25 lat 26 10 0 10 1 10 2 10 3 block size [yr] 10 -4 10 -2 10 0 c.i. width exceedance ratio [%] -0.2 -0.1 0 0.1 0.2 GP sh.p. d) 10 -5 10 -4 10 -3 10 -2 10 -1 exceedance ratio [%] 10 -4 10 -2 10 0 c.i. width 10 0 10 1 10 2 10 3 block size [yr] 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 GEV sh.p. bias and prec. b) 10 -5 10 -4 10 -3 10 -2 10 -1 exceedance ratio [%] 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 GP sh.p. bias and prec. e) 9 11 13 15 17 19 21 23 25 -0.1 0 0.1 GEV shape p. c) 9 11 13 15 17 19 21 23 25 -0.1 0 0.1 GP shape p. f) 0.25 yr 0.5 yr 1 yr 2 yr 4 yr 8 yr 16 yr 32 yr 64 yr 128 yr 256 yr 512 yr Figure 2: GEV and GP shape parameters as w ell as bias and precision estimates in case of the lo cal observ ables, for ∆ T = 133 K. (a) GEV and (d) GP shap e parameter estimates as functions of the blo c k size and exceedance ratio. The circle markers p oin t out shape parameter v alues, for which the p-v alue of the KS-test is ab ov e 0.05. Low er y -axis: zonal mean of the maximum likelihoo d 95 % conﬁdence interv al widths. (b) and (e): Estimates for the bias (dot mark ers) and precision (star markers) of the shape parameter. (c) GEV and (f ) GP shap e parameter estimates as functions of the latitude. The grey , horizon tal line illustrates the theoretical shap e parameter of -0.002. The range of theoretical v alues resulting from taking in to consideration p ossible geometrical degeneracies is invisible in this case. The error bars sho w the 95 % conﬁdence in terv als of the MLE. Diﬀeren t colours represen t diﬀerent latitudes (a, b, d, e) or diﬀeren t blo ck sizes (c) or exceedance ratios (f ). Sec. 2 and in more details in 11) and aﬀects short return p eriods; and the other one comes from the larger uncertaint y in the estimation of the GP parameters compared to the GEV parameters, and results slightly wider conﬁdence in terv als in the case of the GP return lev els. The main message of Fig. 4 (a) – (b) is, ho w ev er, that the GEV and GP return lev el estimates using the chosen n and r ﬁt the empirical data quite well, whic h is in agreement with the results of the KS-test rep orted ab o v e. The 95 % conﬁdence interv als of the estimated return lev els (contin uous lines) contain the empirical return lev els (dot mark ers) or are very near to them, except a few very high extremes at some latitudes. The return lev el is almost linear to the logarithm of the return p erio d, showing the eﬀect of a shap e parameter very close to 0 (see Eq. (3) and (4)). If the GEV distribution is an adequate mo del for extreme even ts for a certain blo c k size, 17 10 -5 10 -4 10 -3 10 -2 10 -1 exceedance ratio [%] -0.05 0 0.05 0.1 GP mod. scale p. Figure 3: GP mo diﬁed scale parameter estimates in case of the lo cal observ ables, for ∆ T = 133 K. The error bars sho w the 95 % conﬁdence interv als of the maximum lik eli- ho od estimation. Diﬀeren t colours represent diﬀeren t latitudes according to Fig. 2. . one exp ects return levels with a certain return perio d not to c hange muc h any more with increasing blo c k size. Figure 4 (c) – (e) shows indeed that, ab o v e a certain blo c k size, the estimated return lev els for three diﬀeren t return p eriods (10 3 , 10 4 , and 10 5 y ears) are stable against further increase of n . But it also shows that the longer the return p eriod, the slo wer the conv ergence. While in case of the 10 3 -y ears return p erio d we obtain stable return lev el estimates already at n K S , in case of 10 5 -y ears the return lev el estimates are still increasing for n > n K S . Here we experience the practical eﬀect of the issue men tioned ab ov e, namely , that the KS-test suggests a go o d ﬁt ev en for ξ K S < ξ ( n ∗ ). This implies that the estimation of return lev els with long return perio ds can b e erroneous even if the KS-test do es not reject the GEV distribution. W e also notice that the return levels are underestimated if the blo c k size is to o small, and this underestimation is more severe in the case of return levels with longer return p erio ds. W e come to the same conclusion by considering the con vergence of the GP return levels (not sho wn), as suggested already by the similarity betw een panels (a) and (b). After having discussed in detail the con vergence in case of the lo cal observ ables, w e pro ceed with the results for the zonally-av eraged observ ables. Figure 5 illustrates the GEV and GP shap e parameters for the zonally-av eraged observ ables (a, c, d, f ) as w ell as the estimated bias and precision of the inferred shap e parameters (b, e). As men tioned ab ov e, in case of the zonally-av eraged observ ables we ha v e shorter time series (9 . 6 × 10 4 instead of 1 . 7 × 10 6 y ears). Because of this, results for the accuracies of estimates cannot b e ‘fairly’ compared to the accuracies found for lo cal observ ables. Nev ertheless, w e pro duce the same t yp e of diagrams suitable to determine the accuracies and show it in Fig. 5 (b) and (e). Clearly , the range of accuracy v alues dep ending on the latitude and the maximal v alue of the accuracies (i.e, of the bias at the optimal blo c k size) are b oth considerably larger than those for the lo cal observ ables. What is fair to compare, ho wev er, is the range of biases for a certain blo c k size where the conﬁdence of the estimates is high, P  T , and the amount of data do es not aﬀect signiﬁcan tly the parameter estimate. In this regard, the zonal observ ables 18 Figure 4: Return lev els for ∆ T = 133 K. (a) Return lev els vs. return p erio ds based on GEV parameters using a block size of 128 y ears and (b) based on GP parameters using an exceedance ratio of 5 × 10 − 4 %. (Dotted lines: estimated return lev els; con tinuous lines: 95 % maximum lik eliho od conﬁdence in terv al limits of the return level estimates; dot mark ers: empirical return levels.) GEV Return levels for (c) 10 3 -y ears, (d) 10 4 -y ears, and (f ) 10 5 -y ears return perio ds as functions of the blo c k size. The error bars sho w the 95 % conﬁdence in terv als of the MLE. The circle markers p oin t out estimates for which the p-v alue of the KS-test is ab o v e 0.05. The colours mark diﬀerent latitudes according to the legend. displa y a muc h larger non-uniformity regarding the shap e parameter estimates. Otherwise, the estimates feature typically a monotonic change to wards the theoretical v alue (up to at least the optimal blo c k size), what can b e seen as con vergence. Our observ ation that the estimated shap e parameters dep end strongly on the considered latitude has to do with the eﬀect of serial correlation on the con vergence to the limiting distribution. W e obtain weak auto-correlations, fast conv ergence to ξ δ , and lo w bias in the middle of the channel, versus strong auto-correlations, slow con vergence, and large bias at the margins of the channel. As already mentioned b efore, the stronger the serial correlation the less the num ber of uncorrelated data in a block, and the larger blo c k sizes are needed in order to approac h asymptotic levels (see also [11]). Th us the latitudinal structure of the GEV shap e parameter estimates (Fig. 5, c) is related to the one of the integrated auto-correlation time scale (Fig. 1, a, dotted line with circle markers). By increasing the blo c k size, this latitudinal structure ﬂattens, the estimated shape parameters seem to approach ξ δ , and the conﬁdence interv als con tain ξ δ in the case of sev eral latitudes, esp ecially the cen tral ones. 19 Nonetheless w e note that, due to the presence of (relatively) large statistical uncertain ty on the shap e parameter, w e cannot mak e more precise statemen ts on the success of the analysis. block size [yr] -0.2 -0.1 0 0.1 0.2 GEV sh.p. a) lat 9 lat 10 lat 11 lat 12 lat 13 lat 14 lat 15 lat 16 lat 17 lat 18 lat 19 lat 20 lat 21 lat 22 lat 23 lat 24 lat 25 lat 26 10 0 10 1 10 2 block size [yr] 10 -4 10 -2 10 0 c.i. width 10 -3 10 -2 10 -1 exceedance ratio [%] -0.2 -0.1 0 0.1 0.2 GP sh.p. d) 10 -3 10 -2 10 -1 exceedance ratio [%] 10 -4 10 -2 10 0 c.i. width 10 0 10 1 10 2 block size [yr] 10 -4 10 -3 10 -2 10 -1 10 0 GEV sh.p. bias and prec. b) 10 -3 10 -2 10 -1 exceedance ratio [%] 10 -4 10 -3 10 -2 10 -1 10 0 GP sh.p. bias and prec. e) 9 11 13 15 17 19 21 23 25 -0.2 -0.1 0 0.1 GEV shape p. c) 9 11 13 15 17 19 21 23 25 -0.2 -0.1 0 0.1 GP shape p. f) 0.25 yr 0.5 yr 1 yr 2 yr 4 yr 8 yr 16 yr 32 yr 64 yr Figure 5: Same as Fig. 2, but for the zonally-av eraged observ ables. W e presen t now the analysis of extremes of the av erage mid-latitude observ able. Figure 6 sho ws the GEV and GP shap e parameter estimates for the a verage mid-latitude observ able and their estimated bias and precision as a function of the blo c k size and exceedance ratio, resp ectiv ely . In case of the av erage mid-latitude energy , w e ha ve the same amount of data as in case of the zonally-a veraged energy . Similarly to the zonally-av eraged observ ables, the estimated GEV and GP shap e parameters seem to approac h the theoretical shape parameter, but, when more stringent deﬁnitions for selecting the extremes are used, the bias is relativ ely large, being ab out 4 × 10 − 2 at the optimal blo c k size in the case of the GEV and ab out 7 × 10 − 2 at the optimal exceedance ratio in the case of the GP shap e parameter. Again, also in this case, our analysis is limited by the amount of av ailable data. In short, our n umerical results do allo w for conclusions regarding the univ ersality of extremes, as predicted b y the theory presen ted in Sec. 2. Ho wev er, considering the most v arious observ- ables one would typically see a non-uniformity in the ﬁnite-size shap e parameter estimates simply b ecause of their distinct con vergence prop erties (not predicted by the theory). The observ ables that we found in our study to ha ve the fastest con verging shap e parameter es- timates are the the lo cal observ ables at every latitude and the zonally-a veraged observ ables at cen tral latitudes, where the auto-correlation has a minimum. How ev er, con vergence is v ery slow, and is additionally slo wed down by the presence of serial correlations in the time series. Th us, the estimated shap e parameters are relatively far from the theoretical v alue 20 10 0 10 1 10 2 block size [yr] -0.3 -0.2 -0.1 0 0.1 0.2 GEV shape p. a) 10 -3 10 -2 10 -1 exceedance ratio [%] -0.3 -0.2 -0.1 0 0.1 0.2 GP shape p. c) 10 0 10 1 10 2 block size [yr] 10 -4 10 -3 10 -2 10 -1 10 0 GEV sh.p. bias and prec. b) 10 -3 10 -2 10 -1 exceedance ratio [%] 10 -4 10 -3 10 -2 10 -1 10 0 GP sh.p. bias and prec. d) Figure 6: GEV and GP shape parameter as well as bias and precision estimates in case of the a verage mid-latitude observ able, for ∆ T = 133 K. (a) GEV and (c) GP shap e parameters as functions of blo c k size and exceedance ratio. The circle markers point out shap e parameter v alues for whic h the p-v alue of the KS-test is ab o ve 0.05. The grey , horizontal line illustrates the theoretical shap e parameter of -0.002. The range of theoretical v alues resulting from taking into consideration p ossible geometrical degeneracies is invisible in this case. The error bars sho w the 95 % conﬁdence in terv als of the MLE. (b) and (d): Estimates of the bias (dot mark ers) and precision (star markers) of the shape parameter. in case of sev eral latitudes of the zonally-a veraged observ ables (esp ecially marginal latitude exhibiting strong auto-correlations) and in case of the a verage mid-latitude observ able. This slo w con v ergence in com bination with the ﬁnite size of the data makes the actual observ ation of the theoretical limit extremely diﬃcult. 4.2 W eak forcing ( ∆ T = 40 K) Before analysing the extreme ev en ts for w eak forcing, w e discuss some statistical (and dynam- ical) prop erties of our observ ables, whic h inﬂuence directly the statistics of extremes. Figure 7 (a) shows the in tegrated auto-correlation time scales for the three observ ables: lo cal, zonally- a veraged, and av erage mid-latitude energy . W e compute the integrated auto-correlation time scale according to the metho d describ ed in Sec. 4.1 for the strong forcing. In the case of w eak forcing, how ever, w e set the time lag n l = 1728 (corresp onding to ab out 400 da ys) as upp er limit for the in tegration, according to the slo w decay of the auto-correlation (especially in case of the zonally-av eraged and a verage mid-latitude observ ables). The in tegrated auto- correlation time scales are substan tially higher than for strong forcing: around 10 da ys in 21 case of the local, ab out 30 – 48 da ys in case of the zonally-av eraged observ ables, and approx- imately 45 days for the a verage mid-latitude observ able. Figure 7 (b) shows the time series of the lo cal observ ables at the central latitude i y = 17 (at tw o diﬀeren t longitudes i x = 4 and i x = 19) and suggests tw o alternating states of our system: one with strong ﬂuctuations and another one with reduced ﬂuctuations. Th us, it seems that our system exhibits a regime b eha viour, which deﬁnitely supp orts the presence of strong correlations. In contrast to the case of strong forcing, the zonal av erages of the lo cal energy observ ables sho w remark able deviations from a Gaussian b eha viour, ev en more than the PDFs of the lo cal energy observ ables (Fig. 7, c – e). One has that the PDFs of the zonally-av eraged ob- serv ables typically ha v e a marked skewness and very strong large kurtosis, and often con tain rather pronounced “shoulders”, where smo othness is basically lost. The presence of large kurtosis indicates that there is signiﬁcant p ositiv e spatial correlation of the energy along a longitude. The presence of skewness indicates that there is asymmetry b et w een the o ccur- rence of anomalies of either sign. Another particular prop ert y of the spatial energy ﬁeld for w eak forcing is the strong anti-correlation (especially in case of the zonally-av eraged observ- ables) b et w een time series at central and marginal latitudes (not sho wn). Accordingly , the “shoulders” app ear in diﬀeren t parts of the PDFs at diﬀeren t latitudes: on the left in the case of central latitudes and on the righ t in the case of marginal latitudes. W e conclude that the regime b eha viour is connected to non-trivial spatial structures, with the system living in a transitional range where one can still distinguish long-liv ed unstable w av es amidst chaos. W e note that such conditions are diﬀerent from what is foreseen by the chaotic h yp othesis, and, therefore, the statistics of extremes migh t not conv erge (according to our ﬁnite-sized data set) to what is predicted b y the theory developed for Axiom A systems. F or the analysis of extreme even ts, we use a similar pro cedure as in the case of strong forcing (∆ T = 133 K), and concatenate the three-mon thly maxima series for every second longitude one after the other in x -direction. Thus, we increase the length of a v ailable data for the lo cal observ ables to about 8 . 7 × 10 6 (from ab out 4 . 8 × 10 5 ) y ears. Although the time series at ev ery second longitude are correlated with each other, the correlation almost v anishes at extreme lev els, b eing b elo w 0.1 for every latitude in the case of the 8-years BM. W e deﬁne the GP exceedance ratios so that the num b er of threshold exceedances corresponds to the num b er of blo c k maxima, as describ ed in Sec. 4.1. In the case of weak forcing, the theoretical shap e parameter is -0.03, sho wn by the grey hor- izon tal line in Fig. 8. The grey shading represents the range of theoretical v alues resulting from taking into consideration p ossible geometrical degeneracies according to the limits de- scrib ed in Sec. 2. W e plot the GEV shape parameter against exp onen tially increasing blo c k sizes of n = 2 i y ears, where i = 0 , ..., 15 for the lo cal observ ables and i = 0 , ..., 11 for the zonally-a veraged and a verage mid-latitude observ ables. F o cusing ﬁrst on the lo cal observ- ables, w e notice a non-monotonic c hange of the shape parameter according to increasing block sizes. F or blo c k sizes smaller than 30 years, the shap e parameter ev en reac hes non-ph ysical, p ositiv e v alues for certain latitudes. This change of sign of the estimated shap e parameters is similar to what has b een observed b y V annitsem [27] in case of lo cal temp erature extremes in a more realistic QG mo del with orograph y . The non-monotonic c hanges and the p ositiv e shap e parameter estimates ha ve to do with the fact that, if the blo c k size is not large enough, w e select ev ents from b oth regimes (more and less ﬂuctuating) thus “con taminating” the statistics of extremes; whereas if the blo c k size is large enough, only extremes from the more unstable regime are selected. Figure 8 (a) also sho ws that the estimated shap e parameter 22 0 5 10 15 Time [mon] 2 4 6 8 10 12 14 e 10 -3 b) 9 11 13 15 17 19 21 23 25 0 10 20 30 40 50 Int. auto-corr. [days] a) c) 0 0.01 0.02 e 0 100 200 Density lat 9 lat 11 lat 13 lat 15 lat 17 d) 2 4 6 8 10 12 e z 10 -3 0 500 1000 Density e) 6.2 6.4 6.6 6.8 7 7.2 7.4 e ml 10 -3 10 0 10 1 10 2 10 3 Density + 1 Figure 7: Statistical prop erties of the total energy for ∆ T = 40 K. (a) Integrated auto- correlation time scales for the lo cal (dotted lines with star mark ers), zonally-av eraged (dashed lines with circle markers), and a verage mid-latitude (con tinuous line) observ ables. (b) Time series of the local energy at latitude i y = 17 and at t w o diﬀeren t longitudes: i x = 4 (red line) and i x = 19 (blac k line). Histograms of the (c) lo cal, (d) zonally-av eraged, and (e) av erage mid-latitude observ ables. In the case of (c) – (e), the con tinuous lines sho w the approximation of the PDFs b y k ernel smoothing ( ksdensity function in MA TLAB); the colours mark diﬀeren t latitudes according to the legend. seems to con v erge at almost ev ery latitude to a v alue which is low er than the theoretical shap e parameter, yet near to the range of v alues obtained taking into consideration p ossible geometrical degeneracies, see Sec. 2. As discussed ab ov e, this is in fact unsurprising giv en the qualitative prop erties of the system in the low forcing regime. In case of the zonally-av eraged and av erage mid-latitude observ ables we cannot detect any con vergence. This is an exp ected result, considering the statistical and dynamical character- istics of our data and the fact that the length of the time series is in this case even shorter than for the lo cal observ ables. As an eﬀect of the “shoulders” in the PDFs, we obtain very uncertain estimates ev en for large block sizes, and the KS-tests reject the h yp othesis of a GEV mo del in these cases. The shap e parameter estimates ha ve a large latitudinal spread due to the v arying form of PDFs according to latitudes. Except the diﬀerences b et ween the GEV (Fig. 8, a – c) and GP (8, d – f ) shap e parameters at small blo c k sizes and high exceedance ratios , b oth metho ds sho w us basically the same picture. The misleading prop ert y of the KS-test p-v alues p is underlined by Fig. 8. Ev en in case of the zonally-a veraged and av erage mid-latitude observ ables, where we cannot detect any con vergence at all, we ﬁnd p > 0 . 05 for a wide range of blo c k sizes and exceedance ratios (circle mark ers). 23 10 0 10 1 10 2 10 3 10 4 block size [yr] -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 GEV shape p. a) lat 9 lat 10 lat 11 lat 12 lat 13 lat 14 lat 15 lat 16 lat 17 lat 18 lat 19 lat 20 lat 21 lat 22 lat 23 lat 24 lat 25 lat 26 10 0 10 1 10 2 10 3 10 4 block size [yr] -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 GEV shape p. b) 10 0 10 1 10 2 10 3 10 4 block size [yr] -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 GEV shape p. c) 10 -5 10 -4 10 -3 10 -2 exceedance ratio [%] -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 GP shape p. d) 10 -5 10 -4 10 -3 10 -2 exceedance ratio [%] -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 GP shape p. e) 10 -5 10 -4 10 -3 10 -2 exceedance ratio [%] -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 GP shape p. f) Figure 8: Shap e parameter for ∆ T = 40 K. GEV shap e parameter for (a) lo cal, (b) zonally- a veraged, and (c) a verage mid-latitude energy . GP shap e parameter for (d) local, (e) zonally- a veraged, and (f ) av erage mid-latitude energy . The circle mark ers p oin t out shape parameter v alues, for whic h the p-v alue of the KS-test is ab o v e 0.05. Grey , horizontal line: theoretical shap e parameter. The grey shading represen ts a p ossible range of the theoretical shap e parameter according to the limits describ ed in Sec. 2. The error bars show 95 % conﬁdence in terv als of the maximum lik eliho od estimation. Diﬀerent colours represen t diﬀerent latitudes 5 Summary and discussion In this pap er we ha v e studied the con vergence of statistically estimated GEV and GP shap e parameters to a a theoretical shap e parameter. The latter is calculated based on prop erties of the attractor [21, 22]. W e analyse a quasi-geostrophic 2-la yer atmospheric mo del. W e study extremes of diﬀerent t yp es of energy observ ables: lo cal, zonally-av eraged, and av erage mid- latitude energy . W e p erform simulations with tw o diﬀerent forcing levels: a strong forcing (∆ T = 133 K), pro ducing a highly c haotic b eha viour of the system, and a w eak forcing (∆ T = 40 K), pro ducing a less pronounced c haotic be ha viour. In the case of strong (w eak) forcing, w e produce time series of about 9 . 6 × 10 4 (4 . 8 × 10 5 ) y ears, represen ting a deterministic equiv alen t to a stationary pro cess. W e estimate the GEV and GP shap e parameters for exp onen tially increasing blo c k sizes and exponentially decreasing exceedance ratios (fractions of ab o ve-threshold ev ents), i.e., increasing thresholds, b y p erforming maximum likelihoo d estimation. F or comparabilit y , we choose the GP thresholds so that the num b er of threshold exceedances corresp onds to the n umber of block maxima. W e tak e adv an tage of the fact that statistics are uniform in the zonal direction, and use the data from every second longitude for the analysis of extreme ev ents, thus increasing the length of av ailable data for the lo cal observ ables to ab out 1 . 7 × 10 6 (8 . 7 × 10 6 ) years in the case of strong (weak) forcing. W e start the discussion of our results with the strong forcing regime. In this case, w e observe 24 a roughly monotonic increase of the estimated GEV (GP) shap e parameters tow ards the theoretical v alue ξ δ = − 0 . 002. The estimated shap e parameters seem to conv erge to ξ δ in case of the lo cal observ ables at every latitude and in case of the zonally-av eraged observ ables at central latitudes. Thus, our numerical results allo w for robust conclusions regarding the univ ersality of extremes, according to the theory presen ted in Sec. 2. Ho wev er, in the case of several (esp ecially marginal) latitudes of the zonally-a veraged observ ables, as w ell as for the av erage mid-latitude observ able, the estimated shap e parameter is relativ ely far from the theoretical one. F or these observ ables the amount of data seems to b e not enough to approac h asymptotic lev els, thus we cannot mak e more precise statemen ts on the success of the analysis. Ev en in this extremely chaotic case, the conv ergence is very slo w, suggesting that customary choices like yearly maxima are not alwa ys the b est option for an accurate mo delling of extremes. Despite the predicted universal asymptotic prop erties of extremes, if we consider a certain blo c k size (threshold), we ﬁnd that the shap e parameter estimates are diﬀerent among the observ ables and latitudes. Th us, on pre-asymptotic level, extremes sho w rather div erse prop- erties. The speed of conv ergence to the asymptotic lev el is not univ ersal. The local observ ables exhibit high-frequency ﬂuctuations, as an eﬀect of b oundary ﬂuxes, and at the same time, the fastest conv ergence of the shap e parameter estimates to the theoretical v alue. Since the energy is transp orted mostly along the zonal direction by the zonal mean ﬂo w, by a veraging along a latitudinal band the highest frequencies are ﬁltered out, and ﬂuctuations with low er frequencies b ecome stronger. In the case of the zonally-a veraged observ ables, w e obtain w eak auto-correlations and fast conv ergence to ξ δ in the middle of the c hannel where the baro clin- icit y is the strongest, v ersus high auto-correlations and slow conv ergence at the margins of the c hannel where instead the baro clinicit y is w eak. The stronger the serial correlation, the less the n umber of uncorrelated data in a blo c k, and the larger blo ck sizes are needed in order to approac h asymptotic levels (see also [11]). By av eraging ov er the mid-latitude area, one merges zonally-av eraged time series exhibiting diﬀeren t auto-correlations. Thereb y , the con- v ergence to ξ δ is faster than in case of the zonally-a veraged observ ables at marginal latitudes. T o sum up, a v ery imp ortant conclusion of our study is the existence of latitude-dep enden t pre-asymptotic diﬀerences, as a coun terpart to the universal asymptotic prop erties. W e assume that the extremely slow con vergence has to do mainly with the fact that ξ δ is negativ e but v ery close to 0. Based on ξ δ and on the estimated GP mo diﬁed scale parameter, one is able to estimate according to Lucarini et al. [22] the absolute maxim um, whic h is the upp er end p oint of the GP distribution, as men tioned in Sec. 4.1. By p erforming a v ery rough estimation (and neglecting the weak latitude-dep endence of the GP mo diﬁed scale parameter), the absolute maximum in case of the lo cal observ ables A max ≈ 12 . 5, whic h is ab out 200 times the mean lo cal energy v alue (see Fig. 1) and 20 times larger than some of the largest estimated return levels obtained for the largest return times considered here (see Fig. 4). This means that extremes are bounded, and an absolute maxim um do es exist, but the tail is extremely stretc hed out, and ultra long simulations are needed to explore this absolute maxim um. Our results p oint out the discrepancy b et ween the existence of a mathematical limit and the actual p ossibilit y of observing it. Note that if the asymptotic shape parameter is lo wer, the absolute maxim um will be muc h closer to the maximum observ ed within a long y et ﬁnite time series, as it is shown in a recen t study on temp erature extremes in Southern P akistan [42]. Our conclusions regarding the conv ergence of the estimated shap e parameter to ξ δ are con- 25 ﬁrmed by results based on the GP mo diﬁed scale and return lev el estimates, in the case of the lo cal observ ables. W e p oin t out, how ev er, that the longer the return p erio d, the slow er the conv ergence of the estimated return lev els to their asymptotic v alues, and the larger the underestimation of the asymptotic return levels if w e consider small blo c k sizes (lo w thresholds). In the case of weak forcing, temporal and spatial correlations are very strong due to a regime b eha viour of our system, which exhibits t w o regimes: a more unstable one with stronger ﬂuctuations and a less unstable one with reduced ﬂuctuations. Due to this regime behaviour the statistics of extreme even ts are “contaminated”: if the blo ck size (threshold) is not large (high) enough, w e select ev ents from b oth regimes, whereas if it is large (high) enough, only extremes from the more unstable regime are selected. This induces non-monotonic c hanges of the estimated shap e parameters by increasing the blo c k size (threshold), and leads to the app earance of p ositiv e, i.e., non-physical, or v ery low shape parameter estimates. In the case of the lo cal observ ables, the estimated shap e parameters seem to conv erge at almost ev ery latitude to a v alue which is low er ( ≈ − 0 . 06) than the theoretical shap e parameter ( ξ δ = − 0 . 03). F urthermore, in the case of the zonally-av eraged and a verage mid-latitude observ ables, we cannot detect any conv ergence at all. The inconsistency of our numerical results with the theory is, in fact, unsurprising given the qualitativ e prop erties of the system in the low forcing regime, whic h do not resem ble characteristics of Axiom A systems (at least on the ﬁnite time scales we are able to explore based on the a v ailable data). Our results show that with increasing blo c k size or threshold the shap e parameters of the GEV and GP distributions are b ecoming more and more similar, according to the asymptotic equiv alence of the tw o mo dels [11, 22]. Both metho ds sho w us basically the same picture regarding the statistical properties of extreme ev ents. Despite the men tioned similarities, we observ e also some diﬀerences betw een the t w o approaches. The con vergence to the limit- ing distribution seems to b e somewhat faster in the case of the POT approach. This is in agreemen t with the well-established fact that the POT approach pro duces often more ac- curate predictions in case of applications [31, 11]. Despite the faster con vergence, how ev er, the b est GP shap e parameter estimates (deﬁned in Sec. 4.1) do not approximate ξ δ more accurately than the b est GEV shap e parameter estimates. Therefore, the adv an tage of the POT approach compared the BM approach is irrelev an t in the case of v ery long time series. W e use the Kolmogoro v-Smirnov test (KS-test) to v erify the ﬁt of the GEV (GP) distribution to the distribution of extremes, selected as blo c k maxima (threshold exceedances). Our results show that the KS-test is merely an indicator of the ﬁt quality , and do es not show whether the conv ergence to the correct GEV (GP) distribution is reac hed or not. The KS test suggests a go o d ﬁt to the GEV (GP) distribution ev en in cases when the distance b et w een the estimated and the asymptotic shape parameter is substantial and even if no conv ergence can b e detected. The misleading prop ert y of p-v alues of the KS and Pearson’s chi-squared tests w as also p oin ted out in previous studies in the case of more simple systems [8, 25]. In this w ork, w e estimate the GEV and GP parameters p erforming maxim um lik eliho od estimation [11], but it would b e relev an t to ﬁnd out to what exten t other estimation pro cedures, like the L-momen ts [43] or probabilit y-w eighted moments methods [44], would change the results. Concluding, we w ould lik e to emphasise some k ey messages one can get from our results: • Indeed, we hav e b een able to ﬁnd the signature of the univ ersal prop erties of the ex- tremes of physical observ ables in strongly c haotic dynamical systems, as predicted in 26 the case of Axiom A systems. Nonetheless, giv en the av ailability of v ery long y et ﬁnite time series, w e ha v e been able to ﬁnd more con vincing results (yet with a relativ ely large uncertain ty) only for sp eciﬁc observ ables, b ecause in the case of observ ables featuring serial correlations it is extremely hard to collect robust statistics of extremes. • W e hav e observ ed that in the case of strong forcing the estimate of the shape parameter increases monotonically tow ards its asymptotic v alue for stricter and stricter criteria of selection of extremes. This corresp onds to the fact that we manage to collect more detailed information on the local prop erties of the attractor near the p oin t of absolute maxim um of the observ able, and thus explore all the dimensions of the attractor. • W e also remark that agreemen t of the results with the theory of extremes of observ ables of dynamical systems developed in the con text of Axiom A ﬂows cannot b e found in the case of the weakly c haotic ﬂow featuring regime b eha viour and strong spatial and temp oral correlations, according to what one would exp ect considering the chaotic h yp othesis. • Finally , we note that the predicted and estimated shap e parameters are extremely small so that the statistics of extremes is virtually indistinguishable, up to ultra long return p eriods, from what w ould be predicted by a Gum bel distribution ( ξ = 0), whic h emerges as the statistical mo del of reference for ph ysical extremes in high dimensional c haotic systems, and suggests in the case of ﬂuids the existence of a well-dev elop ed turbulent state. Ho w ever, if the considered observ able liv es on a lo w er dimensional attractor, it is theoretically p ossible that the asymptotic shap e parameter is muc h low er, and th us distinguishable from 0. In case of applications, one deals with data sets having heterogeneous prop erties (non-smooth observ ables, existence of m ultiple time scales, eﬀects of spatial inhomogenities), and th us the estimated shape parameter can exhibit c hange of signs, ﬂuctuations, or (maybe even several) stable v alues only for a limited range of blo ck sizes or thresholds (some of these features can b e observed in our study in case of weak forcing, but can b e found also in lo wer dimensional models, as sho wn by Bo dai 2017, or ev en in more realistic mo dels with orography , as observed b y V annitsem 2007). In some cases, an extreme v alue la w does not ev en exist or is extremely diﬃcult to b e explored, and in other cases other distributions than the Gumbel, i.e. F rec het or W eibull, can b e more appropriate for the estimation of return levels up to a limited return p erio d. Ac kno wledgemen ts W e w ould like to thank to Sebastian Sc h ub ert, Christian F ranzk e, Maida Zahid and Richard Blender for useful discussions. W e are indebted to Sebastian Sc hubert for his support in p erforming the simulations and pro viding the co de for computing the Lyapuno v Exp onen ts and Kaplan-Y ork dimensions. V alerio Lucarini ackno wledges the many exc hanges on these topics with Davide F aranda, Antonio Sp eranza, and Renato Vitolo. V alerio Lucarini also ac knowledges supp ort received from the Sfb/T ransregio pro ject TRR181. V. M. G´ alﬁ ac- kno wledges funding from the International Max Planck Research Sc ho ol on Earth System Mo delling (IMPRS-ESM) The authors ackno wledge an anonymous referee for suggesting an 27 alternativ e deﬁnition of the accuracy . Conﬂicts of in terest The author(s) declare(s) that there is no conﬂict of in terest regarding the publication of this pap er. 28 References [1] IPCC. Climate Change 2013: The Physical Science Basis. Contribution of W orking Group I to the Fifth Assessmen t Report of the Intergo v ernmental P anel on Climate Change. Cambridge Universit y Press, Cam bridge, United Kingdom and New Y ork, NY, USA, 2013. [Sto c k er, T.F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S.K., Bosc hung, J., Nauels, A., Xia, Y., Bex, V., Midgley ,P .M. (eds.)]. [2] R. A. Fisher and L. H. C. Tipp ett. Limiting forms of the frequency distribution of the largest or smallest mem b er of a sample. Pro ceedings of the Cambridge Philosophical So ciet y, 24:180–190, 1928. [3] B. Gnedenk o. Sur la distribution limite du terme maxim um d’une serie aleatoire. Annals of Mathematics, 44:423–453, 1943. [4] J. Pick ands. Statistical inference using extreme order statistics. The Annals of Statistics, 3:119–131, 1975. [5] A. A. Balk ema and L. de Haan. Residual life time at great age. The Annals of Probability, 2:792–804, 1974. [6] M. Ghil, P . Yiou, S. Hallegatte, B. D. Malam ud, P . Nav eau, A. Soloviev, P . F riederichs, V. Keilis-Borok, D. Kondrasho v, V. Kossob ok ov, O. Mestre, C. Nicolis, H. W. Rust, P . Shebalin, M. V rac, A. Witt, and I. Zaliapin. Extreme ev en ts: dynamics, statistics and prediction. Nonlin. Pro cesses Geoph ys., 18:295–350, 2011. [7] V. Lucarini, D. F aranda, A. C. M. F reitas, J. M. F reitas, M. Holland, T. Kuna, M. Nicol, M. T o dd, and S V aienti. Extremes and Recurrence in Dynamical Systems. Wiley , 2016. [8] D. F aranda, V. Lucarini, G. T urchetti, and S. V aienti. Numerical conv ergence of the blo c k-maxima approach to the Generalized Extreme Value distribution. J. Stat. Phys., 145(5):1156–1180, 2011. [9] M. R. Leadb etter, G. Lindgren, and H. Ro otz ´ en. Extremes and Related Prop erties of Random Sequences and Pro cesses. Springer V erlag, 1983. [10] H. W. Rust. The eﬀect of long-range dep endence on mo delling extremes with the gen- eralised extreme v alue distribution. Eur. Phys. J. Sp ecial T opics, 174:91–97, 2009. [11] S. Coles. An Introduction to Statistical Mo deling of Extreme V alues. Springer V erlag, 2001. [12] F. J. Jr. Massey . The Kolmogorov-Smirno v test for goo dness of ﬁt. Journal of the American Statistical Asso ciation, 46:68–78, 1951. [13] T. W. Anderson and D. A. Darling. A test of go o dness of ﬁt. Journal of the American Statistical Asso ciation, 49:765–769, 1954. [14] A. Agresti. An In tro duction to Categorical Data Analysis. John Wiley and Sons, 2 edition, 2007. [15] J.-P . Ec kmann and D. Ruelle. Ergo dic theory of c haos and strange attractors. Reviews of Mo dern Physics, 57:617–656, 1985. [16] D. Ruelle. Chaotic Ev olution and Strange A ttractors. Cambridge Univ ersity Press, 1989. I [17] G. Galla v otti and E. Cohen. Dynamical ensembles in stationary states. Journal of Statistical Physics, 80:931–970, 1995. [18] G. Gallav otti. Nonequilibrium and irrev ersibilit y. Springer, 2014. [19] V. Lucarini, D. F aranda, and J. W outers. Univ ersal b eha viour of extreme v alue statistics for selected observ ables of dynamical systems. J. Stat. Ph ys., 147:63–73, 2012a. [20] V. Lucarini, D. F aranda, G. T urchetti, and S. V aienti. Extreme v alue theory for singular measures. Chaos, 22:023135, 2012b. [21] M. P . Holland, R. Vitolo, P . Rabassa, A. E. Sterk, and H. W. Bro er. Extreme v alue laws in dynamical systems under ph ysical observ ables. Physica D, 241(5):497–513, 2012. [22] V. Lucarini, D. F aranda, J. W outers, and T. Kuna. T ow ards a general theory of extremes for observ ables of c haotic dynamical systems. J. Stat. Phys., 154:723–750, 2014. [23] Luis Barreira, Y ako v Pesin, and Jorg Schmeling. Dimension and pro duct structure of h yp erb olic measures. Annals of Mathematics, 149(3):755–783, 1999. [24] E. Ott. Chaos in dynamical systems. Cam bridge Universit y Press, 1993. [25] T am´ as B´ odai. Extreme v alue analysis in dynamical systems: Tw o case studies. In C. F ranzke and T. O’Kane, editors, Nonlinear and Sto c hastic Climate Dynamics. Cam- bridge Universit y Press, Cam bridge, 2017. [26] E. Lorenz. Predictabilit y – a problem partly solved. In Seminar on Predictability, v olume I. ECMWF, 1996. [27] S. V annitsem. Statistical properties of the temp erature maxima in an in termediate order quasi-geostrophic mo del. T ellus, 59A:80–95, 2007. [28] M. F elici, V. Lucarini, A. Sp eranza, and R. Vitolo. Extreme v alue statistics of the total energy in an in termediate-complexit y mo del of the midlatitude atmospheric jet. Part I: Stationary case. Journal of Atmospheric Sciences, 64:2137–2158, 2007a. [29] S. Sch ubert and V. Lucarini. Co v arian t Ly apunov v ectors of a quasi-geostrophic baro- clinic mo del: analysis of instabilities and feedbacks. Quarterly Journal of the Ro yal Meteorological So ciet y, 141:3040–3055, 2015. [30] S. Sc hubert and V. Lucarini. Dynamical analysis of blo c king even ts: spatial and temp oral ﬂuctuations of cov ariant Ly apunov vectors. Quarterly Journal of the Ro y al Meteorological So ciet y, 142:2143–2158, 2016. [31] A. C. Da vison and R. L. Smith. Mo dels for exceedances ov er high thresholds. J. R. Statist. So c. B, 52:393–442, 1990. [32] M. R. Leadb etter. On extreme v alues in stationary sequences. Z. W ahrsc heinlichk eitstheorie verw. Geb., 28:289–303, 1974. [33] M. R. Leadb etter, I. W eissman, L. de Haan, and H. Ro otzen. On clustering of high v alues in statistically stationary series. In J. Sansom, editor, Proceedings of the 4th In ternational Meeting on Statistical Climatology, pages 217–222, W ellington, 1989. New Zealand Meteorological Service. I I [34] R. L. Smith. Extreme v alue analysis of en vironmental time series: An application to trend detection in ground-level ozone. Statistical Science, 4:367–393, 1989. [35] C. A. T. F erro and J. Segers. Inference for clusters of extreme v alues. J. R. Statist. So c. B, 65:545–556, 2003. [36] N. A. Phillips. The general circulation of the atmosphere: a n umerical exp erimen t. Q. J. R. Meteorol. So c, 82:123–164, 1956. [37] T. F risius. A mechanism for the barotropic equilibration of baro clinic wa ves. Journal of A tmospheric Sciences, 55:2918–2936, 1998. [38] J. R. Holton. An introduction to dynamic meteorology. Elsevier Academic Press, 2004. [39] G. Benettin, L. Galgani, A. Giorgilli, and J.-M Strelcyn. Lyapuno v characteristic expo- nen ts for smo oth dynamical systems and for Hamiltonian systems; a metho d for com- puting all of them. Part 1: Theory . Meccanica, 15:9–20, 1980. [40] P . F. Dunn. Measuremen t and data analysis for engineering and science. CR C Press, 2010. [41] C. F ranzke, A. J. Ma jda, and E. V anden-Eijnden. Low-order stochastic mo de reduction for a realistic barotropic mo del climate. Journal of the A tmospheric Sciences, 62:1722– 1745, 2005. [42] M. Zahid, R. Blender, V. Lucarini, and M. C. Bramati. Return lev els of temp erature extremes in Southern Pakistan. Earth Syst. Dynam. Discuss., in review, 2017. [43] R. M. Hosking. L-moments: Analysis and estimation of distributions using linear com- bination of order statistics. J. R. Statist. Soc. B, 52:105–124, 1990. [44] R. M. Hosking, J. R. W allis, and E. F. W o o d. Estimation of the Generalized Extreme- Value distribution b y the metho d of Probabilit y-Weighted Moments. T ec hnometrics, 27: 251–261, 1985. I II

Convergence of extreme value statistics in a two-layer quasi-geostrophic atmospheric model

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment