Complexity

Volume 2017, Article ID 5340858, 20 pages

https://doi.org/10.1155/2017/5340858

## Convergence of Extreme Value Statistics in a Two-Layer Quasi-Geostrophic Atmospheric Model

^{1}Meteorological Institute, CEN, University of Hamburg, Hamburg, Germany^{2}IMPRS-ESM, Max Planck Institute for Meteorology, Hamburg, Germany^{3}Department of Mathematics and Statistics, University of Reading, Reading, UK^{4}Centre for Environmental Policy, Imperial College London, London, UK

Correspondence should be addressed to Vera Melinda Gálfi; ed.gpm.temipm@iflag.adnilem-arev

Received 14 April 2017; Accepted 10 July 2017; Published 6 September 2017

Academic Editor: Gabriele Messori

Copyright © 2017 Vera Melinda Gálfi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We search for the signature of universal properties of extreme events, theoretically predicted for Axiom A flows, in a chaotic and high-dimensional dynamical system. We study the convergence of GEV (Generalized Extreme Value) and GP (Generalized Pareto) shape parameter estimates to the theoretical value, which is expressed in terms of the partial information dimensions of the attractor. We consider a two-layer quasi-geostrophic atmospheric model of the mid-latitudes, adopt two levels of forcing, and analyse the extremes of different types of physical observables (local energy, zonally averaged energy, and globally averaged energy). We find good agreement in the shape parameter estimates with the theory only in the case of more intense forcing, corresponding to a strong chaotic behaviour, for some observables (the local energy at every latitude). Due to the limited (though very large) data size and to the presence of serial correlations, it is difficult to obtain robust statistics of extremes in the case of the other observables. In the case of weak forcing, which leads to weaker chaotic conditions with regime behaviour, we find, unsurprisingly, worse agreement with the theory developed for Axiom A flows.

#### 1. Introduction and Motivation

The investigation of extreme events is extremely relevant for a range of disciplines in mathematical, natural, and social sciences and engineering. Understanding the large fluctuations of the system of interest is of great importance from a theoretical point of view, but also when it comes to assessing the risk associated with low probability and high impact events. In many cases, in order to gauge preparedness and resilience properly, one would like to be able to quantify the return times for events of different intensity and take suitable measures for preventing the expected impacts. Prominent examples are weather and climate extremes, which can have a huge impact on human society and natural ecosystems. The present uncertainty in the future projections of extremes makes their study even more urgent and crucial [1].

In practical terms, the main goal behind the study of extreme events is to understand the properties of the highest quantiles of the variable of interest. A fundamental drawback comes from the fact that extreme events are rare, so that it is difficult to collect satisfactory statistics from the analysis of a time series of finite length. Additionally, in the absence of a strong mathematical framework, it has its limits to make quantitative statements about the probability of occurrence of events larger than observed. Therefore, statistical inference based on empirical models tends to suffer from the lack of predictive power. A theoretical framework for analysing extreme events is provided by extreme value theory (EVT). After the early contributions by Fisher and Tippett [2], EVT was introduced by Gnedenko [3], who discovered that under rather general conditions the extreme events associated with stochastic variables can be described by a family of parametric distributions. As a result of this, the problem of studying the extremes of a stochastic variable can be reduced to estimating the parameters of a known probability distribution. The presumably most important parameter of such a distribution is called the shape parameter, which determines the qualitative properties of the distribution and, in particular, whether it has a finite or infinite upper point or, more concretely, whether extremes are bounded by an absolute finite maximum or not.

The two most important properties of extremes—their rare occurrence and their unusually high or low magnitude—constitute the basis of two popular methods of EVT, the block maxima (BM) and the peak-over-threshold (POT) approaches. The BM approach aims at finding the limiting distribution of maxima of independent identically distributed random (i.i.d.r.) variables separated into blocks of size , as . Under rather general conditions and for many different parent distributions, the limiting distributions of block maxima belong to the Generalized Extreme Value (GEV) distribution family; this is the point of view originally proposed by Gnedenko [3]. The POT approach aims at understanding the statistical properties of the exceedances of the random variable above a given threshold . Under the same conditions as for the BM approach, if is very high, the distribution of the threshold exceedances belongs to the Generalized Pareto (GP) distribution family [4, 5]. The existence of well-defined functional forms for the distributions describing extreme events provides predictive power: one can in principle compute the return time of events of yet unobserved magnitude. It is remarkable that these two points of view on extremes (which lead to different selection procedures and different choices of the events classified as extremes) are, in fact, equivalent. For a given random variable , the GEV and GP distributions describing its extreme events have the same asymptotic behaviour. Additionally, it is possible to express the parameters characterizing one distribution as a function of the parameters characterizing the other distribution; in particular, the shape parameter is the same for both.

