Introduction

Nowadays, skillful numerical weather prediction is limited by about 10 days due to the chaotic nature of atmospheric dynamics [1]. Skillful seasonal forecasting typically relies on the predictability of slow-varying components of the climate system, such as sea surface temperature, sea ice, snow cover, and soil moisture. For instance, the predictability of El Niño Southern Oscillation (ENSO) phenomenon is a remarkable example o f high skills of the seasonal forecast system [2,3]. However, recent studies demonstrated that the maximum seasonal forecast skills have not yet been achieved, and pointed to the stratosphere as a potential source for enhancing seasonal predictability [1,4].

Before [5] the stratosphere was considered playing a passive role in the stratosphere-troposphere coupling, that is, it does not influence troposphere dynamics. Baldwin and Dunkerton showed [5] that the strength of the polar vortex affects the main features of the paths of the cyclones propagation. Fro m that time, interest in the stratospheric dynamics and its climatic changes have constantly been growing.

Polar stratosphere became extremely cold during the polar night, and meridional temperature gradient becomes strong, which leads to the formation of the polar vortex. Planetary waves propagating from the troposphere to the stratosphere disturb and sometimes destroy polar vortex. These events are known as Sudden Stratospheric Warming events (SSW), as the temperature near the pole dramatically increase (up to 40° per 4 -7 days) when the vortex is destroyed. Two types of SSW are classified nowadays: "displacement" with the center of the vortex shifted significantly towards the equator, and "split" with the vortex split into two vortices. There are periods of the extremely strong vortex as well. Variability of the PV intensity is the most influential factor of intraseasonal variability in the winter stratosphere. During springtime as the polar day is coming to high latitudes, temperature gradient decrease, so the polar vortex disappears. Summertime variability of the stratosphere dynamics is low, and there is no source of the stratosphere-troposphere coupling.

It has been shown that weakenings of the polar vortex precede the shift of the storm tracks (main p ath of midlatitude cyclones propagation) to the south, which may cause cold outbreaks in the North Atlantic -Europe region [6][7][8][9][10]. These anomalies may act in the troposphere up to 2 months [6,11]. Amplitudes of these anomalies are comparable to the effect of the ENSO [3]. Therefore, one can extend and improve long-term weather forecast [12][13][14][15][16]. However, relatively little attention is paid to strong polar vortex events compared to SSW. One of the main reasons for that is the difficulty of their identification classification.

There are a few metrics for the describing states of the PV [3,17]. Widely used key features of PV for applying different types of analysis are the aggregated and diagnostic parameters like maximum pressure anomaly, zonalmean zonal winds, etc. These parameters do not preserve spatial characteristics like PV center shift from the North Pole, vortex shape parameters or various anomalies duration.

Over last decades machine learning methods demonstrated spectacular results in research of climatic changes of atmospheric circulation [18][19][20][21][22]. Researchers mostly rely on state-of-the-art clustering methods. Using the resulting atmosphere states grouping one can assess atmospheric circulation characteristics and trends within each cluster. Most of this kind of works are focused on troposphere states research [18][19][20]22]. However, machine learning techniques are shown recently to be fruitful for stratosphere states clustering [21]. Mostly Kohonen selforganizing maps [23] is applied in these works for clustering. However, its ability to preserve the topology of a dataset is rarely used. Hierarchical agglomerative clustering is the less frequently used method [18,21]. For this type of clustering, the source data of geophysical fields are usually a ggregated following a researcher sense of a particular operation ability to preserve the informational content. As a result, each PV state is represented with vectors. This procedure leads to the loss of crucial information about the spatial configuration of objects being researched.

