Predicting Slowdowns in Decadal Climate Warming Trends With Explainable Neural Networks

The global mean surface temperature (GMST) record exhibits both interannual to multidecadal variability and a long‐term warming trend due to external climate forcing. To explore the predictability of temporary slowdowns in decadal warming, we apply an artificial neural network (ANN) to climate model data from the Community Earth System Model Version 2 Large Ensemble. Here, an ANN is tasked with whether or not there will be a slowdown in the rate of the GMST trend by using maps of ocean heat content (OHC) at the onset. Through a machine learning explainability method, we find the ANN is learning off‐equatorial patterns of anomalous OHC that resemble transitions in the phase of the Interdecadal Pacific Oscillation in order to make slowdown predictions. Finally, we test our ANN on observed historical data, which further reveals how explainable neural networks are useful tools for understanding decadal variability in both climate models and observations.


Introduction
One of the most recognizable indicators of anthropogenic climate change is the positive trend in global mean surface temperature (GMST; Hansen et al., 2010;Johnson et al., 2020). GMST also exhibits interannual to multidecadal variability with periods of accelerations and slowdowns in the rate of decadal trends (Dai et al., 2015;Maher et al., 2020;Thompson et al., 2009;Trenberth et al., 2002). A notable example of one of these GMST slowdowns occurred in the early 2000s (Flato et al., 2013;Fyfe et al., 2013). This temporary warming slowdown ended in the mid-2010s (Mann et al., 2017;Zhang et al., 2019), and more recently, 2020 was one of the three warmest years in the observational record (Dunn et al., 2021). Although the early 2000s was commonly described as a "hiatus" or "pause" in global warming within scientific studies and popular media (Boykoff, 2014;Lewandowsky et al., 2016), we will refer to it here as a "slowdown in decadal warming" (Fyfe et al., 2016), which is more consistent with our understanding of internal variability in the climate system.
Numerous mechanisms have been proposed to explain the cause of the early 2000s slowdown, as reviewed in Medhaug et al. (2017) and Xie and Kosaka (2017), but it was likely a combination of factors ranging from uncertainties in the observational data record (e.g., Cowtan & Way, 2014;Karl et al., 2015), fluctuations in radiative forcing (Schmidt et al., 2014), cooling in the eastern Pacific associated with a negative phase of the Interdecadal Pacific Oscillation (IPO; England et al., 2014;Meehl, Hu, et al., 2013;Roberts et al., 2015), anthropogenic aerosol and volcanic forcing (Santer et al., 2014;Smith et al., 2016), changes in deep ocean heat uptake (Watanabe et al., 2013), top-of-atmosphere (TOA) energy imbalance (Hedemann et al., 2017;Meehl et al., 2011), and interactions between modes of climate variability (W. Liu & Xie, 2018). Motivated by the increasing body of Abstract The global mean surface temperature (GMST) record exhibits both interannual to multidecadal variability and a long-term warming trend due to external climate forcing. To explore the predictability of temporary slowdowns in decadal warming, we apply an artificial neural network (ANN) to climate model data from the Community Earth System Model Version 2 Large Ensemble. Here, an ANN is tasked with whether or not there will be a slowdown in the rate of the GMST trend by using maps of ocean heat content (OHC) at the onset. Through a machine learning explainability method, we find the ANN is learning off-equatorial patterns of anomalous OHC that resemble transitions in the phase of the Interdecadal Pacific Oscillation in order to make slowdown predictions. Finally, we test our ANN on observed historical data, which further reveals how explainable neural networks are useful tools for understanding decadal variability in both climate models and observations. Plain Language Summary Long-term observations reveal that Earth's average temperature is rising due to human-caused climate change. Along with this warming trend are also variations from year-toyear and even over multiple decades. This temperature variability is often tied to regional patterns of heat in the deep ocean, which can then modulate weather and climate extremes over land. In an attempt to better predict temperature variability on decadal timescales, we use a machine learning method called artificial neural networks and data from a climate model experiment, which was designed to compare climate change and variability. Here, our artificial neural network (ANN) uses maps of ocean heat to predict the onset of temporary slowdowns in the rate of global warming in both the climate model and in real-world observations. We then use a visualization technique to find which areas of ocean heat that the ANN is using to make its correct predictions, which are found to be mainly across the Pacific Ocean. In agreement with recent research, our study finds that new data science methods, like machine learning, can be useful tools for predicting variations in global climate. literature on the causes and impacts of the early 2000s slowdown, we aim to investigate the predictability of similar temporary GMST slowdowns occurring in a warming climate. While decadal predictability has been explored using other statistical methods (e.g., Mann et al., 2016;Sévellec & Drijfhout, 2018), sensitivity experiments (e.g., Kosaka & Xie, 2013), and hindcasts with initialized state climate modeling frameworks (e.g., Boer et al., 2016;Fyfe et al., 2011;Guemas et al., 2013;, we explore this problem through the lens of a machine learning prediction task. Deep learning methods, such as neural networks, have the ability to extract and leverage nonlinear patterns across data-intensive spatial fields, which make them promising tools for revealing new insights and sources of predictability in climate science (Barnes, Mayer, et al., 2020;Irrgang et al., 2021;Reichstein et al., 2019;. Recent work has demonstrated the utility for neural networks in identifying climate modes, teleconnections, and forecasts of opportunity for a wide variety of timescales (e.g., Gibson et al., 2021;Gordon et al., 2021;Ham et al., 2019;J. Liu et al., 2021;Mayer & Barnes, 2021;Nadiga, 2021;Tang & Duan, 2021;Toms et al., 2021;Wu & Hsieh, 2004). Further, a growing number of explainable artificial intelligence (XAI) methods have been adapted for applications in weather and climate science (McGovern et al., 2019;Toms et al., 2020), which can retrospectively trace the decisions of neural networks and assist scientists in comparing the attribution of input features to known physical mechanisms in the Earth system. Besides evaluating trust and credibility to the machine learning prediction, XAI methods can also be used for physics-guided scientific discovery and hypothesis testing (Ebert-Uphoff & Hilburn, 2020;Toms et al., 2020).
In this study, we use an artificial neural network (ANN) to explore the predictability of decadal warming slowdowns due to variability in the upper ocean within a new large ensemble experiment and real-world observations. In addition to slowdown predictability, we also use a complimentary XAI method to investigate the oceanic patterns that may provide insight to these temporary warming slowdowns.

