The Use of Cluster Analysis in Climate Science

Cluster analysis is used in the analysis of weather and climate to identify distinct, discrete groups of atmospheric circulation states, and sequences of states (also called trajectories), from observations, reanalyses, and numerical model simulations and predictions (Ghil & Robertson, 2002; Hannachi, Straus, Franzke, Corti, & Woollings, 2017; Straus, Molteni, & Corti, 2017). (Reanalyses are sequences of full three-dimensional atmospheric and oceanic states constructed by combining observations with short-range model forecasts.)

One of the goals of cluster analysis is to provide physical insight into the states that are preferred and also possibly unusually persistent (Kimoto & Ghil, 1993a), when such states can be identified and distinguished from the continuous background distribution of geophysical variables. Preferred states (and their sequences) are those that are significantly more likely to occur than would be predicted by a suitable background distribution, such as a multivariate Gaussian distribution (Toth, 1993), while persistent states are those with lifetimes distinctly longer than those of the background states (Dole & Gordon, 1983; Vautard & Legras, 1988; Toth, 1992).

Preferred states arise due to mechanisms that operate on multiple timescales, as reviewed by Hannachi et al. (2017). On seasonal and inter-annual timescales these states are related to the interaction of Rossby waves forced by topography (mountain heights) and the distribution of land and ocean surface temperatures, and in midlatitudes to the remote forcing by large anomalies in seasonal mean tropical diabatic heating, such as those due to the El Niño Southern Oscillation. On intra-seasonal timescales, preferred states arise due to the interaction of rapidly propagating transient eddies associated with midlatitude storms with regional jet streams, and with Rossby waves forced by transient tropical diabetic heating.

A separate goal of cluster analysis is to provide a physical classification of weather and climate patterns for several purposes. One is to provide a good classification of the basic surface climates on earth (e.g., Netzel & Stepinski, 2016). Another is to better understand and forecast the occurrence of weather extremes (storms, floods), climate extremes (heat waves, cold snaps), and other factors relating to human health (e.g., Cassou, Terray, & Phillips, 2005; Coleman & Rogers, 2007; Polo, Ullmann, Roucou, & Fontaine, 2011). In these applications, the emphasis is on reducing the very large nominal dimensionality of the atmospheric dynamic and thermodynamic space to a manageable, discrete set of weather and climate states that have predictive value in terms of how the climate system affects human beings. On weather timescales, these states often correspond to patterns familiar to weather forecasters (James, 2006; Ferranti & Corti, 2011; Ferranti, Corti, & Janousek, 2015).

The Goals of This Article

This article is meant to be an in-depth introduction to the many techniques used in clustering techniques in the field of climate and weather analysis. The techniques most commonly used are described in enough detail so that the reader can better understand the literature, while the interested user can make an intelligent decision as to the best technique for a given problem. The resources necessary to further explore a given technique are provided in the bibliography. In cases where the basic algorithms can be described succinctly, and shed some light on the clustering technique, the algorithms are described. In other cases, the algorithms and their underlying statistical theory are too complex to lie within the scope of this article, but appropriate references are given.

The article does not attempt to follow a historical chronology, as this would not have been very useful in terms of understanding the clustering techniques themselves or their use in climate science.

Space-Time Objects of Clustering

In climate science applications, the various states on which clustering are carried out are usually geophysical fields, or sequences of these fields. Examples of fields are geopotential height, temperature, and horizontal winds on a suitable surface. These fields may be global in extent, but are most often hemispheric or regional. In some special applications, the states are not fields, but the locations of well-defined features, such as midlatitude or tropical cyclones, and it is their trajectories that are subject to clustering. The clustering of trajectories is discussed further in the section “Dynamic Clustering Techniques.”

The states (fields) may arise from time averages (such as annual, seasonal, monthly, or even five-day averages, also called pentads), or they may represent instantaneous (“snapshot”) states taken once per day or more often. In the majority of applications, the annual cycle of the particular field is assumed to be known, and it is the departures from this annual cycle (“anomalies”) that are of interest. However, this is not always the case. When the annual cycle interacts very strongly with the phenomena of interest, clustering may be carried out on the full fields (see, e.g., Qian, Robertson, & Moron, 2010; Moron, Robertson, Qian, & Ghil, 2014).

