Non-linear modes of global sea surface temperature variability and their relationships with global precipitation and temperature

Sea surface temperature (SST) modes are comprised of variability that arises from inherently nonlinear processes. Historically, a limitation arises from applying linear statistics to define these modes. Accurate depiction of the complex, non-linear nature of SST modes of variability necessitates the specification of a model capable of producing nonlinear patterns. In this study, we apply an artificial neural network algorithm integrated with autoencoders to analyze the seasonal non-linear global SST modes allowing for improved characterization of the modes and their large-scale temperature and precipitation teleconnections. Our results show that during boreal summer, SST cooling over the central to eastern tropical Pacific co-occurs with the Arctic amplification. In recent decades, the negative SST trend in the central to eastern tropical Pacific, combined with the positive trend in the western tropical Pacific is linked to an increase in the amplitude of SST modes associated with the Arctic warming, resulting in warmer temperatures over large portions of the global land, particularly over Greenland. In boreal winter, El Niño Southern Oscillation (ENSO) is the prominent global SST mode. The distinct spatiotemporal patterns of ENSO modes are associated with unique effects on regional land temperature and precipitation. The central Pacific El Niño is more associated with the combination of warm and dry conditions over Western Australia, and the northern part of South America. Conversely, the central to eastern El Niño is more associated with the combination of warm and dry conditions over parts of Southern Africa, and the northern part of South America. The spatiotemporal patterns and trends in the amplitude of the analyzed non-linear global SST modes alongside their regional influences on temperature and precipitation are discussed. The broader impact of this study is on the potential of neural networks in effectively delineating non-linear global SST modes and their associations with regional climates.


Introduction
Sea surface temperature (SST) modes are patterns of recurrent spatial and temporal variability existing at different time scales.The SST modes play a critical role in regulating the global climate system, influencing extreme weather events, global weather patterns, and the longer-term evolution of climate [1,2].A challenge in studying SST modes arises from finding the appropriate methods to represent their inherent nonlinearity arising from the relationships and interactions between oceanic variables, such as temperature, salinity, and currents, that might not be linearly related and can exhibit complex behavior [3].
El Niño Southern Oscillation (ENSO) is one of the most prominent modes of variability in the tropical Pacific.El Niño and La Niña phenomena are products of intricate nonlinear ocean-atmosphere feedback mechanisms [4][5][6].Studies such as [7,8] have highlighted the nonlinearity in ENSO, particularly in its amplitude, asymmetry, and the interaction between the zonal advection feedback and thermocline feedback.In the same paradigm, other studies have reported on the non-linearity of SST modes [9][10][11][12][13][14][15] and the effects of tropical SST variability and changes on regional climates [16,17].
Until recently, linear techniques, such as principal component analysis (PCA), have been a favored tool to define ENSO and other SST modes [18].Despite its usefulness in denoising a data set and extracting leading patterns of variability, PCA can depict only a linear space-time approximation of a non-linear process and might fail to depict important nonlinear interactions [19].Artificial neural networks (ANNs) offer an alternative to linear techniques for modeling complex, nonlinear systems [20].ANNs can approximate any linear or nonlinear function, provided they are defined with sufficient architecture, potentially making them well-suited for capturing the nonlinear nature of SST modes.
Despite the relatively recent use of ANNs to study SST modes, they have been applied mainly to limited spatial domains with encouraging results.For instance, [21] found that the application of a neural network to create a nonlinear PCA approach was successful in extracting periodic modes in the tropical Pacific.References [22][23][24][25] used convolutional neural networks to forecast the time amplitude and type of ENSO.Labe and Barnes [26] applied ANN to predict the onset of slowdowns in decadal warming trends of global mean surface temperature.Gordon et al [27] applied ANN to predict Pacific decadal oscillation (PDO) transitions with over 12 months lead time.Further, [28] showed that the results of ANN applied to the tropical Pacific Ocean region can be interpreted physically.Motivated by such findings, we extend the dataset domain to global oceans and apply an autoencoder type ANN to investigate if we can improve our understanding and characterization of the spatial aspects of global SST modes.This is crucial as the identification of non-linear SST modes can enhance our ability to predict climate variables related to such modes.This can have practical implications in areas such as seasonal climate predictions that are important for agriculture, water resource management, and disaster risk reduction.Additionally, a more accurate characterization of SST patterns can facilitate the improvement of state-of-the-art climate models.Hence the focus of this study is to explore the capability of autoencoder ANN in characterizing the non-linear spatiotemporal features of seasonal global SST modes and their influences on global temperature and precipitation.To our knowledge, this is the first application of autoencoder ANN to diagnose global modes of SST variability.
Several studies have utilized self-organizing maps (SOMs), a type of unsupervised neural network, for clustering and visualizing high-dimensional SST data, capitalizing on their ability to preserve topological properties and ease of interpretation [29][30][31].Although SOMs excel in pattern identification and simplicity, owing to their static architecture, SOMs have limited capabilities in representing hierarchical relations of the input data and may face challenges in effectively mapping and interpreting data when the number of dimensions or features in the dataset is very large and might also have difficulties capturing the deeper, multi-layered structures that often exist in high-dimensional climate data [32].Conversely, autoencoders can adapt their layers and neurons during the training process to better represent the data.Moreover, with their deep learning architecture, autoencoders are adept at encoding (global SST) data by reconstructing inputs from compressed forms.This feature enables them to capture both linear and non-linear relationships within data and to learn hierarchically structured features, presenting a more detailed representation of SST patterns.However, they can be prone to overfitting and require considerable computational resources.Despite SOMs offering straightforward, interpretable results, autoencoders can provide a more holistic understanding of complex and non-linear interactions in SST data.