Climate Model Large Ensemble
For climate model data, we use a large ensemble experiment conducted by the Community Earth System Model Version 2 (CESM2; Danabasoglu et al., 2020; see Supporting Information for more details). Specifically, we use simulations from the CESM2 Large Ensemble Community Project (CESM2-LE; Rodgers et al., 2021), which includes 100 ensemble members branched from the fully coupled CESM2 preindustrial control (1850 radiative forcing conditions) using different atmospheric and oceanic initial states. CESM2-LE members follow historical Coupled Model Intercomparison Project Phase 6 (CMIP6) forcing from 1850 to 2014 and thereafter follow the SSP3-7.0 future radiative forcing (high emissions scenario) until 2100 O'Neill et al., 2016). We consider the first 50 ensemble members (1-50), which are prescribed with biomass burning emissions following CMIP6 protocol (Van Marle et al., 2017). In contrast, the second set of 50 ensemble members follow temporally smoothed biomass burning fluxes (51-100). As discussed in Fasullo et al. (2022) and Rodgers et al. (2021), this difference in biomass burning forcing has been shown to affect large-scale climate features, including the GMST record in present day.
Due to limited data availability at the time of our analysis, we analyze only 40 ensemble members within the first subset of CESM2-LE (1-50). From these 40 members, we use monthly outputs of near-surface air temperature (T2M) and sea surface temperature (SST). We also utilize monthly ocean heat content (OHC), which is derived as the vertical heat content integral between three distinct depth layers (0-100 m, OHC100; 0-300 m, OHC300; 0-700 m, OHC700); although we focus on maps of OHC100 for the actual training of our ANN. We then apply a bilinear interpolation to all variables so that they share a common (slightly coarser) latitude by longitude grid (1.9° × 2.5°). We calculate annual means from the monthly data and use the period from 1990 to 2099 to classify slowdowns in decadal warming. To focus on warming slowdowns driven by internal variability, we remove the 40-member ensemble mean from each individual ensemble in every year and grid box for SST and OHC (Maher et al., 2021;Phillips et al., 2020).