In our work, we focused on polar vortex states research with the use of geopotential height at 10hPa level (see section 2.1 "Data and preprocessing"). We applied hierarchical agglomerative clustering method on the embedded representations of PV states that are generated by variational convolutional autoencoder with the sparsity constraint (hereafter SpCVAE) [24,25]. SpCVAE as a feature extractor and a tool for dimensionality reduction preserves features of PV spatial configuration, yet it is capable of decreasing the computational costs of clustering. Purposes of a SpCVAE here are dimensionality reduction and extraction of significant information based on the whole PV states dataset. Being a neural network trained end-to-end the SpCVAE avoids manual feature engineering. Thus there is no need to rely on a researcher`s sense of the importance of PV aggregated parameters.

2

Data and Methods

Data and Preprocessing

We analyzed geopotential heights (HGT) and potential vorticity (PVt) fields at 10 hPa level from JRA-55 (Japanese 55-year Reanalysis) in this study [26]. We considered wintertime period December-February (DJF) for 1958 -2014 years. The spatial resolution of the data is 1,25° х 1,25°, the upper level of the model is 0.1 hPa, which is crucial for the analysis of stratospheric processes. JRA-55 shown to be in good coherence with all modern reanalysis data (S-RIP [27]). The main advantage of JRA-55 is the extended period from the 1958 year in comparison with ERA-Interim (European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis data) and MERRA (NASA Modern Era Reanalysis for Research and Applications) that start from 1979.

Source data for describing PV states were considered to be field values to the north from 40 °N. At the preprocessing stage, timestep-defined snapshots of HGT and PVt fields were projected using North-polar Lambert azimuthal projection and interpolated to form a two-dimensional flat matrix of size 256x256. For each date of a year, we calculated the climatological median and subtracted it from each snapshot of this date (e.g., from all snapshots within January 25th of each year). Since we use the North-polar Lambert azimuthal projection, only central rounded part of each sample is informative, so during all calculations, we used the mask 𝑀 𝑖𝑗 (Fig. 1b).

We further normalized these snapshots, so all their values are limited and take values between 0.0 and 1.0:

𝑥 𝑖𝑗 = 𝑥 𝑖𝑗 −min 𝑖,𝑗 𝑋 max 𝑖,𝑗 𝑋−min 𝑖,𝑗 𝑋 ,(1)

where 𝑋 denotes the whole dataset snapshots, 𝑥 is a particular snapshot, 𝑖 and 𝑗 are 𝑥 matrix indices, min(𝑋) and max(𝑋) are calculated taking the mask 𝑀 𝑖𝑗 into account. This kind of normalization applied to both HGT and PVt datasets separately. With this preprocessing procedure, the PV states dataset is represented with two fields containing 21476 matrices of size 256 x 256 which values fit the range between 0.0 and 1.0. Some examples of this dataset are shown in Fig. 1a.

Sparse Variational Convolutional Autoencoder

Convolutional autoencoder (CAE) is a particular type of neural autoencoders [24,25] which purpose is to deal with two-dimensional fields when it is crucial to take into account spatial relations between data anomalies. The fundamental purpose of autoencoders is the approximate copying of input data with some constraints imposed on the network structure. These constraints are often the restriction on the maximal dimensionality of the embedded representation obtained on the encoder output (see Fig. 2), and the CAE is called undercomplete. In this sense, the whole undercomplete CAE is essentially a transformation: ℱ: 𝐴 → 𝐻 → 𝐴 ̃, where 𝐴 is an input example matrix, 𝐴 ∈ 𝔸, 𝐻 is its hidden representation vector such as 𝐻 ∈ ℋ where ℋ is a hidden representations space, and 𝐴 ̃ is the reconstructed example matrix, 𝐴 ̃∈ 𝔸. Here 𝔸 is the space ℝ 𝑛 where 𝑛 is the number of pixels of input examples. In our study, 𝑛 = 2 * 256 * 256 since HGT and PVt projected examples are matrices 256x256. The transformations 𝐴 → 𝐻 and 𝐻 → 𝐴 ̃ are referred hereafter as encoder and decoder respectively. An autoencoder is trained with the loss function defined according to the similarity definition suitable to the problem. Mean squared error is a widely used loss function that is suitable for most tasks that imply the processing of geophysical fields. In a case of limited source data values, one may normalize them accordingly to use binary cross-entropy loss (BCE) (2). We normalized source data, so its values are limited (see Section 2.1) therefore we use BCE loss ℒ 𝑏𝑐𝑒 :

ℒ 𝑏𝑐𝑒 (𝑋, 𝑝) = − ∑ ∑ 𝑀 𝑖𝑗 (𝑥 𝑖𝑗 ln ℱ(𝑥 𝑖𝑗 )) 𝑛 𝑗=0 𝑚 𝑖=0 ,(2)

It was shown [24,25,28] that successful training of a neural network implies tuning its weights the way that the trainable part of the network extracts latent parameters distribution of the training dataset. With the trained CAE, the embedded representation [29,30] of input samples preserve enough information for the network to be able to reconstruct it with the appropriate quality. We use this feature of CAEs to perform dimensionality reduction with minimum loss of information about latent parameters variability and the spatial structure of PV states.

We applied commonly used techniques of convolutional neural networks quality improvement called Transfer Learning (hereafter TL) [28,[31][32][33][34][35] and Fine Tuning (hereafter FT) [36,37]. In practice, TL implies the construction of a new neural network based on a subset of layers of the network that was previously trained on a huge dataset, e.g., ImageNet [38]. Moreover, these layers are set to be "frozen," that is, their weights are not optimized during training. This approach is fruitful when one uses a dataset which statistical characteristics are similar to ImageNet. In practice, it means that the new dataset images should contain visual patterns similar to ones that are frequently met in ImageNet. TL application significantly decreases the computational costs of new models training. Fine Tuning approach implies turning off the "frozen" state for some top layers of the transferred neural network. With this approach, one can tune the whole CAE taking into account the peculiar properties of the dataset.

We applied the TL technique while building the encoder part of the CAE. We used pre-trained VGG-16 [39] as a transferred network. We used the convolutional part of VGG-16 and a set of new fully-connected layers attached to it. VGG-16 convolutional sub-network is denoted as "convolutional core" in fig. 2. We also applied several regularizations to prevent overfitting. Particularly we used the dropout [40] approach and L2 regularization which penalizes high-value weights.

In our model, the VGG-16 output is reshaped to a vector which is then the input for the encoder fully-connected part. This fully-connected part consists of three layers FC1, FC2, and FC3, which are alternating with dropout layers. FC3 is an intermediate one between encoder and decoder. The output of FC3 is the hidden representation of the input sample and the input vector for the decoder.

Autoencoders usually tend to be symmetric. With this approach, we constructed the decoding part to be mirrored to the encoder. All the weights of the decoder are trainable. Following the best practices of composing convolutional autoencoders, we used two-dimensional upsampling layers which mirror max-pooling layers of the encoder. Upsampling here is an operation of repeating the layer`s input rows and columns by the specified number of times. Decoder outputs are the reconstructed HGT and PVt fields of input example which has to be similar to the input in a sense defined by the loss function of the network which is BCE loss (2).