In addition to removing the annual cycle, further filtering of the states to be clustered is also appropriate for many applications. In order to identify preferred states of low frequency variability on intra-seasonal timescales (roughly 10- to 90-day periods), a “low-pass” time filter is often applied, as in Vautard (1990). Such a filter typically removes fluctuations with timescales less than roughly 6 to 10 days, fluctuations that are characteristic of rapidly moving baroclinic systems in midlatitudes.

In some applications the states may not be time filtered, but rather prechosen to correspond to special situations. These may include persistent periods (Toth, 1992), periods of extreme precipitation (Zhao, Deng, & Black, 2017a), extreme dryness (Zhao, Deng, & Black, 2017b), extreme temperature (Xie, Black, & Deng, 2017), or periods of particular synoptic developments such as surges of cold air (Park, Ho, & Deng, 2014).

The two most common mathematical choices for representing states are in terms of grid point values, or in terms of principal components (PCs) and empirical orthogonal functions. The transformation to the PC representation is defined by the covariance matrix of the original field, whose eigenvectors yield a set of fixed patterns, each associated with its own time series (the PC). The PCs are obtained by projecting the original data onto the fixed patterns. More than one type of field can be combined into the definition of a single state; if the types have different units then each is divided by its own temporal standard deviation prior to the computation of the covariance matrix.

The individual state can be represented by an M-dimensional vector $z$, where M is either the number of grid points, or the number of PCs retained. The major advantage of using the PC representation is that M can be greatly reduced while still retain most of the space-time variability of the set of fields. Because the great majority of studies use the PC representation as a starting point, it will be assumed throughout this article unless otherwise stated.

Definitions of Distance (Dissimilarity)

A critical quantity in most clustering methods is the distance between two states (or sequences of states), also called the dissimilarity. We denote the N states in a data set as the set $\left\{z_{\alpha }\right\}(1\le \alpha \le N)$, where the PC coordinates of each M-dimensional vector (denoted by i or j) are Cartesian. The scalar distance $d_{\alpha \beta }$ between two states $z_{\alpha }$ and $z_{\beta }$ can be defined in a number of ways.

Euclidean Distance: This is the M-dimensional version of the Euclidean geometric definition of length:

This is the most often used measure of distance. When applied to meteorological maps, the quantity $d_{\alpha \beta }$ is called the root-mean-square difference. This distance measure is very often used in conjunction with hierarchical and partitioning methods of cluster analysis (e.g., Cheng & Wallace, 1993; Michelangeli, Vautard, & Legras, 1995).

Cosine Distance: This distance can be defined as:

where the M-dimensional dot product is $z_{\alpha }\cdot z_{\beta }={\displaystyle {\sum }_i{z}_{\alpha i}{z}_{\beta i}}$ and ${\alpha }_{\alpha \beta }$ is the angle between the two vectors. If all states are projected onto the M-dimensional unit hypersphere, this is a natural measure of the distance on that sphere. Such a projection is used, for example, by Crommelin (2004) and Mo and Ghil (1988). Note that if two physical fields are represented using principal components and empirical orthogonal functions, the (un-centered) spatial correlation is just the $\cos {\alpha }_{\alpha \beta }$.

Mahalanobis Distance: This is a measure of the dissimilarity between two vectors $z_1$ and $z_2$, which are drawn from a distribution D of states $\left\{z_{\alpha }\right\}(1\le \alpha \le N)$ characterized by the mean $\mu$ and the covariance matrix $S_{ij}$:

The Mahalanobis distance is given by:

Several properties of this distance measure are worth noting: it is preserved under invertible linear coordinate transformations. If the mean vanishes, and the covariance matrix is diagonal, one can scale the coordinates so that this measure of distance becomes the Euclidean measure. The Mahalanobis distance appears in the definition of a multivariate Gaussian probability distribution, and in this context is used in the mixture modeling approach to probabilty density estimation (Stephenson, 1997; Stephenson, Hannachi, & O’Neill, 2001).

