Southern Ocean Dynamics Under Climate Change: New Knowledge Through Physics-Guided
Machine Learning

Complex ocean systems such as the Antarctic Circumpolar Current (ACC) play a central role in the climate through mechanisms such as heat transport via powerful currents, large-scale overturning, and carbon exchange with the atmosphere [1, 2, 3]. Such systems are known to exhibit a wide array of responses to anthropogenic forcing, but such changes and their underlying physics are poorly constrained [3, 4, 5]. This is, in part, a direct consequence of difficulties in both characterizing ocean physics and tracking its shifts under climate change. Model intercomparison projects, such as the Coupled Model Intercomparison Project (CMIP6) [6, 7, 8], have set up a vital framework for understanding differences between climate models in both of these areas. A key target for the sixth and current phase of CMIP was to increase the horizontal resolution, which in the ocean partially allows for mesoscale turbulence to be expressed. However, modern models have only exacerbated challenges in characterizing and tracking ocean physics due to their complexity and size. As a result, key systems such as the ACC as well as their variability among models are often described using bulk metrics, such as, the overturning strength integrated over global longitude [1, 9]. Bulk metrics have limited scope for elucidating the dynamics driving the changes, as these largely vary in space. For the first time, we apply the method Tracking global Heating with Ocean Regimes (THOR) proposed by [10] to the ocean of a mesoscale turbulence-permitting model, the Coupled Model 4 [11, 12], to both understand and track ocean dynamics under climate change. THOR is a three-step process beginning with unsupervised classification of ocean grid cells into dynamical regimes based on their physics. For this step, we use Native Emergent Manifold Interrogation (NEMI) [13]. The second component of THOR trains a deep ensemble of neural networks to predict these dynamical regimes from the previous step using readily accessible fields as inputs (e.g., sea surface height, depth, wind stress, and mass transport). Finally, the third step of THOR applies this trained deep ensemble to a new climate change scenario or entirely new model of interest to understand changes in dynamical regimes over time. This step also incorporates uncertainty quantification and eXplainable Artificial Intelligence (XAI) methods to guide further exploration of newly discovered ocean phenomena. Here we use THOR to investigate changes in the the region where the ACC meets the Pacific-Antarctic Ridge (PAR). THOR specifically reveals a shift in dynamical regime in this region under climate change, and XAI methods additionally show that these are related to changes in wind stress and interactions with the bathymetry. Using this new information provided by THOR, we further explore the ACC under climate change and find that the shifts in regional dynamical regime is explained by northward movement of the ACC driven by changes in wind stress. This ACC movement brings it into a new, less variable bathymetric region where its interactions with the sea floor are less strong, thus leading to a flow distributed over a larger area and concentrated at the surface. The remainder of this manuscript is structured as follows. Section 2 provides a more detailed overview of the THOR method. Section 3 details our exploration of the ACC guided by new knowledge revealed by THOR. Lastly, Section 4 concludes the paper with potential directions for future work.

In this manuscript we apply THOR (Figure 2) to the Modular Ocean Model version 6 (MOM6) [14, 12, 15], a component of the Coupled Model version 4 (CM4) [11]. CM4 operates at an ocean resolution of $0.25^{\circ}$ which allows for mesoscale eddies, differentiating itself from the $1^{\circ}$ ECCO model [16] studied by [10]. The first step of THOR addresses the key issue of unlabeled data in ML applications for climate sciences by leveraging unsupervised learning to create a labeled dataset from high dimensional ocean model data. Specifically, ocean grid cells are clustered into groups, called dynamical regimes, based on their mean balance of the barotropic vorticity (BV) equation, a characterization of local circulation of fluid flow, throughout a preindustrial control (piControl) run, illustrating the state of the ocean before the start of anthropogenically induced climate change. For more details, see [10, 17].