Observations
To evaluate our ANN trained on CESM2-LE for predicting the early 2000s warming slowdown in the historical record, we use SST and T2M from the European Centre for Medium-Range Weather Forecasts ERA5 reanalysis (Hersbach et al., 2020) and OHC from the Institute of Atmospheric Physics (IAP) ocean gridded product (Cheng et al., 2017;Cheng & Zhu, 2016; both data sets referred to here as "observations"). SSTs from ERA5 are an interpolated product between HadISST2 (Titchner & Rayner, 2014) from January 1979 to August 2007 and OSTIA (Donlon et al., 2012) from September 2007 to present. Overall, both regional and global mean time series of SST and T2M are consistent with other observational data sets (Bell et al., 2021;Hersbach et al., 2020), including decadal trends ( Figure S1 in Supporting Information S1). Gridded upper OHC from IAP also compares well with in situ measurements and is based on temperature data from the World Ocean Database (Boyer et al., 2013), which is then further bias-corrected, interpolated, and quality controlled (Cheng et al., 2017;Li-Jing et al., 2015).
In all observations, we use monthly output and bilinearly interpolate these fields onto the same 1.9° × 2.5° grid as CESM2-LE before calculating annual means. We linearly detrend each grid point (over 1979-2020) for SST and OHC predictors to remove long-term warming signals and thus focus on patterns of interannual variability for slowdown predictions. Figure 1 shows an example of how we define warming slowdown events in CESM2-LE and observations. While there have been numerous definitions and data sets used for identifying warming slowdown (or so-called hiatus/ pause) events (e.g., Risbey et al., 2018;Wei et al., 2021), they are often classified as a near-zero or negative 10-20 yr linear trend of the GMST. Recent studies show that the frequency of slowdown events in CMIP5 models decreases substantially by the end of the 21st century using a negative 10 yr linear trend definition (e.g., Li & Baker, 2016;Maher et al., 2014;Sévellec et al., 2016). Their frequency could also decrease due to increasing climate sensitivity (Modak & Mauritsen, 2021). Yet, internal variability is still projected to affect regional and global climate trends even under higher future emission scenarios (Cassou et al., 2018;Easterling & Wehner, 2009;Li & Baker, 2016;Maher et al., 2020). Here, we take a GMST slowdown threshold using decadal (10 yr) trends (e.g., Maher et al., 2014;Meehl et al., 2011), which is on the shorter end of previously used timescales and focus on the influence of oceanic internal variability.

Defining Slowdowns in Decadal Warming
First, to classify slowdown events in observations, we compute the area-weighted GMST and calculate 10 yr moving linear trends beginning in 1990. We start our analysis in 1990 to avoid any multidecadal slowdown events earlier in the 20th century when the influence of the forced climate change signal may not have fully emerged (Delworth & Knutson, 2000;Hawkins et al., 2020;Papalexiou et al., 2020). We then calculate the mean of all 10 yr trends between 1990 and 2020 and take one standard deviation below this mean as our threshold for slowdown events in observations (equating to about +0.01°C/yr, or 0.44 of the mean observational trends; black dashed line in Figure 1b). We identify four consecutive slowdown events in observations, which begin in 2002. These years are consistent with previous studies  and are similarly classified in other data sets with our definition ( Figure S1 in Supporting Information S1).
For CESM2-LE, we first compute the area-weighted GMST for the ensemble mean from all 40 members through 2099 (Figure 1a). We then calculate 10 yr moving linear trends, which begin in 1990 for consistency with observations. The climate model threshold for a slowdown (red dashed line in Figure 1b) is defined by multiplying the fraction of the mean trend from the observations (0.44) times each trend period from the ensemble mean (thick gray line in Figure 1b).
Then, we calculate the GMST for each ensemble member and the associated moving 10 yr linear trends (thin gray lines in Figure 1b). We define a warming slowdown event when these 10 yr trends fall below the climate model threshold. That is, we define a slowdown event as a fraction of the forced response. However, our definition still leads to a reduction in the number of slowdown events after 2040 in CESM2-LE (Figure 1c), as shown in several past studies (e.g., Sévellec et al., 2016).