Dynamic Time Warping Dissimilarity: This measures the dissimilarity between two sequences of states $\left\{z_{\alpha }\right\},1\le \alpha \le N$ and $\left\{y_{\beta }\right\},1\le \beta \le N^{\prime }$ which are not necessarily of equal length. This method does not assume a one-to-one association between each member of the two sequences, but rather considers the $N\times N^{\prime }$ matrix whose elements $d_{\alpha \beta }=D(z_{\alpha },y_{\beta })$ are given by the dissimilarity between the two states $z_{\alpha }$ and $y_{\beta }$ (using any measure of dissimilarity). If each element of this matrix is considered a small box, a one-dimensional path through this matrix would consist of choosing a set of connecting boxes whose beginning is at one corner, and whose end is at the opposite corner. More formally, a warping path is defined as a contiguous set of matrix elements connecting $d_{11}$ to $d_{N{N}^{\prime }}$. The optimal warping path is then defined as that path that minimizes the sum of $d_{\alpha \beta }$ along the path. This sum is the dynamic time warping dissimilarity. Its use in climate science is for measuring the alignment between two time series that may be locally out of phase. As an example, Netzel and Stepinski (2016), who are concerned with defining global climate types, use this technique to compare the similarity between the annual cycles of temperature and precipitation in different regions.

Kernel Density Estimation

In one dimension (say that of the leading PC), the raw data can be represented by a simple histogram, in which, however, the bin width is somewhat arbitrary. In two or more dimensions $M\ge 2$, a histogram becomes unwieldy, and a continuous estimate is helpful. The concept behind the kernel density estimate is that the contribution of each data point (essentially a delta-function in the phase space) gets smoothed out to a local pdf (with unit area under the curve). The final (normalized) estimate of the pdf is then given by:

where the vector z is a general point in the PC coordinate space, the vectors $z_{\alpha }$ are given in the PC coordinates, the index $\alpha$ labels the coordinate, and h is a smoothing parameter (Silverman, 1986; Corti et al., 1999). The kernel function K can be a range of functions. One commonly used kernel function is the standard normal density function:

In fact the equation for $\widehat{f}(z)$ is over-simplified, because in more than one dimension a separate value of h in each dimension is appropriate; usually h is proportional to the standard deviation of the data in that dimension. Furthermore, it is very helpful to choose the smoothing in an adaptive manner, so that the parameters vary according to the value of the local pdf at each sample point in order to achieve more reliable estimates in data-sparse regions. See Kimoto and Ghil (1993a) for details on how this is done.

For a sample drawn from an underlying Gaussian (or normal) distribution one would expect the pdf estimate to also be normal, with a maximum at the point representing the overall mean (usually the origin). Significant deviations from this suggest preferred (and also unpreferred) regions of phase space. It is in this context that the choice of the smoothing parameters is of importance: if too large, local maxima in $\widehat{f}$ that represent significant preferred anomalies will be smoothed out, if too small, spurious local maxima due to sampling error may appear. An illustration of how increasing the smoothing parameters changes the appearance of the pdf in two-dimensions is given in Figure 1.

Open in new tab

Figure 1. Kernel density estimation of the probability distribution function in the space of the leading two principal components of filtered 500 hPa geopotential height for winter (November 16–March 31 for 1980–1981 through 2011–2012). The filter removed periods of 20 days or fewer. An adaptive algorithm was used. The kernel width is given by h times the value of the standard deviation, with the value of h indicated in each panel.

One simple approach to determining the set of parameters h is to use synthetic (or surrogate) data sets to define the background pdf. The purpose in using these synthetic data sets is to get some estimate of how likely multiple peaks in the pdf are due only to sampling error. Each of these synthetic data sets has the same size as the true data set, but the states in PC-space are generated by stochastic processes. We give a particular example here in some detail. In each realization of the synthetic data set, each PC is generated by an independent stochastic Gaussian process, whose variance and auto-correlation are similar to that of the corresponding original PC. The underlying probability distribution P of each PC is then Gaussian, hence has a single maximum. Considering for example two dimensions, the underlying joint probability distribution of variables x (value of PC-1) and y (value of PC-2) would be simply $P(x,y)=P_1(x)P_2(y)$, which would have one peak in the x/y plane. While the true $P(x,y)$ has only one peak, estimates of $P(x,y)$ obtained from finite samples may have multiple peaks.