Previous applications of THOR to the ECCO model [10, 18] used k-means clustering to reveal dynamical regimes. However, [13] demonstrates that k-means fails to converge for several information criteria metrics on the higher resolution MOM6 data. As such, for our first step of THOR (Figure 2a) we substitute k-means for Native Emergent Manifold Interrogation (NEMI), which has been shown by [13] to successfully cluster MOM6 BV data into a user-defined number of interpretable dynamical regimes. Here, we choose six labels, or dynamical regimes, following [10]. The global regimes found by NEMI for the piControl run of CM4 are shown in Figure 1. Notice the oceanographic features which are reflected in the dynamical regimes such as large, wind stress-driven gyres and the the mid-Atlantic ridge. Additionally, note that the vast majority of grid cells are clustered into Regimes 0-4. This indicating that these capture the dominant dynamical regimes in the ocean [13] and that Regime 5 represents residual nonlinear regimes, as in [19].

While the first step of THOR is useful for revealing areas of the ocean characterized by similar physics, it requires that the BV budget be closed for each new ocean model and experiment in order for NEMI clustering to be done. This is infeasible for a wide array of models and impossible with the data routinely made available through CMIP. As such, the second step of THOR (Figure 2b) trains a neural network (NN) to predict the clustered dynamical regimes from above from more accessible fields from an ocean model. Specifically, the inputs for each grid cell are the sea surface height above the geoid (ZOS), its $x$- and $y$-gradients, the depth relative to sea level (bathymetry), its $x$- and $y$-gradients, the curl of surface wind stress torque ($\nabla\times\tau_{s}$), the Coriolis parameter ($f$), and the depth-summed zonal and meridional mass transport (umo_2d and vmo_2d). The notable additions to this dataset (differentiating it from that of [10]) are the mass transport fields. These were added after preliminary testing showed that neural networks were unable to accurately classify pairs of regimes whose dominant physics were correlated with mass transport. Each of these inputs was chosen because they directly impact the values of BV equation terms, meaning that there is some physical relation between the chosen inputs and the dynamical regimes found in the first step of THOR. Thus, eXplainable artificial intelligence (XAI) methods may be applied to attribute NN predictions to specific physical inputs, allowing us to gain insight into the drivers of regime shift as demonstrated in the next section. During initial testing XAI also revealed spurious correlations between the input fields and output regimes which guided our choice to add the mass transport variables to the input.

We use the same MLP architecture described in [10] to predict the dynamical regimes, and an ensemble of 50 such MLPs are trained in order to regularize the latent loss space and gauge uncertainty in regime predictions. Specifically, we use entropy across the ensemble’s predictions for each grid cell to quantify uncertainty following [20]. More detailed results from training can be found in Section A.1. Importantly, we find that though the NN does not achieve exceptionally high accuracy, the average entropy of incorrectly classified grid cells is higher than that of correctly classified grid cells. This indicates that the NN is making well-calibrated predictions, correctly represents its own uncertainty, and is not erroneously overconfident in its predictions.

The third and final step of THOR (Figure 2c) applies the trained NN ensemble from the previous section to predict dynamical regimes for entirely new experiments or ocean models. Here, we apply the ensemble to the Historical (post industrial revolution to 2015) and SSP585 (aggressive anthropogenic warming post 2015) runs of CM4 in order to track regime shifts in the Southern Ocean under climate change. We analyze the 1992-2013 window of the Historical scenario following [10] and the 2080-2100 window of the SSP585 scenario to elucidate the effects of anthropogenic warming. Figure 3 shows the deep ensemble’s mean Southern Ocean regime predictions for the Historical and SSP585 scenarios as well as the entropy for the SSP585 predictions as a proxy for the ensemble’s uncertainty in its predictions under climate change. Acting as a tool to guide new discovery, THOR reveals a regime change where the ACC meets the Pacific-Antarctic Ridge (PAR), a divergent tectonic plate boundary characterized by rough bathymetry at around $60^{\circ}$S, $166^{\circ}$W. This region of interest is shown in the insets of Figure 3. The notable shift revealed by THOR from Regime 4 (light green) to Regime 0 (dark blue), combined with the high entropy motivates further investigation.