Artificial Neural Network
For this analysis, we adopt a neural network architecture that is designed to receive input maps of OHC100 anomalies and output whether the next 10 yr will observe a decadal warming slowdown. In other words, if we input a map of OHC100 anomalies for the year 2000, the ANN will output whether the decade from 2000 to 2009 will be a slowdown event or not. A schematic of our ANN can be found in Figure S3, and the architecture parameters are outlined in the Supporting Information S1.
In addition to seeing if warming slowdown events are predictable, we are also interested in the sources of predictability in fields of anomalous OHC100. To attempt to understand the ANN's decision-making process, we use a method of XAI called layer-wise relevance propagation (LRP;Bach et al., 2015;Montavon et al., 2017Montavon et al., , 2018. The utility of LRP has been demonstrated in a wide range of weather and climate applications (e.g., Davenport & Diffenbaugh, 2021;Gordon et al., 2021;Labe & Barnes, 2021;, and an overview for the geosciences can be found in Toms et al. (2020). In short, prior to the softmax, a single prediction output is propagated backward through the ANN after freezing the model weights and biases. LRP then returns a vectorized spatial map, which shows the feature relevance for every input sample's latitude and longitude pixel. Therefore, we have a unique LRP heatmap for every input sample of OHC100. Throughout this study, regions of higher relevance can be interpreted as more important for the ANN's prediction. We implement the LRP z rule for back propagation, which was found by Mamalakis et al. (2021) to be a well performing XAI method using a benchmark climate data set similar to ours. To improve interpretation and reduce the amount of noise in the LRP heatmaps, we only focus on positive areas of relevance, which are features that contribute positively to the ANN's prediction output. Figure 2a shows the results of our ANN for each CESM2-LE ensemble member in the testing data set from 1990 to 2090 (i.e., 2090-2099 is the last complete decade of data). Given the large class imbalance, we focus on the F1 score (balancing precision and recall), rather than categorical accuracy, to evaluate the performance of our ANN for correctly identifying slowdown events. Figure S7 in Supporting Information S1 provides a collection of skill metrics for our testing data (Accuracy = 0.87, Precision = 0.39, Recall = 0.41). Overall, the network achieves an F1 score of 40% and performs better than random chance (10.4%). While our ANN sometimes struggles with correctly classifying slowdown events, especially those that occur simultaneously in a row, it generally classifies at least one 10 yr period during these extended events. This skill suggests that the ANN is learning information from OHC100 anomalies that corresponds to future slowdown periods in CESM2-LE.

Predicting Slowdown Trends in a Large Ensemble
We test the robustness of our results by training 100 ANNs with unique random initialization seeds and different combinations of ensemble members used for training, validation, and testing data. The F1 score of our single seed ANN falls around the ≈85th percentile of this distribution, and additional metric scores are shown in Figure  S8 in Supporting Information S1 for the 100 ANNs. The spread between this distribution can be attributed to uncertainties from random initialization states of the ANNs and different combinations of ensemble members. This suggests that differences in the skill of slowdown predictions across the ANN distribution could be related to individual realizations of internal variability as simulated per each ensemble member. Although we found that there is no relationship between the accuracy of testing predications compared to the number of training slowdown events each ANN learned for the 100 iterations ( Figure S9 in Supporting Information S1).