Problems in applying EVT to actual time series result from the fact that, typically, the observed data feature a certain degree of serial correlations [6]. Note that the setting behind the construction of EVT can be extended by relaxing the hypothesis of independence of the random variables, as long as correlations are decaying sufficiently fast. This is of clear relevance when trying to use EVT for studying observables of deterministic dynamical systems. In this case, in fact, the underlying dynamics determines the existence of correlations between the values of observables at different times, and one can easily guess that, when the dynamical system is chaotic, there is good hope of deriving EVT for its observables [7]. Obtaining the true limiting EVT can be extremely hard, even in simple dynamical systems [8], or, in that case, it can be argued that they do not even exist [9]. When analysing finite time series, the convergence of the estimated GEV or GP shape parameters to near-zero asymptotic values can be very slow. The speed of convergence depends on the type of parent distribution [10] and can be additionally slowed down by correlations [11, 12]. Due to the fact that the dataset size is always limited, there is typically a difference between the asymptotic GEV or GP parameters and the estimated ones; finite-size estimates are generally biased. For example, the GEV shape parameter of a simple Gaussian process is 0, but, for any finite time series, we would estimate typically a negative shape parameter [2].

When performing statistical inference using the BM or POT method (fitting the GEV or GP model, resp., to data), it is crucial to have an appropriate protocol of selection of “good" candidates for extremes [12]. On the one hand, if the chosen blocks (for the BM method) are too short or the threshold (for the POT method) is too low, the approximation of the limit model is likely to be inadequate. Hence, the verification of the agreement between the statistical model and the available data is essential, which can be done based on goodness-of-fit tests, like the Kolmogorov-Smirnov [13], Anderson-Darling [14], or Pearson’s chi-squared tests [15]. On the other hand, if the blocks are too large or the threshold is too high, the number of extremes may be insufficient for a reliable estimation of the parameters, and the uncertainty becomes very high. As discussed later in the paper, Coles [12] shows how to derive an optimal choice for the value of the block size or the threshold, in such a way as to verify that we are close to the asymptotic level as required by EVT but we use the available data as efficiently as possible. In this regard, we introduce a new way of choosing an optimal block size or threshold based on a concept of accuracy of estimation, which can be actually carried out only when the true shape parameter is known from theory.

Classical EVT has been extended and adapted to analyse extremes of observables of chaotic dynamical systems, where the sensitive dependence on initial conditions is fundamentally responsible for generating a de facto stochastic process. The reader is referred to Lucarini et al. [7] for a detailed overview of the field of EVT for dynamical systems. It is possible to establish an EVT for observables of chaotic dynamical systems when one considers Axiom A dynamical systems. These are rather special chaotic dynamical systems that are uniformly hyperbolic on their attractor (so that stable and unstable directions are well separated), which supports a so-called Sinai-Ruelle-Bowen (SRB) measure. Such an invariant measure has physical relevance because it is stable against weak stochastic perturbations [16, 17]. One of the great merits of Axiom A dynamical systems is that they allow for deriving rigorous and robust statistical mechanical properties for purely deterministic background dynamics. Despite having deterministic dynamics, when looking at their observables, they behave just like generators of stochastic processes. While Axiom A systems are rather special and indeed not generic, they have great relevance for applications if one takes into account the chaotic hypothesis, which indicates that high-dimensional chaotic systems behave at all practical purposes as if they were Axiom A [18, 19].