Computing the sample pdf of many such synthetic data sets with a given set of h, one finds that p% of the samples show distinct maxima (not at the origin), maxima that must be due to sampling error. The smoothing parameters must be increased until p is reduced to an acceptable level, usually 5% or 1%. For other choices of the kernel function K, and alternative methods for choosing the values of h, see Silverman (1986).

Mixture Modeling

Many of the clustering methods described here have the property that individual states are either assigned to a single cluster uniquely, or may lie within a preferred region of the pdf. The mixture modeling method expands this framework to include the possibility that an individual state may be associated with more than one cluster, but with different probabilities. The approach is to model the pdf of the data set as a linear combination of k multivariate normal (Gaussian) distributions, or components, where each of the k components is proportional to:

where M is the dimension of the state vectors z (number of PCs), $\mu$ is the mean for that component, and $\Sigma$ the (symmetric) covariance matrix. The vertical bars denote the determinant, and the superscript T denoted vector transpose. Note the use of the Mahalanobis distance. We will denote this distribution for the jth component as:

emphasizing that this component depends on the actual state vectors $\left\{z_{\alpha }\right\}$ on the one hand, and the unknown parameters ${\mu }_j,{\Sigma }_j$ on the other. The final pdf can then by written as a sum of these component distributions, with weights that add to unity:

where ${\displaystyle {\sum }_{j=1}^k{w}_j=1}$. The quantity p_j (z) is used to indicate the contribution of component j to the overall probability at point z.

Each of the component distributions $g_j$ can then be thought of as a generalized cluster. The weight $w_j$ can be thought of as the probability that a particular state $z$ is in the cluster j. The parameters (mean ${\mu }_k$ and covariance ${\Sigma }_k$ and the weight w_j of each component) are determined as those that maximize the likelihood (or logarithm of the likelihood) given the data. In other words, we treat the observed state vectors as fixed, and take the partial derivatives of $p(z)$ with respect to each parameter and set that to zero. This yields a set of coupled nonlinear equations, not usually amenable to a direct solution.

The expectation maximization (EM) procedure can be adapted to this problem. This iterative procedure starts out with an initial guess for the parameters, and then estimates the probability that each data point $z_{\alpha }$ belongs to the jth component (cluster) as:

where the $p_j$ are evaluated with initial guesses for the parameters. Then at the next iteration we estimate:

The procedure is repeated until convergence is achieved.

Detailed explanations of this method are given in the context of its application to climate by Haines and Hannachi (1995), Hannachi (1997), Smyth, Ide, and Ghil (1999), Hannachi and O’Neill (2001), and Hannachi (2007). The optimum number of component distributions (clusters) k have been determined by a variety of methods: cross-validation using independent subsamples of the data (as in Branstator & Selten, 2009), cross validation as discussed by Smyth et al. (1999), or the use of order statistics (Hannachi, 2007). Note that the parameters to be chosen consist of the M coordinates for the mean of each component, and the M(M+1)/2 parameters for the covariance matrix, thus $kM(M+3)/2-1$ parameters in all. This large number of parameters tends to keep the number k of mixture components relatively low in order to avoid over-fitting of the sample data.

Hierarchical Method (Agglomerative)

The hierarchical method sequentially merges individual clusters (groups) of states pairwise, leading to fewer clusters as the process proceeds. At each step of this iterative method, a number of clusters $k$ have been defined, each containing a number of states. Each cluster is characterized by its own internal variance $E_k$ or error sum of squares, which is simply the sum of the squared distance of each state about the centroid of that cluster, the centroid being defined as that point whose coordinates are the average of the coordinates of the points in that cluster. Each step of the clustering is thus characterized by its total error sum of squares $E={\sum }_k{E}_k$. To proceed to the next step, two of the $k$ clusters are merged, thereby increasing $E$. The two clusters to be merged are chosen so as to minimize the increase in$E$.

To initiate the iteration, each state is considered its own cluster. The final step consists of merging two clusters into a single cluster containing all the states. The entire procedure can be represented by a tree (also called a dendogram), with the various branches merging as you move upward toward the top. The height at which two branches merge, thus defining a new cluster, can be chosen to be the increase in $E$ that characterizes the merging. The entire procedure is known as agglomerative hierarchical clustering. The final number of clusters can be chosen in several ways, including identifying that step that causes a kink, or sudden upturn in $E$, and determining which clusters which are most consistent in randomly chosen half-length data sets. A straightforward algorithm is described by Cheng and Wallace (1993), who apply this method to find clusters of the geopotential height field. Applications of hierarchical clustering to precipitation and temperature extremes can be found in Park et al. (2014), Xie et al. (2017), and Zhao et al. (2017a, 2017b).