The dimensionality constraint mentioned above is the dimensionality of ℋ . The common approach of the clustering involving autoencoders relies on their capability of projecting the input examples to the hidden representation space ℋ. This transformation was shown to be trained so that the examples which are close to each other in ℋ are similar [29]. However, the key feature of a dataset should be the opposite for the reliable, stable, and reproducible clustering, that is, similar examples should be located close to each other in ℋ . The ordinary undercomplete CAE does not guarantee this property of the projection 𝔸 → ℋ. This issue might be addressed with the variational autoencoder (VAE) [41] which was shown to produce continuous latent variables space ℋ. That is, with VAE involved, similar examples are located close to each other in ℋ. However, since the distribution of the features of 𝐻 is normal in case of VAE, the clustering cannot produce valuable results in the generated feature space ℋ. In our study, this issue is addressed with the constraint of sparsity, that is, features of the vector 𝐻 are forced to be Bernoulli-distributed. With this constraint, the hidden representation vectors 𝐻 tend to be sparse, that is, only a few features are non-zero for a particular example 𝐴. The undercomplete CAE with the mentioned constraints imposed is referred hereafter as sparse variational convolutional autoencoder (SpCVAE). Its structure is presented in Fig. 2. As shown in Fig. 2, the fields of the input example are processed separately, similar to the approach applied in [42]. In our work, the only hyperparameter of the proposed SpCVAE is the number of nodes of the FC3 layer. This number at the same time is the number of features of the hidden representation 𝐻 (see Fig. 2) and the dimensionality of ℋ (hereafter 𝐻𝐷𝑖𝑚). There is a trade-off between reconstruction quality and the 𝐻𝐷𝑖𝑚. We use the multiscale structural similarity (𝑀𝑆𝑆𝑆𝐼𝑀 [43]) as a measure for reconstruction quality. For optimization reasons, we use the metric (1 − 𝑀𝑆𝑆𝑆𝐼𝑀) with the rule "less is better" (Fig. 3a). We have conducted the research of the reconstruction quality versus the 𝐻𝐷𝑖𝑚 (see Fig. 3a). There is a reasonable value of 𝐻𝐷𝑖𝑚 = 96 since after this value the reconstruction quality stops improving significantly. There are more candidates for the best choice of 𝐻𝐷𝑖𝑚 , however, only starting with the 𝐻𝐷𝑖𝑚 = 96 the clustering results become stable and reproducible. Taking this result into account, we further used the SpCVAE with 96 neurons of layer FC3. We have trained this SpCVAE using the data described in section 2.1. We use then the encoder output in the inference mode of the trained SpCVAE as PV states representation of reduced dimensionality.

Hierarchical Agglomerative Clustering