Several studies dealing with EVT for dynamical systems reveal a link between the statistical properties of the extremes and geometric (and possibly in turn global dynamical) characteristics of the system producing these extremes [8, 20–24]. The main findings are that when suitable observables are chosen for the dynamical system of interest, it is possible to relate the GEV or GP parameters describing the extremes to basic properties of the dynamics and especially to the geometry of the measure supported by the attractor. In particular, depending on the choice of the observable, one can associate the most important parameter of the GEV or GP distribution, the shape parameter, with the information dimension of the measure supported on the attractor or with the partial information dimension along the stable and unstable directions of the flow [23]. These partial dimensions are well defined everywhere on the chaotic attractor, possibly with a variation with location, and also for nonuniformly hyperbolic systems [25], beside Axiom A systems. However, Axiom A systems possess an ergodic SRB measure which lends itself to a universality of the shape parameter for all sufficiently smooth observables; the local or pointwise (partial) dimensions take the same value almost everywhere [26]. In this case, the uniform shape parameter can be related to the (partial) Kaplan-Yorke dimension(s), which is (are) defined by the global dynamical characteristic numbers, the Lyapunov exponents. Clearly, this is an asymptotic result, and one must expect that differences emerge on preasymptotic levels when different observables are studied. In other words, this theory does not make predictions regarding the convergence of shape parameter estimates, the analysis of which is the main objective of this paper. Via the connection with fractal dimensions, it can be said that the analysis of extremes acts as a microscope able to assess the fine scale properties of the invariant measures.

Some preliminary numerical tests show that there is no obvious convergence to the predicted asymptotic shape parameter in low-dimensional cases [23]. Bódai [9] examined the convergence to the GEV distribution in the case of extremes of site variables in the Lorenz 96 model [27], investigating separately a range of cases extending from weak to strong chaos. He found that when considering configurations supporting weak chaos with a low-dimensional attractor, the theoretical results obtained in the context of the Axiom A hypothesis are hard to verify. For lower dimensions, up to a dimension of about 5, shape parameter estimates fluctuate greatly rather than converge, while block maxima data can be shown not to conform to a GEV model, and for somewhat larger dimensions, up to 9 in the study, estimates could diverge from the predicted value while data already conform to a GEV model. Good agreement with the theory was found only in the highly turbulent case possessing a higher-dimensional attractor, about 30, supporting the basic idea behind the chaotic hypothesis. Also in this case, nonetheless, very slow convergence was found.

In a previous analysis performed on higher-dimensional, intermediate complexity models with degrees of freedom, very slow (if any) convergence to EVT distributions could be found in the case of extremes of local temperature observables [28]. In another analysis of a similar model [29], the agreement of the distribution of global energy extremes with a member of the GEV family was indeed good, yet large uncertainty remained on the value of the shape parameter, and no stringent test was made to make sure that the estimate was stable against changes in the block size considered in the BM analysis. Clearly, the specific choice of the observable and the degree of chaoticity of the underlying dynamics is of primary relevance regarding the convergence to the limiting GEV or GP distribution.

In this work, we use a quasi-geostrophic (QG) atmospheric model of intermediate complexity, featuring 1056 degrees of freedom, to analyse extremes of different types of observables: local energy (defined at each grid point), zonally averaged energy, and the average value of energy over the mid-latitudes. Our main objective is to compare the estimated GEV and GP shape parameters with a shape parameter derived, based on the theory referred to above, from the properties of the attractor and of the measure supported on it along the stable, unstable, and neutral directions. We refer to this as the “theoretical shape parameter.” Thus, we explore numerically the link between the purely statistical properties of extreme events based on EVT and the dynamical properties of the system producing these extremes. We perform simulations applying two different levels of forcing: a strong forcing, producing a highly chaotic behaviour of the system, and a weak forcing, producing a less pronounced chaotic behaviour. The dimensionality of the attractor is much larger in the former than in the latter case. This work goes beyond the previously mentioned studies, based on more simple dynamical systems, in a sense that with our model we can study the convergence for observables being different physical quantities or representing different spatial scales/characteristics of the same physical quantity. Additionally, compared to previous studies also performed on intermediate complexity models, we consider longer time series and a variety of observables. Our model is simple compared to a GCM (General Circulation Model) but contains two of the main processes relevant for mid-latitude atmospheric dynamics: baroclinic and barotropic instabilities. Hence, we contribute to bridging the gap between the analysis of extremes in simple and very high-dimensional dynamical systems, as in the case of the GCMs used for atmospheric and climate simulations, by using a model that simulates to a certain degree Earth-like atmospheric processes and allows also for computing with feasible computational costs some dynamical system properties, like Lyapunov exponents or Kaplan-Yorke dimensions. The properties of the model have been extensively studied by Schubert and Lucarini [30, 31].