Partitioning Methods

Clustering Trajectories

An adaption of the mixture modeling approach allows for the clustering of trajectories of tropical or midlatitude cyclones (Gaffney, Robertson, Smyth, Camargo, & Ghil, 2007; Camargo, Robertson, Gaffney, & Smyth, 2007a, 2007b). Here the state vectors $z_i$ now refer to the successive positions of the cyclone center in discrete time. To be concrete, let the n discrete times t be described by the indices $t_j$:

and let the longitudes ${\lambda }_j$ and latitudes ${\varphi }_j$ of a single cyclone be described by the following second order polynomial regression $(p=2)$ model in time:

where ${ϵ}_1$ and ${ϵ}_2$ represent Gaussian noise terms with zero mean and variances ${\sigma }_1^2$ and ${\sigma }_2^2$. The state vector $z$ for this cyclone becomes an $n\times 2$ matrix that describes successive positions of the cyclone in time, and can be written as:

where $T$ is the $n\times 3$ matrix whose elements $T_{jm}$ are given by $T_{jm}=t_j^{m-1}$, $\beta$ is the $(3\times 2)$ matrix of regression coefficients:

and $ϵ$ is the $(n\times 2)$ noise matrix. To generalize to a higher order regression model $(p>2)$ is straightforward: $T$ becomes an $n\times (p+1)$ matrix and $\beta$ an $(p+1)\times 2$ matrix.

This model of cyclone trajectories can then be used to describe the full data set of trajectories (all with length n) as a linear combination of k component pdfs, each of the form:

where the matrix of regression coefficients $\beta$, the covariance matrix $\Sigma$ for each component, and the weights given each component are to be determined, while the length of the trajectory n is fixed. Here the vertical bars denote the determinant, the superscript T the transpose, and $Tr$ means the trace. This is very similar to the component pdf written in the section on the “Mixture Modeling” approach, with some modifications necessary due to the matrix nature of the states $z$.

Hidden Markov Models

This approach seeks to determine persistent states in the atmosphere (Franzke, Crommelin, Fischer, & Majda, 2008; Franzke, Woollings, & Martius, 2011). The method identifies a pre-specified number of underlying multivariate Gaussian pdfs (called output distributions), and associates each pdf with an underlying “hidden” discrete state X. The evolution between these discrete states is governed by a Markov (first-order autoregressive) chain. This corresponds to the assumption that the hidden state transitions are governed by a stationary process.

The identification of persistent states (clusters) depends on the character of the matrix $A$, which describes the transition between the discrete states $X$. This transition probability matrix has elements $A_{ij}$: if the system resides in state $X_i$, the probability of transitioning to hidden state $X_j$ in a fixed time interval is $A_{ij}$. Note that the diagonal elements give the probability of remaining in a certain state. From this matrix it may be possible to identify certain groups of hidden states for which transitions are very likely within the group, but unlikely between a state in the group and one outside of it. Such groups are termed “metastable” and correspond to persistent states. Figure 2 gives a simple illustration of this, while figure 3 illustrates how a number of Gaussian pdfs, each associated with a hidden state, combine to give a distinctly non-Gaussian pdf for the observed system. Note the similarity here with the mixture modeling method.

Open in new tab

Figure 2. Schematic of Hidden Markov model of clusters. Each dark oval is a hidden state, corresponding to a cluster (or persistent regime). Each regime is associated with a probability distribution function (pdf) of an observed variable, which is given in the curve below the regime. The black numbers indicate the probability that each regime will persist for one time unit. The colored numbers indicate the probability of a transition between regimes.

Open in new tab

Figure 3. The three probability distribution functions shown are each associated with a hidden state, or cluster. The final pdf of the observed variable is a combination of the three component pdfs.