We applied Lance-Williams hierarchical agglomerative clustering [44][45][46] to define groups of stable PV states. This method is frequently used for atmosphere and stratosphere states clustering [18][19][20][21][22]. We applied this clustering method to PV states objects described with low-dimensional hidden representations generated by SpCVAE (see Fig. 2). We considered the Euclidean metric as a distance between objects in this feature space. We used Ward minimal inter-cluster distance [47] as a criterion for clusters union. Ward inter-cluster distance between clusters U and V is the following: (6) c. for each 𝐶 ∈ 𝑆 𝑡+1 , calculate inter-cluster distances to 𝐷 𝑊 (𝐶, 𝑊) using ( 3). Agglomerative hierarchical clustering procedure with Ward inter-cluster distance definition ( 3) is the one demonstrated most valuable results in a set of synthetic clustering problems [48]. With this procedure, the only hyperparameter is the target clusters number 𝐾. We use the empirical "elbow rule" to choose the best clusters number. Additionally, we considered topic-specific reasoning: we wanted the clustering method to be capable of discriminating weak PV states of types "displacement" and "split" yet to be capable of discriminating various shifted states that were demonstrated recently to be present [21]. We considered mean Silhouette score (hereafter 𝑆𝑠𝑐𝑜𝑟𝑒) as a measure for clustering quality and used it for "elbow rule". We calculated 𝑆𝑠𝑐𝑜𝑟𝑒 using Euclidean metric of the hidden representation feature space ℋ. For each 𝑖-th object of the 𝑗-th cluster 𝑥 𝑖 ∈ 𝐶 𝑗 : 𝑎(𝑥 𝑖 ) is the mean distance between 𝑥 𝑖 and all other objects of cluster 𝐶 𝑗 , 𝑏(𝑥 𝑖 ) is the mean distance between 𝑥 𝑖 and all the objects of all other clusters:

𝐷 𝑊 (𝑈, 𝑉) = ‖𝑈‖‖𝑉‖ ‖𝑈‖+‖𝑉‖ 𝜌 2 (∑ 𝑥 𝑈 ‖𝑈‖ 𝑈 , ∑ 𝑥 𝑉 ‖𝑉‖ 𝑉 ) ,(3)𝑎(𝑥 𝑖 ) = 1 ‖𝐶 𝑗 ‖ ∑ 𝐷(𝑥 𝑖 , 𝑥 𝑘 ) 𝑥 𝑘 ∈𝐶 𝑗 ,𝑘≠𝑖(7)𝑏(𝑥 𝑖 ) = 1 ‖𝒯‖−‖𝐶 𝑗 ‖ ∑ 𝐷(𝑥 𝑖 , 𝑥 𝑘 ) 𝑥 𝑘 ∈𝒯\𝐶 𝑗 , (8)

where 𝒯 is the whole dataset, ‖𝒯‖ is the number of its elements, 𝐶 𝑗 is the cluster that 𝑥 𝑖 is assigned to, ‖𝐶 𝑗 ‖ is the number of its elements, 𝐷(𝑥 𝑖 , 𝑥 𝑘 ) is the function defining the distance between 𝑥 𝑖 and 𝑥 𝑘 . 𝐷(𝑥 𝑖 , 𝑥 𝑘 ) is the Euclidean distance in our study. With these notations 𝑆𝑠𝑐𝑜𝑟𝑒 for one object 𝑥 𝑖 is given by: 𝑠(𝑥 𝑖 ) = 𝑏(𝑥 𝑖 )−𝑎(𝑥 𝑖 ) max(𝑏(𝑥 𝑖 ),𝑎(𝑥 𝑖 )) ,

and mean 𝑆𝑠𝑐𝑜𝑟𝑒 is given by:

𝑆(𝒯) = 1 ‖𝒯‖ ∑ 𝑠(𝑥 𝑖 ) 𝒯 .(10

) 𝑆𝑠𝑐𝑜𝑟𝑒 for each 𝑥 𝑖 is the measure of its similarity to the cluster 𝐶 𝑗 and yet its dissimilarity to all the other elements outside 𝐶 𝑗 . So the higher 𝑆𝑠𝑐𝑜𝑟𝑒, the more 𝑥 𝑖 is similar to 𝐶 𝑗 and the less similar to other clusters. Therefore, mean Silhouette score (10) is considered as a measure of clustering quality with the rule "higher is better". Even though the maximum mean 𝑆𝑠𝑐𝑜𝑟𝑒 is achieved with two clusters, more reasoning should be involved when one is choosing the number of clusters. First, there should be observed "split" and "displacement" SSW events. There are also should be The Use of Modern Information Technology for Research of Technical and Natural Systems ________________________________________________________________________________________________ 57 observed at least one strong pole-centered state, and some other shifted states which were presented recently [21]. For each number of clusters more than 3 we inspected maps of mean geopotential heights at 10hPa level. We have chosen the minimal number of clusters that let us discriminate weak states of PV of types "displacement" and "split". This discrimination is observed starting from 𝐾 = 7. With this reasoning, the group "L" of clustering results (see Fig. 3b) should not be considered as an option. There is also group "N" of clustering results which are characterized by low mean Sscore, that is, low clustering quality. The results within the "M" group are characterized by too high clusters number or low clustering quality. In our study, the group "region of interest" is considered as a group of promising clustering results (Fig. 3b). The numbers of clusters which produces high average 𝑆𝑠𝑐𝑜𝑟𝑒 within this group are 12 and 13. In our study, we use 𝐾 = 12.

Summarizing the proposed method for clustering states of PV here is its general structure: 3).

Results and Discussion

We applied the approach presented in Section 2 to the dataset described in Section 2.1, "Data and preprocessing." PV states were clustered using 𝐾 = 12 number of clusters. We calculated the map of mean geopotential heights at 10hPa level for each cluster. These maps are shown in Fig. 4. We also selectively inspected individual PV states represented by HGT fields. Visual inspection of these examples shows that PV states grouped by the proposed method are similar to each other within each cluster. SSW events are clearly seen in the composites in Fig. 4: cluster 2 for "split" type and clusters 1 and 3 for "displacement" type. We have conducted a more detailed analysis of the clusters 1-3. The calendar in fig. 6 presents the periods associated with the states of the clusters 1-3 and concurrently the known SSW events identified by experts. Central

The Use of Modern Information Technology for Research of Technical and Natural Systems ________________________________________________________________________________________________ 58 dates of the "split" and "displacement" SSWs obtained by Charlton and Polvani for period 1958 -2002 [49], and non-classified SSWs from [50] for period 2003 -2013 are shown. Almost all expert-defined SSW events are clearly colocated in time with the segments of clusters 1-3. Most of "split" SSW events are co-located with the cluster 2 or with the segment consisting of states associated with two clusters including cluster 2. The only one SSW event in 2002 missing corresponding states of clusters 1-3 is a subject for further research.

Central dates of the SSWs are defined as the dates when zonal-mean zonal winds at 10hPa and 60°N fall below zero m/s (became easterly). In Fig. 5b, we present zonal wind averaged along 60°N for each cluster. In this figure, clusters 1-3 are distinctively characterized by winds close to zero or even easterly winds. This behavior is consistent with the nature and definition of SSW events [17]. Since there was no expert-level knowledge involved during clustering, we may consider the proposed method to be capable of objective weak vortex clustering.

It is also clear that clusters 10 and 12 represent strong vortex centered on the pole, and only a slight shift is observed for cluster 12. In Fig. 5a the frequency histogram is shown for the clusters obtained in this study. Cluster 10 is the most frequent state of PV, which is consistent with the current understanding of the nature of PV. Clusters 4-9 and 11 represent shifted PV states of different intensity, which may be estimated by the zonal mean zonal wind at 60°N shown in Fig. 5b. In accordance with the recent study [21], there are different shift directions: towards the Atlantic (clusters 7 and 11), towards Eurasia (clusters 4-6 and 8) and North America (cluster 9).

Conclusion

We propose a sparse convolutional variational autoencoder as a dimensionality reduction tool. With this model, we extracted valuable features of PV states initially represented by geopotential height and potential vorticity at 10hPa level. Using the representation of reduced dimensionality, we applied the Lance-Williams hierarchical agglomerative clustering with Ward inter-cluster distance definition. This method for the first time is capable of discriminating weak PV states of types "displacement" and "split," which is crucial for the analysis of the stratosphere-troposphere interactions. The proposed method is also capable of classifying stable states of strong PV characterizing by different directions of its center shift. This classification is found to be physically valid and consistent with recent studies. The presented clustering method the first time provides an opportunity to analyze the influence of the strong PV characterized by various shift directions on characteristics of tropospheric circulation.

The proposed clustering method provides an opportunity of researching climatic changes of wintertime stratosphere, which is crucial for improving seasonal forecast skill and assessing the long-term variability of the climate system.

Results of this work can be used as a basis for new stratosphere-troposphere interactions research for existing and future climate scenarios.