Data
Gridded global SST data were obtained from the NOAA Extended Reconstructed SST V5 [33].The temporal resolution of the SST data is monthly from 1979 to 2022; the spatial resolution is 2.0 • longitude and latitude comprising 16 020 grid points.Global gridded observational temperature and precipitation data sets are obtained from the Climate Research Unit [34] at monthly temporal resolution and 0.5 • longitude and latitude.The global climate data sets are stratified into the boreal winter (January to March: JFM) and boreal summer (June to August: JJA).Hence the time series dimension is 132 per season.The SST data are de-seasonalized at the preprocessing stage by subtracting the long-term (1979-2022) monthly mean value of each month from the corresponding data values of the respective month.

Method
To derive non-linear global seasonal SST modes of variability, we utilized autoencoder-based ANN.Autoencoders are a specialized category of ANN designed for unsupervised learning and efficient data representations [35].Autoencoders are adept at extracting non-linear patterns from complex highdimensional data without the need for labels [36].Figure 1 shows the architecture of the autoencoder model utilized for the non-linear encoding of global SST.
In the pre-processing stage, before inputting the SST data into the autoencoder, the data are normalized to have values between 0 and 1, ensuring that all SST observations have a similar scale.This helps in speeding up the training process and faster convergence [37].Eighty percent of the total data was used to train the model and the remaining 20% was used for testing the model's performance on data that it has not seen during training, thereby providing an estimate of its generalization ability.Twenty percent of the training data was used as the validation subset-for hyperparameter tuning.
The autoencoder is defined as a neural network with an encoder and a decoder.The encoder learns to compress the input data into a lower-dimensional representation.In this case, it consists of a single dense layer with rectified linear unit activation [38].The decoder learns to reconstruct the original data from this compressed representation.The reconstruction error between the original input data and the reconstructed data provides the metric of how well the trained model has learned the crucial patterns of the global SST data.
The input and output neurons correspond to the 16 020 grid points of the input global SST data.We examined architectures of several neurons of different sizes (2, 4, 16, 32, 64, 128,…) aiming to choose the number of neurons that create a sufficient model balancing complexity against the risk of overfitting.This led to the selection of 64 neurons as we observed equivalent magnitudes of validation loss beyond 64 neurons suggesting that adding more neurons does not significantly improve the model's performance.Moreover, using additional neurons might lead to overfitting, whereas fewer might lead to underfitting.The trained encoder part of the autoencoder is used to transform the original SST data into a lowerdimensional space, aiming to capture the underlying non-linear modes of SST variability.Therefore, the autoencoder will learn to encode the 16 020 dimensional input into a 64 dimensional representation and then decode it back to 16 020 dimensions.
The next stage involves training the autoencoder to minimize the mean squared error between the original SST data and its reconstructed version, using the Adam optimizer [39].To avoid overfitting, the training was stopped at the epoch number (10) when the validation loss increases, which also corresponds to the instance when the rate of mean squared error decrease becomes insignificant in the training loss (indicating a learning plateau), implying that the model is not ascertaining new significant features [40].The Adam optimizer was used for training the model because of its adaptive learning rates, efficiency, and documented performance [39].
Additionally, the encoded SST patterns were compared to the patterns derived from PCA using the congruence coefficient (Ibebuchi and Richman [41]) to observe if the autoencoder reproduced the linear patterns, or captured new non-linear patterns, altering linear patterns that PCA might have failed to resolve properly.Congruence coefficients values that range from 0.98 to 1.00 represent an excellent match; 0.92 to <0.98 represent a good match; 0.82 to <0.92 represent a borderline match; 0.68 to <0.82 represent a poor match and values <0.68 represent a terrible match [41,42].
To obtain the time series associated with the encoded patterns, they were projected onto global SST anomaly data.Additionally, the time series were analyzed for periodicities using wavelets [43].After such analyses, the time series were regressed onto seasonal global temperature and precipitation to establish the teleconnective importance of these seasonal encoded patterns.