The main statistical problem in constructing a Hidden Markov model is the estimation of the elements of the transition matrix as well as the means and variances of the multivariate Gaussian output pdfs. This is accomplished by use of the expectation maximization algorithm. In addition, because each hidden state is associated with a range of possible observed states, one would like to estimate the most probable hidden state sequence. The technical approach to the statistical approaches, including estimates of statistical significance, are covered at length in Franzke, Crommelin, Fischer, and Majda (2008) and Franzke, Woollings, and Martius (2011).

This approach yields the possibility of a non-Gaussian (even multimodal) pdf of all the states even though each hidden state is associated with a Gaussian pdf (see figure 2 for an illustration). This feature is also seen in mixture modeling. However, the time sequencing of the observed states, which plays no role in the mixture modeling, is very important in the Hidden Markov model.

Non-Stationary Clustering; Finite Element Clustering

The great majority of approaches to cluster analysis assume that the statistics of the variables considered are stationary in time, so that there is no underlying slow change in such basic quantities as first and second-order moments. However, long-term changes in the climate system are known to occur, due to for example the changes in greenhouse gases or volcanic eruptions. Non-stationary clustering methods are designed to take into account these possible influences of time-dependent forcing functions on the preferred states. In the finite element, vector autoregressive factor approach (Horenko, 2010), the sequence of observed states $z$ is modeled by a discrete process in time t, in general form:

Ignoring for the moment the subscript $X$, this equation represents a general, discrete auto-regressive model, plus a term $u_t$ representing time-dependent external forcing, plus noise $ϵ$.

The autoregressive model uses multiples of the time lag $\tau$ up to $m$, the memory depth. The term $\mu$ represents the mean. However, the set of parameters $(\mu ,A_q,B,C)$ depends on which cluster $X$ the system is in, and this affiliation of the system with a particular cluster, thus a particular set of parameters, depends in general on time. This is illustrated in the schematic in figure 4.

Open in new tab

Figure 4. Schematic of the finite element vector autoregressive model. The time evolution of the state z(t) is given schematically in the black line. Two clusters are illustrated by the red and blue lines; z(t) is described by either the upper or lower equation, depending on which cluster the system is in at a given time.

Determination of the parameters for each cluster involves the adaptive finite element technique, while the order of the model (values of $m$ and the number of clusters) can be determined via the Akaike information criterion (Horenko, 2010).

Compared to the Hidden Markov model, this technique includes a more general treatment of the evolution of the observed states, including the effects of external forcing, which may have a timescale that is very long compared to the shortest timescale sampled in the data. The disadvantage of this technique is that more parameters have to be estimated compared to the Hidden Markov model.

Statistical Significance

One general approach to the assessment of statistical significance of multiple maxima (multimodality) in the direct estimate of a probability distribution function, or in the results of a cluster analysis, is to rely on the use of synthetic (also called surrogate) data sets (Nitsche, Wallace, & Kooperberg, 1994; Hsu & Zwiers, 2001; Christiansen, 2007; Straus et al., 2007; Franzke et al., 2008). For example, assume the original (physical) data set is expressed in terms of principal components. Then each of the N original maps (or states of the system) is represented as a vector in the PC-space of dimension M, where M is the number of PCs retained. The entire data set thus consists of M time series, each of length N. Each synthetic data set is constructed to have the same number N of original states as the physical data set, but with the time series of each of the M PC components replaced by a stochastically generated Gaussian time series, generated so that the lag-correlation statistics are close to the observed statistics for that PC component. By construction the time series of the various components are statistically independent of each other. Such synthetic data sets are constructed so that any deviations from a multi-normal distribution can be due only to sampling error.

The use of many such synthetic data sets can help establish the confidence of either the occurrence of multiple maxima (multimodality) in the pdf. For example, estimating the pdf in two dimensions with the kernel density estimation, one can increase the kernel smoothing h until only an acceptable percentage of the synthetic series show multiple maxima (Corti et al., 1999). On the other hand, one may rely on more stringent tests to distinguish the entire pdf from a multi-normal pdf, such as testing the marginal distributions against a normal distribution with the Kolmogorov-Smirnov test (Stephenson et al., 2001).