Based on numerical results (Sebastian Schubert, personal communication), the model is expected to be nonhyperbolic, but one can assume that the chaotic hypothesis applies to it (in the case of sufficiently high forcing levels inducing a strongly chaotic behaviour of the system), and so in analysing the convergence of shape parameter estimates one can take the predicted theoretical shape parameter as reference. We also assume that the symmetry of the model with respect to longitude, which introduces a central direction besides the directions of expansion and contraction in phase space, does not alter the ergodicity of the system at a practical level, and hence the true shape parameter remains uniform.

Although we use an idealised model, our results are transferable to time series obtained from more realistic model simulations or from measurements. By understanding the differences among the analysed observables, we gain insight into the statistical properties of extremes of geophysical observables with different spatial scales. By using two forcing levels, we are able to study the convergence to theoretical shape parameters related to different chaotic systems: one exhibiting fast decaying correlations and another one characterised by slower decaying correlations. These aspects are relevant in the case of geophysical applications where one deals also with time series on several spatial scales and with different degrees of correlations.

The structure of this article is as follows. Section 2 gives a theoretical overview, describing the block maxima and the peak-over-threshold approaches. In Section 3, we present our model, the performed simulations, and the applied methods. In Section 4, we discuss our results regarding the statistics of extremes for strong forcing and for weak forcing. We summarise and discuss our results in Section 5.

#### 2. Elements of Extreme Value Theory

Let us consider , where is a sequence of i.i.d.r. variables with common distribution function . The extremal types theorem [2, 3] states that if there exist sequences of constants and , so that the distribution of normalised , that is, , converges for to a nondegenerate distribution function , then is one of three possible types of the so-called extreme value distributions, having the cumulative distribution function where , , for , and for . Note that the form of on the first line of (1) is the generic form, valid for [12], as its limit as is the form written out on the second line. Similar comments apply to (2), (3), and (4).

represents the GEV family of distributions with three parameters: the location parameter , the scale parameter , and the shape parameter . The shape parameter describes the tail behaviour and determines to which one of the three types of extreme value distributions belongs. If , the tail decays exponentially, and is a type I extreme value distribution or Gumbel distribution. If , the tail decays polynomially, and belongs to type II or Fréchet distribution. If , the domain of the distribution has an upper limit and is referred to as a type III or Weibull distribution.

Under the same conditions, for which the distribution of converges to the GEV distribution, the exceedances above a threshold reaching the upper right point of the distributions of , given that , are asymptotically distributed according to the Generalized Pareto (GP) distribution family [12]where for , , and . has two parameters: the scale parameter and the shape parameter . The shape parameter describes again the tail behaviour and determines to which one of the three types of GP distributions belongs. If , the tail of the distribution decays exponentially; if , the tail decays polynomially; and if , the distribution is bounded [4, 5, 32]. If convergence to the GEV and GP distributions is realised, and . As a result, once we estimate the parameters for the GEV, we can derive the corresponding GP parameters (including the threshold ) and vice versa [12].

From the values of the GEV or GP parameters, it is possible to infer the expected return levels or extreme quantiles. Return levels are obtained from the GEV distribution by inverting (1): where and represents the return period. In the case of the GP distribution, the -observation return level (i.e., the level that is exceeded on average every observations) can be derived from (2): where represents the probability of an individual observation exceeding the threshold [12]. By plotting the GEV (GP) return level () against the return period () on a logarithmic scale, the plot is linear if (), convex if (), and concave if ().

In the case of a correlated stationary stochastic process, the same GEV limit laws apply as for i.i.d.r. variables if certain conditions, regarding the decay of serial correlation, are fulfilled [7, 33, 34]. By stationary, we refer to a sequence of correlated variables whose joint probability distribution is time-invariant. However, an important restriction is that, as an effect of serial correlation, we have to introduce an effective size for the block, which is smaller than the actual number of observations contained in it. This can enhance the bias in the parameter estimation, appearing as a slower or delayed convergence of the block maxima distribution to the limiting GEV distribution [11, 12]. Another possible effect of serial correlation is the appearance of extremes at consecutive time steps (clusters). If an extreme value law does exist in this case, then , where is called the extremal index and ( denotes the limiting distribution of BM from the correlated sequence and the one from an uncorrelated sequence, having the same marginal distribution). Clusters of extremes represent a problem especially when applying the POT approach. A widely adopted method to get rid of the correlated extremes is declustering, which basically consists of identifying the maximum excess within each cluster and fitting the GP distribution to the cluster maxima [34–36].