Results
The compressed seasonal spatial patterns captured by the first ten nodes (i.e. the encoded patterns) of the autoencoder, presented as standardized maps are presented in figure 2, showcasing the spatial variability of each node.Figure 3 shows the time series of the encoded patterns in figure 2. The patterns exhibit a range of spatial scales, revealing the underlying complexity of the SST data.Utilizing 64 neurons for the non-linear encoding of the SST (figure 1) also implies the potential of 64 encoded patterns in our study.However, since our study is exploratory-aiming to examine the potential of autoencoders in encoding global SST data-we concentrate on describing the patterns and teleconnective properties of the first ten encoded patterns, for brevity.Had this been a confirmatory study, with a target climate pattern or patterns, it would be crucial to scan through each of the 64 encoded patterns and further analyze only those patterns that optimally match with the target [44].
We assessed the degree of non-linearity of the autoencoder patterns by comparing them to SST patterns derived from rotated PCA, using the congruence coefficient.Only the first two encoded summer and winter patterns in figure 2 were identified to have limited similarity with the rotated PCA patterns, given their congruence matches in the range of 0.70-0.78(figures S1(a) and (b)), values associated with a poor match [41].This indicates that when the constraint of linearity is removed, these encoder patterns change from their linear counterparts.
The very limited similarity between the encoded and the principal component (PC) patterns is closest for node 1 and PC1-both resembling one of the ENSO variabilities, during JFM (figure S1(a), top panels).However, the node 1 encoder representation of El Niño is characterized by a stronger linkage to the South Pacific convergence zone than its PC counterpart, potentially making the encoder pattern useful to climate modelers in addressing the double inter-tropical convergence zone problem [45].Despite these differences concentrating on just the ENSO region, by projecting the patterns of node 1 and PC1 onto global SST anomaly data during JFM, respectively, the encoded patterns, as well as the PCAbased patterns, correctly identify ENSO events (figure S2).The ability of the autoencoder patterns to correctly discriminate ENSO events is quite promising and indicates that they can be interpreted physically.
Utilizing the congruence coefficient, figure 4 links traditional climate modes to the encoded patterns in figure 3. Node 1 is associated with SST anomalies over the central to eastern tropical Pacific (figure 2(a)) and is closest to the bivariate ENSO index (BET) [46] with The encoded patterns in figure 2(a), such as nodes 1, 2, 4, 5, 7 8, and 10 identify the different flavors of ENSO mode [47], marked by zonal dipole patterns in the tropical Pacific.However, figure 4(a) shows that nodes 1, 5, and 10 are highly related to traditional ENSO indices but not so much related to the other encoded patterns (figure 4(a)).Nodes 1 and 5 have the highest congruence matches with the BET, the eastern tropical Pacific (ETP), and the inter-decadal Pacific oscillation (IPO).Since BET and ETP are different ENSO indices, and IPO is an inter-decadal mode of ENSO, it is plausible that nodes 1 and 5 are identifying the underlying signals that are common to these phenomena.This could be due to the nodes detecting overarching patterns that are characteristic of both the short-term ENSO cycle and the longer-term IPO cycle.Compared to the PCA representation of ENSO, the capability of the autoencoder to discriminate such differences in ENSO modes adds value in the autoencoder representations.
Despite the gross visual similarity of the autoencoder patterns, the finer differences drive unique temporal characteristics (figure 3(a)).Node 3 has a strong significant positive trend with a highfrequency periodicity of 3 years (table 1), whereas node 6 has a strong negative significant trend with a high-frequency periodicity of 3-5 years (table 1).By regressing the time series of figure 3(a) onto global SST, La Niña-like patterns are derived, for the (positive phase of) node 3 (figure 5).The positive trends in node 3 suggest that the western Pacific is warming faster compared to the eastern Pacific during JFM (consistent with the findings of [23]).Trend analysis on global SST data also supports the significant warming over some western parts of the tropical Pacific during JFM (figure S3, top panel).
In general, the encoded JFM patterns have significant periodicities within 3-5 years (table 1), suggesting ENSO variability.Some of the patterns, such as nodes 1, 3, 5, and 9 have statistically significant periodicities in the inter-decadal scale, in their local wavelet spectra, consistent with the IPO (figure 4(a)) [48].Figure 4(a), nodes 5 and 9 have congruence matches of 0.90 with the IPO, supporting the interpretation of the wavelets.Some of the other encoded patterns may be weakly related to other traditional climate modes, with congruence matches above 0.68.For example, during JFM and JJA, node 7 has a congruence match of 0.7 with the PDO.As some climate modes, such as the North Atlantic Oscillation, are defined using a different data set (i.e.sea level pressure) focused on a specific region (i.e. the North Atlantic region), the global SST data used in this study did not capture their variability (figure 4).
We examined the associations between the encoded JFM patterns in figure 2 1).Nodes 1, 2, and 7 have positive significant trends, whereas the remaining nodes have negative significant trends (table 1).Other than the long-term trend, table 1 indicates the periodicity of the analyzed JJA patterns is generally high frequency (between 2 and 7 years, with a few decadal).3 and periodicities in the seasonal mean time series  based on wavelet analysis.± show positive/negative trends.Asterisks ( * ) show the statistical significance of the trend at a 95% confidence level.Dash (-) implies the periodicity is not statistically significant at a 95% confidence level.