For the k-means method, one can similarly establish the confidence level by the percentage of synthetic data sets for which the variance ratio is larger than that associated with the physical data set (Straus et al., 2007). Often one finds that the confidence level is high for all values of k greater than some threshold, because the variance ratio generally increases with k (becoming infinite when $k=N$ [each state its own cluster])! There is no accepted definition of the best or optimal value of k, which may depend on the size of the data set and the use to which the clusters are put. For example, the reproducibility of results to subsampling tends to decrease with increasing k, and could be used to determine an upper limit.

One source of difficulty in assessing significance is the relevance of a Gaussian distribution as the relevant background against which to test, because it is known that many geophysical variables have skewed distributions (White, 1980). Skewness can be taken into account in generating the stochastic surrogate data sets (Christiansen, 2007), and generally lowers the significance.

An important issue in the significance testing of clusters is the size of the sample of data. Often the available data record is just too short to rigorously confirm significance (Stephenson et al., 2001), yet the clusters obtained may be consistent with other statistical analyses or may provide a useful coarse-graining of the original data for various purposes. This is discussed in the section on “Physical Significance and Usefulness of Clusters.”

A second, often critical, issue is the choice of how the physical data is chosen and prepared for clustering. To the extent that the user has some a priori knowledge of a particular set of structures that dominate fluctuations on a certain timescale, filtering of the data to remove much fast timescales (which may represent noise) will boost the significance of the clusters. For example, in seeking the preferred states associated with intra-seasonal fluctuations (roughly 10- to 90-day periods), the failure to filter out rapidly evolving synoptic scale, baroclinic weather systems can degrade the significance (e.g., Fereday et al., 2008). For static clustering techniques, significance can also be boosted if only quasi-stationary periods in the data record, periods in which the states change more slowly than usual, are first isolated (e.g., Toth, 1992; Straus et al., 2007).

Physical Significance and Usefulness of Clusters

The existence of preferred regions of the state space of the atmospheric circulation is suggested by studies of a hierarchy of dynamical models truncated to retain only a few modes. In the most highly truncated models, steady states were identified and associated with preferred states of the real atmosphere (Charney & DeVore, 1979; Charney & Straus, 1981; Reinhold & Pierrehumbert, 1982). States that characterize minima in the tendency of large-scale fields were identified in less highly truncated barotropic (Anderson, 1992; Branstator & Opsteegh, 1989) and two-level baroclinic models (Hannachi, 1997). In the intermediate three-level model of Itoh and Kimoto (1996), remnants of unstable periodic cycles could be identified in the full dynamical attractor, suggesting a physical basis for preferred states.

The clusters of intra-seasonal variability for the more complete (primitive equation) general circulation models (GCMs) with full suites of physical parameterizations generally compare well with those obtained from reanalyses, lending support to the physical nature of the results (Dawson & Palmer, 2015; Swenson & Straus, 2017). Further, the importance of the sample size in establishing formal significance is underlined by cluster analysis of large ensembles of GCM simulations and forecasts: the significance of the k-means clusters of intra-seasonal variability tends to be much higher for GCMs compared to reanalyses because the sample size is much larger. This indicates that the GCM results are robust across ensemble members.

Cluster analysis has also been proven to be a useful tool in a number of topics of weather and climate research, including not only understanding the midlatitude response to tropical forcing, but also in downscaling and extremes, model evaluation, weather forecasting, and synoptic climatology, all of which are reviewed in Hannachi et al. (2017). Downscaling is required for climate impact assessment. Statistical downscaling is, in some sense, the inverse of parameterization used in numerical weather prediction, in that it uses large-scale flow features and circulation regimes (clusters) as predictors to infer smaller-scale processes, notably precipitation and surface temperature at regional and local scales.

The large-scale atmospheric circulation patterns obtained from clustering can also be used for ocean downscaling. In particular, winter and summer circulation regimes over the North Atlantic–European region have been used as predictors to reconstruct ocean surface variables consisting of surface winds and temperature, which can then be used to force ocean models (Barrier, Cassou, Deshayes, & Treguier, 2014).

The link to surface weather and climate extremes is another application of clusters (Park et al., 2014; Xie et al., 2017; Zhao et al., 2017a, 2017b; Amini & Straus, 2018). The consequent persistence of these large-scale and synoptic patterns yields, through a cumulative process of weather fluctuations, an amplitude increase leading to extremes. A related application of regimes (clusters) to local weather is the endeavor known as synoptic climatology (James, 2006; Coleman & Rogers, 2007). Here it is the synoptic weather patterns themselves that are classified into clusters representing preferred synoptic structures.