We use two XAI methods, the NN specific layer-wise relevance propagation (LRP) [21] and model agnostic Shapley Additive exPlanations (SHAP) [22], to determine which inputs help THOR’s deep ensemble make its predictions where the ACC meets the PAR. The two values generally agreed on the relevance of the surface wind stress curl and the bathymetry, as well the sea surface height which is closely related to the wind (see Section A.2). Given that the ACC is largely driven by westerly winds surrounding Antarctica, the XAI methods may be indicating a shift in the wind stress and the interactions between the ACC and the bathymetry. Motivated by this, we investigate lateral (east-west) mass transport at $166^{\circ}$W, which represents the portion of the ACC which cuts through the PAR, as shown in the bottom two panels Figure 4. Here, we see a clear shift in the core of the ACC northward away from the PAR. Specifically, in the Historical scenario the eastward jets are closely connected to the bathymetric features of the PAR and concentrated south of diffuse circulation without a jet-like structure (intense red) at $53^{\circ}$S. In the SSP585 scenario, however, the area of intense flow around $50^{\circ}$S strengthens markedly and is associated with a strengthening of the wind stress. As such, the core of the ACC moves into a region of less dynamic bathymetry, allowing it to flow more freely and increase in intensity. Aligning with the new knowledge of wind stress relevance, the local wind stress curl strengthens just north of the PAR as shown in the top panel of Figure 4. In response, the ACC shifts in the same direction, pulling it away from the bathymetric influence of the PAR and fundamentally changing its governing physics, as reflected in the regime shift initially revealed by THOR in Figure 3.

We extend the machine learning method THOR to a mesoscale eddy-permitting climate model in order to gain insight into the drivers underlying changes in ocean physics under climate change. Importantly, we use THOR not as an oracle, but rather a tool to discover new knowledge about ocean dynamical regimes and guide further investigation of shifts in the driving forces behind ocean circulation. Under this framework, THOR reveals that a fundamental shift in Southern Ocean physics occurs under climate change where the ACC meets the PAR. Using this new knowledge combined with XAI methods as a guide, we find that the wind stress curl increases in strength north of the PAR, pulling the ACC with it and effectively moving the ACC from a bathymetrically locked state over the PAR to one dominated by the wind stress where its flow increases in strength. It is important to note that while such conclusions may theoretically be reached by analyzing the raw ocean model data, such analysis is largely infeasible due to the high resolution and dimensionality of the model output. THOR addresses these key issues and enables us to easily identify areas of interest for further investigation such as those discussed in Section 3.

There are several possible directions for future work. First, different CMIP models have different representations of ocean physics, and the high-resolution version of THOR presented in this work could provide insight into how these differences affect predictions for various ocean systems under climate change. Such insight is at the heart of lowering model spread and uncertainty of future projections. Additionally, since NEMI can cluster the MOM6 data into any given number of regimes, it is possible to train the NNs of THOR to identify more than 6 dynamical regimes discussed in Section 2. Future work will investigate additional oceanographic insights revealed by THOR when analyzing a greater number of regimes. Lastly, the NN ensemble applied in this work classifies the dynamical regimes of grid cells based on local information to that grid cell. Spatially-aware NN architectures may be applied, but the large continents and difficulties applying XAI methods to more complex architectures pose exciting future challenges. Regardless of the supervised technique used to classify dynamical regimes, we maintain that its predictive power should not trade off with its interpretability and transparency.

The authors would like to thank Stephen Griffies, Tarun Verma, and Isaac Held for their helpful discussions and feedback. WY acknowledges funding from the NOAA Ernest F. Hollings Undergraduate Scholarship. MS funding: Cooperative Institute for Modeling the Earth System, Princeton University, under Award NA18OAR4320123 from the National Oceanic and Atmospheric Administration, U.S. Department of Commerce.

The data and code are publicly available at https://github.com/yikwill/THOR-MOM6.