Sources of Predictability for Slowdowns
To understand the sources of skill for the ANN's correct slowdown predictions in CESM2-LE, we turn to composite maps of LRP. Recall that LRP traces the decision-making process of a neural network, where higher relevance corresponds to greater importance for the ANN to make its final prediction. While we have LRP heatmaps for every input of annual-mean OHC100, we focus on correct predictions by the ANN in the testing data set. Figure 2b shows the LRP composite for all correct slowdown predictions. We find higher relevance in the off-equatorial regions of the eastern Pacific, especially in the regions of the North/South Pacific Meridional Modes (Amaya, 2019). There are also patches of higher relevance across portions of the Indian Ocean, South Atlantic, and South Pacific, which suggests that the ANN is leveraging other regional patterns of OHC anomalies to make predictions. Notably, there is no relevance for a thin band along the equator in the area of the El Niño-Southern Oscillation. Figure 2d shows the corresponding LRP composite for no slowdowns in decadal warming, which shows a similar spatial pattern of relevance across the equatorial Pacific as in Figure 2b, but higher relevance in other ocean basins. Likewise, we also find comparable LRP composites for the incorrect slowdown predictions, although there are notable differences over portions of the Southern Ocean, lower Arctic, and near New Zealand (Figure S10b in Supporting Information S1). These XAI results are a product of the setup of our binary classification problem, and therefore the LRP maps reveal the regions that the ANN is using to make this determination (i.e., yes or no slowdown), but these patterns may not always necessarily correspond to an actual slowdown events driven by OHC variability.
We compare these LRP maps to composites of the raw (normalized) OHC anomalies that were input to the network for correct slowdown predictions (Figure 2c) vs. correct no slowdown predictions (Figure 2e). Now we find striking differences between the two OHC patterns. The composite of OHC100 for the slowdown predictions reveal an IPO-like spatial pattern with cold pools in the west-central North Pacific and west-central South Pacific and warm anomalies in the Southern Ocean and eastern Pacific. We also see the signature of a positive Indian Ocean Dipole (IOD; Saji et al., 1999) and a dipole pattern of OHC100 anomalies between the southern Atlantic and north-central Atlantic. Some studies have shown that a positive IOD can be a precursor for a rapid transition to a cooler equatorial Pacific by modulating the strength of the Walker Circulation (Izumo et al., 2010;Le et al., 2020;Yoo et al., 2020). Figure S11 in Supporting Information S1 shows maps of OHC anomalies at other vertical depth levels for the slowdown predictions compared to 5-10 yr after the start of the slowdown decade. We find a similar spatial pattern of SSTs ( Figure S11a in Supporting Information S1), but a stronger cold pool at deeper depths, which appears to be propagating eastward in the equatorial Pacific (Figures S11c-S11d in Supporting Information S1). In contrast, we find a negative IPO-like pattern for the composites at the end of the slowdown decade (Figures S11e-S11h in Supporting Information S1). This finding is in agreement with earlier studies that showed slowdowns in decadal warming often correspond to trends toward a negative phase of the IPO within CMIP5 models (e.g., Maher et al., 2014). Given this evolution of events and the patterns of LRP relevance, it is feasible that the ANN is learning OHC anomalies associated with transitions in the state of the IPO.
To directly assess the IPO in the maps of OHC100, we compute the unfiltered IPO Tripole Index (normalized) following Henley et al. (2015) using annual-mean SSTs from CESM2-LE ( Figure S12 in Supporting Information S1). As expected from the composite analysis in Figure 2, we find that correct predictions of slowdowns generally correspond to highly positive phases of the IPO index. To demonstrate this point, we select one ensemble member and compare its annual IPO index to the frequency of the slowdowns classifications in the distribution of 100 unique ANNs ( Figure S13 in Supporting Information S1). We find slowdown predictions often correspond to a positive IPO index in this ensemble member, but also importantly, not every positive IPO results in the prediction of a slowdown event.
To further confirm that the ANN is learning additional spatial information beyond a simple reflection of the IPO-like pattern of OHC anomalies, we set up a logistic regression problem by inputting only the value of the IPO index in CESM2-LE to predict whether a slowdown event will occur over the next 10 yr (Accuracy = 0.75, Precision = 0.2, Recall = 0.46, F1 score = 0.28). Thus, we find that using global maps of OHC100 as inputs to the fully connected ANN provides more skillful predictions of warming slowdown events.