The very nature of the preferred (and often persistent) large-scale circulation regimes obtained by cluster analysis suggests a link to weather forecasting and predictability. Operational forecast centers thus use clustering to summarize the wealth of dynamical information available from large ensembles of forecasts to indicate scenarios, or alternate distinct possibilities suggested by the forecasts (Ferranti & Corti, 2011; Ferranti, Corti, & Janousek, 2015).

Hierarchical Method

This is a very basic method that yields a hierarchy of representative patterns or maps. This method is very good for getting a feel for the way in which preferred patterns emerge from the data, and how they combine to form a few distinct clusters. Moreover, the dendogram (graphical tree representation of how clusters merge, and at what level of dissimilarity) provides some guidance as to the choice of how many clusters to retain. However, other than establishing reproducibility among many subsets of the data, overall significance is hard to establish.

Partitioning Methods

These methods include both k-means and a variant, clustering around medoids. These approaches also yield distinct preferred patterns, and assignments of each original state (or map) to one of the k patterns. Because for each value of k there is a well-defined metric of how clustered the data are (the variance ratio), synthetic data sets can be used in a straightforward manner to establish the statistical significance. In addition, the user can decide the null hypothesis (nature of the background pdf) against which the data are to be tested. For example, is the null hypothesis represented by a multivariate Gaussian pdf, or a multivariate pdf that includes skewness (or other properties)? Similar to the hierarchical method, the choice of the final number of clusters is not unique. This method is useful for clustering ensemble forecasts, by picking k to be the lowest value for which a threshold significance is achieved, thus allowing for the possibility that there are not multiple preferred states in the forecasts.

Fuzzy Clustering

This method is similar to the partitioning methods, but assigns a probability of each original state belonging to more than one cluster. Because this approach identifies states that are not particularly close to any single cluster centroid, it provides a method for filtering out those states that might be considered ambiguous with regard to the clusters. The control the user has for the degree of smoothing in assigning probabilities is both an advantage (sensitivity can be established) and a disadvantage, because it is difficult to establish a statistically optimal level of smoothing.

Self-Organizing Maps

The goal of this approach is quite distinct from other methods: it seeks to determine a relatively large (of order 25) set of representative patterns that are ordered by their similarity, and so provides a coarse-grained view of the continuum of possible states. This is useful if one is interested in studying the evolution of a particular regional field (such as sea-level pressure over the Pacific) during the annual cycle, or as function of tropical forcing on intra-seasonal and inter-annual time scales (Johnson & Feldstein, 2010). Because the goal is not the identification of just a few preferred patterns, but the sense of evolution of the entire field, the arbitrariness in defining the number of representative patterns is not of concern. However, the assessment of statistical significance is problematic.

Mixture Modeling

This technique is useful if a statistical model for the entire pdf is sought. The pdf is modeled as a sum of multivariate Gaussian component distributions, with each component distribution interpreted as a cluster. Every observed state is associated with multiple components (clusters), each with a certain probability. Mixture modeling is a particularly efficient approach to modeling the entire pdf, in the sense that a relatively few components can describe a variety of pdfs (given that each component has a number of parameters to fit). This strength in pdf modeling, however, suggests a certain caution in using this technique as a cluster analysis technique: because each component does have many parameters, the number of clusters (components) is likely to be quite small in order to avoid over-fitting.

Hidden Markov Model

This approach can be very useful in identifying preferred states that are also persistent, as well as the most likely transition paths between these states. As with many other approaches, programming packages on various popular platforms are available. However, because the underlying assumptions made in identifying the number of preferred states and their transition properties, as well as significance testing, are not transparent, it is recommended that the user have some prior knowledge and physical intuition with regard to the system under study before application of this method.

Finite Element, Vector Autoregressive Factor Model

This is the only approach discussed that can explicitly take into account non-stationarity, and so may be essential in some applications. Because it, like the Hidden Markov model, is a relatively new approach, the climate community has not had as much experience in its application. Thus it is recommended that the user obtain some some prior knowledge and physical intuition with regard to the system under study before application of this method.