Node
Tau-statistic Low-frequency periodicity High-frequency periodicity

Discussion and conclusions
Our work suggests that ANN can serve a crucial role in modeling the non-linear characteristics of the climate system by capturing complex interactions, learning from historical data, and uncovering underlying patterns that may not be easily discernible through traditional linear modeling approaches, making it an important tool for improving climate predictions [26,27,50].In this study, we explore generating non-linear global SST modes using autoencoders.The encoded patterns reveal several unique non-linear SST modes that are not captured by PCA.Consistent with our results, [51] applied convolutional autoencoders to encode temperature field datasets from pre-industrial control runs in the CMIP5 ensemble.The authors found that convolutional autoencoders surpass PCA in nonlinear dimensionality reduction for climate data, offering more accurate reconstructions and potential in building surrogate climate models.Typically, linear modes identified by methods like PCA are based on linear relationships and can process only the linear relationships passed through by a linear correlation matrix.Rotated PCA excels in capturing localized sparse linear patterns where changes in one variable are highly correlated to changes in another.Nonetheless, being a linear method applied to capture non-linear modes, figure S1 (top and bottom panels) shows that rotated PCA merges the different ENSO signals within and between the various El Niños and La Niñas into a single PC mode.Non-linear modes, as identified by techniques like autoencoders, go beyond this by capturing complex, often non-proportional relationships between variables.These relationships can include interactions where the effect of one variable on the other changes depending on the level or context, which linear methods may not detect.Thus, figure S1 indicates that non-linear modes can reveal deeper, more intricate relationships within the data that linear modes might overlook.This capability is particularly valuable in complex systems such as those measured by climate data, where interactions are rarely strictly linear.
The nonlinear analysis allows for a more holistic identification of several spatiotemporal variabilities of the ENSO mode during JFM, as well as their unique global influences on precipitation and temperature.The encoder ANN analyses indicate the western tropical Pacific is warming faster than the ETP (figure S3, top and bottom panels), as evidenced by an asymmetric warming pattern with positive trends in the patterns associated with warmer (cooler) tropical western (eastern) Pacific.The asymmetric SST trend in the tropical Pacific is possibly linked to the cooling influence of aerosols [52].In recent decades (figure S4, right bottom panel), during JJA, the amplified warming of the western tropical Pacific has been linked to Arctic amplification.According to [16], enhanced warming of the tropical western Pacific, relative to the ETP, strengthens the Pacific-Arctic teleconnection, which is a prominent internal mode linking the Arctic to lower latitudes and amplifies Arctic warming and sea ice melting.The Pacific-Arctic teleconnection is characterized by an out-of-phase relationship between tropical central (to eastern) Pacific SSTs and temperature variability over the Arctic (Baxter et al and Sun et al) [16,53].
Studies have shown that during JJA, events like the central Pacific El Niño contribute to summer Arctic cooling [17].Alternatively, [54,55] claim Arctic amplification increases midlatitude stationary waves via the reduction of the equator-to-pole thermal gradient, resulting in the shift of the jet stream.From our results, the localized warming of the western Pacific is more amplified in the summertime (figures S3 and S4), and thus might contribute to why the Pacific-Arctic teleconnection is relatively stronger during summer, consistent with the findings of [16].
Moreover, other non-linear SST modes are shown to have regional temperature and precipitation teleconnections, associated with societal impacts.These include the warming of the Mediterranean and its linkage to increased rainfall over the Sahel, as captured by node 8 during JJA [56].Although we selected to analyze ten encoded patterns out of the 64 encoded representations for the exploratory work herein for brevity, there is the possibility of analyzing additional encoded patterns specific to a study goal.
Our analysis of the ten autoencoder nodes adds to the body of evidence that ANN algorithms have the potential to advance the understanding of the nonlinear climate system and its variability, beyond the knowledge gained from traditional linear approximation techniques.