Predicting Slowdown Trends in Observations
Lastly, we test the utility of our neural network for capturing the observed early 2000s slowdown by inputting maps of OHC100 from observations, which are first linearly detrended and then normalized by their own mean and standard deviation at every grid point. Figure 3a shows the slowdown prediction from our ANN for each input map of annual-mean OHC100 using observations. During the overlapping period with the actual early 2000s slowdown, the decades from 2003 to 2012 and 2004 to 2013 are classified as warming slowdowns. The range in observational predictions across the distribution of 100 ANNs is also shown in Figure S14 in Supporting Information S1.
To understand the patterns of anomalous OHC100 that the ANN is using to make its prediction for observations, we evaluate a LRP composite map from the single seed ANN (which correctly predicted two slowdown events in the early 2000s) in Figure S15a in Supporting Information S1. Similar to the LRP composites using CESM2-LE (Figure 2), we find areas of higher relevance in the equatorial western Pacific, south-central Atlantic, and patches in the Indian Ocean. The ANN also predicts the onset of slowdown events mainly during positive phases of the IPO (Figure 3b). Although this correlation is not always the case (e.g., during the positive IPO event in the early 1990s), which again suggests that our ANN is leveraging additional spatial information than 10.1029/2022GL098173 8 of 14 simply the IPO pattern to make predictions. This is also supported by our interpretation of the LRP maps, which show higher relevance regions across the western Pacific and not necessarily the canonical IPO/PDO patterns (Newman et al., 2016;Parker et al., 2007; Figure S15a in Supporting Information S1).
At the time of our analysis, the last complete decade of GMST observations covers the decade of 2011-2020 ( Figure 1a). However, since we only need OHC prior to predicting the future 10 yr, we can also explore warming slowdown events extending beyond 2020 (Figure 3a). For these future predictions, 2016-2025 and 2017-2026 are classified as warming slowdown events by the ANN. 2016 was characterized by the dissipation of an extreme El Niño event into a weak La Niña state (Santoso et al., 2017), and the GMST also set a new record high for that respective year (Aaron-Morrison et al., 2017). Similarly, the IPO index also shows a transition from a highly positive phase in 2015 to a neutral or negative phase in the following years through 2020 (Figure 3b). Composites of normalized SST and OHC for 2016 and 2017 show anomalously warm subsurface waters just off the equator in the eastern Pacific and cold pools in the tropical Indo-Pacific and north-central Pacific ( Figure S16 in Supporting Information S1). Comparing the LRP composite map over 2016 and 2017 with the raw OHC100 anomalies (Figures S15b and S16b in Supporting Information S1), we find higher relevance outlining the warm anomalies in the eastern Pacific and patches of relevance in the Indian Ocean and southern Pacific. The LRP composite for the future slowdown prediction in Figure S15b in Supporting Information S1 is more similar to those outlined in

Summary and Conclusions
Due to an increasing need from decision makers and other community stakeholders for near-term climate predictions, there has been a coordinated effort to increase the availability of decadal outlooks from initialized climate models and other operational forecast systems (e.g., Boer et al., 2016;Graham et al., 2011;Hewitt et al., 2021;Kushnir et al., 2019;Meehl, Goddard, et al., 2013;Meehl et al., 2021;Merryfield et al., 2020;Smith et al., 2019). However, these simulations can be computationally expensive to run. Alternatively, recent progress in machine learning has shown promising results for decadal climate applications, especially when combined with explainability methods (e.g., Gordon et al., 2021;G. Liu et al., 2021;Toms et al., 2021). Motivated by this new line of research, we explore the utility of a relatively shallow ANN for predicting temporary decadal warming slowdowns of GMST using upper OHC variability. Although our ANN is trained on climate model data from CESM2-LE, we find that it also produces skillful predictions of the early 2000s warming slowdown in observational data. We further compliment our ANN with a machine learning explainability method (LRP) to attempt to understand what information the ANN is using to make its correct predictions. The LRP maps reveal that the ANN is mainly using off-equatorial anomalies of OHC100 to predict the onset of a decadal warming slowdown. These patterns suggest that the ANN may be learning precursors for transitions to a negative phase of the IPO, although this topic remains an active area of research (Cai et al., 2019;Power et al., 2021).
Finally, we note a few important considerations when interpreting these results. First, the causal mechanisms related to the early 2000s slowdown event remain uncertain (Hedemann et al., 2017;Medhaug et al., 2017;von Känel et al., 2017), and we note that we have only considered one potential predictor for warming slowdowns (i.e., upper OHC). Slowdowns can also occur due to external forcing (e.g., aerosols) or other modes of climate variability (Medhaug et al., 2017). Further, there are a number of different definitions for slowdown-like events . Future work could explore the predictability of slowdowns using ANNs with other climate predictors, such as considering a TOA energy imbalance approach (Hedemann et al., 2017), or taking into account longer duration events. It may also be valuable to combine maps of OHC at different lead times, which was recently demonstrated by Gordon et al. (2021) for predicting transitions in the phase of the PDO. Lastly, we train our ANN on a large ensemble from only one climate model. Thus, our results may be influenced CESM2's inherent model biases and prescribed external forcing, including the SSP3-7.0 emissions scenario and the protocol for biomass burning (Rodgers et al., 2021). The value of using multi-model large ensembles and adding more complexity to the neural network, such as designing a convolutional neural network to evaluate regional OHC patterns, will be left for future exploration. Importantly, even our simple ANN demonstrates that temporary warming slowdowns may have some predictability from Pacific climate variability and demonstrates a new application of machine learning for climate science.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.