As mentioned before, several studies on EVT for observables of dynamical systems relate the GEV and GP shape parameters to certain properties of the system itself. In the case of the so-called “distance” observables, one can relate the GEV and GP parameters to basic geometrical properties of the attractor [8, 20, 21, 24]. The distance observables are functions of the Euclidean distance between one point on the attractor and the orbit . The function is chosen in a way to have a global maximum for , so that large values of correspond to recurrences of the orbit near . Depending on the choice of the function , the extremes of the distance observables can have positive, negative, or vanishing value for the shape parameter. In particular, when is chosen to have the form for , the shape parameter is negative, and it is proportional to the inverse of the Kaplan-Yorke or information dimension of the SRB measure supported by the attractor; see below for further details [20, 21, 23].

While recurrence properties are indeed important for characterising a system, distance observables are not well suited for studying some basic physical properties, such as, in the case of fluids, energy or enstrophy. Hence, Holland et al. [22] studied the extremes of smooth functions which take their maximum on the attractor in a point where the corresponding level surface of is tangential to the unstable manifold of the attractor, referring to them as “physical” observables. They found a relationship between the GEV shape parameter and some geometrical properties of the measure supported by the attractor dealing with the properties of the unstable and stable directions in the tangent space. The results of Holland et al. [22] were reexamined by Lucarini et al. [23], using the POT approach for physical observables of Axiom A systems. They considered the time-continuous time series of physical observables and found that for all nonpathological physical observables the shape parameter can be written as with defined as where , , and are the partial information dimensions of the SRB measure supported by the attractor restricted to the stable, unstable, and neutral (i.e., central) directions. As mentioned in Section 1, these local or pointwise (partial) dimensions take the same value almost everywhere on the attractor if one considers smooth observables of systems possessing an ergodic measure. is equal to the number of positive Lyapunov exponents [26], is equal to the number of zero Lyapunov exponents, and , where is the Kaplan-Yorke dimension with denoting the Lyapunov exponents of the system, in a descending order, and is such that is positive and is negative. We remark that a more general point of view, taking into consideration possible geometrical degeneracies, suggests that , and, additionally, [7].

According to (5) the shape parameter is always negative (due to the compactness of the attractor), and it is close to zero in the case of systems having large values of the Kaplan-Yorke dimension. Furthermore, it shows a universal property of extremes, which does not depend on the chosen observable but only on the geometry of the measure. In what follows, we will focus on comparing (5) with statistically inferred GEV and GP shape parameters in the case of energy extremes of the model investigated here, which is described next.

#### 3. Model Description and Methods

We consider a spectral quasi-geostrophic (QG) 2-layer atmospheric model similar to the one introduced by Phillips [37]. Specifically, our model, including the simulation code, is the same as in [30] and is a modified version of the one presented in [38]. The model represents synoptic scale mid-latitude atmospheric dynamics based on the quasi-geostrophic approximation, which assumes hydrostatic balance and allows only small departures from the geostrophic balance. The model features baroclinic conversion and barotropic stabilisation processes and simulates a turbulent jet-like zonal flow when suitable values are chosen for the parameters of the system. The reader is referred to [39] for a detailed physical and mathematical description of quasi-geostrophic approximation for mid-latitude atmospheric dynamics.

##### 3.1. Model Description

The model domain is a rectangular channel with latitudinal and longitudinal coordinates . represents the Equator and corresponds to the north Pole. We assume periodicity along the -direction, so that corresponds to the length of the parallel at 45°N. The vertical structure of the model atmosphere consists of only two discrete layers: this is the minimal vertical resolution needed to represent baroclinic processes [39]. Five vertical pressure levels define the two layers with boundaries at hPa (surface level), hPa, hPa, hPa, and hPa (top level). The geostrophic stream function is defined at levels and , and , where the quasi-geostrophic vorticity equation for the mid-latitude -plane (7)-(8) is applied, while the vertical velocity is specified at level , where the thermodynamic energy equation (9) is valid.

The model is described by the following equations in terms of the barotropic stream function , baroclinic stream function , and temperature : In the above, we expressed the advection in terms of the Jacobian operator defined as . represents the static stability parameter [39]. We define the stability parameter , where is the Rossby radius of deformation. The name and values of model parameters are listed in Table 1.