Figure 1 .
Figure 1.Architecture of the autoencoder model utilized for the seasonal non-linear encoding of the global SST data.

Figure 2 .
Figure 2. Z-score standardized spatial patterns of the first ten encoded patterns of the autoencoder during boreal winter, JFM (a), and boreal summer, JJA (b).The patterns were standardized to ease the interpretation of regions characterized by positive and negative SST anomalies.

Figure 3 .
Figure 3. Annual average time series of the encoded patterns in figure 2 during JFM (a) and JJA (b) from 1979 to 2022.The red dashed line is the 5 year moving average.The time series were derived by projecting the spatial patterns in figure 2 onto monthly global SST anomaly data.

Figure 4 .
Figure 4. Absolute value of the congruence match between the monthly time series of the encoded patterns in figure 2 and climate indices available at https://psl.noaa.gov/data/climateindices/list/during JFM (a) and JJA (b).The values in the map are the congruence coefficients.The analysis period is 1979-2022.
(a) on global temperature and precipitation (figures 6(a) and 7(a)).Given that nodes 1 and 10 (positive phase) are associated with extreme warming of the central to ETP (figure 5(a)), and nodes 5 and 8 (positive phase) are associated with extreme cooling of the ETP during JFM (figure 5(a)), the duo set of patterns has opposing influences on temperature over large landmass regions globally (figure 6(a)), consistent with extreme El Niño (La Niña) events triggering anomalous average global temperature increase (decrease).Compared to temperature, the influence of the SST modes on precipitation is sparser spatially (figure 7(a)), owing to the smaller scale and complexity of its associated processes.Node 1 closely reproduces the Niño 3.4 index (figure 4(a)) and is associated with temperature increase over the northern part of Australia, whereas node 10 more closely reproduces the Niño 4 index (figure 4(a)), affects temperature in Western Australia (figure 6(a)).However, the ocean temperature off Western Australia is warm in node 10 and cold in node 1 (figures 2(a) and 5(a)), illustrating the modulating influence of local SST anomalies.Additionally, node 10 is dry in Western Australia (figure 7(a)), suggesting an association with drought [43].Our results indicate that the different flavors of ENSO detected by the autoencoder can have regionally unique effects on global temperature and precipitation (figures 6(a) and 7(a)) and therefore, are worthy of further exploration.

Figure 5 .
Figure 5. Regression of the seasonal mean time series of the encoded patterns in figures 1 and 2 onto global SST anomaly data from 1979 to 2022 during JFM (a) and JJA (b).Only statistically significant values at a 95% confidence level are plotted.A false discovery rate[49] was applied at a 95% confidence level to control false positives.
Nodes 1 and 7 (positive phase) are associated with enhanced cooling of the central to eastern tropical Pacific (figure 5(b)).
Figure 6(b) shows that the encoded pattern associated with node 1 (positive phase) is linked to the average global temperature increase, with the largest increase occurring in Greenland.Node 7 corresponds with the IPO with a 0.92 congruence match (figure 4(b)) and has a larger association with the temperature of Western Asia.Nodes 1, 4, and 8 are associated with anomalous rainfall in the Sahel, (figure 7(b)).From figure 5(b), nodes 4 and 6 (positive phases) are associated with the warming of the ETP.

Figure 6 (
b) shows that (the positive phase of) nodes 4 and 6 are linked to an average global temperature decrease.Also, the positive phase of node 8 is associated with the Arctic cooling (figure 5(b)), and average global temperature decrease (figure 6(b)).The JJA patterns with a warming signal (from figure 6(b)) have a positive trend (table 1 and figure 3(b)), whereas those with a cooling signal have a negative trend.Generally, unlike in JFM, for the JJA SST modes, the cooling (warming) of eastern tropical Pacific favors an average global temperature increase (decrease).

Figure 6 .
Figure 6.Regression of the seasonal mean time series of the encoded patterns in figures 1 and 2 onto global temperature during JFM (a) and JJA (b).Statistically significant regions at a 95% confidence level are shaded.A false discovery rate[49] was applied at a 95% confidence level to control false positives.

Figure 7 .
Figure 7. as figure 6 but for global precipitation.

Table 1 .
Tau statistics from the modified Mann-Kendall test for the trends in the encoded patterns in figure