Deep Learning of Sea Surface Temperature Patterns to Identify Ocean Extremes

We perform an out-of-distribution analysis of ~12,000,000 semi-independent 128x128 pixel^2 SST regions, which we define as cutouts, from all nighttime granules in the MODIS R2019 Level-2 public dataset to discover the most complex or extreme phenomena at the ocean surface. Our algorithm (Ulmo) is a probabilistic autoencoder, which combines two deep learning modules: (1) an autoencoder, trained on ~150,000 random cutouts from 2010, to represent any input cutout with a 512-dimensional latent vector akin to a (non-linear) EOF analysis; and (2) a normalizing flow, which maps the autoencoder's latent space distribution onto an isotropic Gaussian manifold. From the latter, we calculate a log-likelihood value for each cutout and define outlier cutouts to be those in the lowest 0.1% of the distribution. These exhibit large gradients and patterns characteristic of a highly dynamic ocean surface, and many are located within larger complexes whose unique dynamics warrant future analysis. Without guidance, Ulmo consistently locates the outliers where the major western boundary currents separate from the continental margin. Buoyed by these results, we begin the process of exploring the fundamental patterns learned by Ulmo, identifying several compelling examples. Future work may find that algorithms like Ulmo hold significant potential/promise to learn and derive other, not-yet-identified behaviors in the ocean from the many archives of satellite-derived SST fields. As important, we see no impediment to applying them to other large, remote-sensing datasets for ocean science (e.g., sea surface height, ocean color).


INTRODUCTION
Satellite-borne sensors have for many years, been collecting data used to estimate a broad range of meteorological, oceanographic, terrestrial and cryospheric properties. Of significance with regard to the fields associated with these properties is their global coverage and relatively high spatial (meters to tens of kilometers) and temporal (hours to tens of days) resolutions. These datasets tend to be very large, well documented and readily accessible making them ideal targets for analyses using modern machine learning techniques. Based on our knowledge of, interest in and access to global sea surface temperature (SST) datasets, we have chosen one of these to explore the possibilities. Specifically, inspired by the question of "what lurks within" and also the desire to identify complex and/or extreme phenomena of the upper ocean, we have developed an unsupervised machine learning algorithm named ulmo 1 to analyze the nighttime MODerate-resolution Imaging Spectroradiometer (MODIS) Level-2 (L2) 2 SST dataset obtained from the NASA spacecraft, Aqua, spanning years 2003-2019. The former (the unknown unknowns) could reveal previously unanticipated physical processes at or near the ocean's surface. Such surprises are, by definition, rare and require massive datasets and semi-automated approaches to examine them. The latter type (extrema) affords an exploration of the incidence and spatial distribution of complex phenomena across the entire ocean. Similar 'fishing' expeditions have been performed in other fields on large imaging datasets (e.g., astronomy Abul Hayat et al. 2020). However, to our knowledge, this is the first application of machine learning for open-ended exploration of a large oceanographic dataset, although there is a rapidly growing body of literature on applying machine learning techniques to the specifics of SST retrieval algorithms Saux Picart et al. (2018), cloud detection Paul & Huntemann (2020), eddy location Moschos et al. (2020), prediction Ratnam et al. (2020); Zhang et al. (2020); Yu et al. (2020), etc. and, more generally, to remote sensing Ma et al. (2019).
Previous analyses of SST on local or global scales have emphasized standard statistics (e.g., mean and RMS) and/or linear methods for pattern assessment (e.g., FFT and EOF). While these metrics and techniques offer fundamental measures of the SST fields, they may not fully capture the complexity inherent in the most dynamic regions of the ocean. Motivated by advances in the analysis of natural images in computer vision, we employ a probabilistic autoencoder (PAE) which utilizes a Convolutional Neural Network (CNN) to learn the diversity of SST patterns. By design, the CNN learns the features most salient to the dataset, with built-in methodology to examine the image on a wide range of scales. Further, its non-linearity and invariance to translation offer additional advantages over Empirical Orthogonal Function (EOF) and like applications.
The ulmo algorithm is a PAE, a deep learning tool designed for density estimation. By combining an autoencoder with a normalizing flow, the PAE is able to approximate the likelihood function for arbitrary data while also avoiding a common downfall of flow models: their sensitivity to noisy or otherwise uninformative background features in the input Nalisnick et al. (2018). By first reducing our raw data (an SST field) to a compact set of the most pertinent learned features via the non-linear compression of an autoencoder, the PAE then provides an estimate of its probability by transforming the latent vector into a sample from an equal-dimension isotropic Gaussian distribution where computing the probability is trivial. We can then select the lowest probability fields as outliers or anomalous.
Our secondary goal of this manuscript, is to pioneer the process for like studies on other large earth science datasets in general and oceanographic datasets in particular including those associated with the output of numerical models. A similar analysis of SST fields output by ocean circulation models is of particular interest as an adjunct to the work presented herein. As will become clear, we understand some of the segmentation suggested by ulmo by not all of it. The method has also identified some anomalous events for which the basics physics is not clear. Assuming that the analysis of model-derived SST fields yields similar results, the additional output available from the model, the vector velocity field and salinity, as well as a time series of fields, will allow for a dynamic investigation of the processes involved.
This manuscript is organized as follows: Section 2 describes the data analyzed here, Section 3 details the methodology, Section 4 presents the primary results, and Section 5 provides a brief set of conclusions. All of the software and final data products generated by this study are made available on-line https://github.com/AI-for-Ocean-Science/ulmo.

DATA
With a primary goal to identify regions of the ocean exhibiting rare yet physical phenomena, we chose to focus on the L2 SST Aqua MODIS dataset (https://oceancolor.gsfc.nasa.gov/data/aqua/). The associated five minute segments, each covering ≈ 2000 × 1350 km of the Earth's surface and referred to as granules, have ≈ 1 km spatial resolution and span the entire ocean, clouds permitting, twice daily. For this study, we examined all nighttime granules from 2003-2019. The SST fields, the primary element of these granules, were processed by the Ocean Biology Processing Group (OBPG) at NASA's Goddard Space Flight Center, Ocean Ecology Laboratory from the MODIS radiometric data using the R2019 retrieval algorithm Minnett et al. (1-4 June 2020) and were uploaded from the OBPG's public server (https://oceancolor.gsfc.nasa.gov/cgi/browse.pl?sen=amod) to the University of Rhode Island (URI).
The method developed here requires a set of same-sized images. When exploring complex physical phenomena in the ocean, one is often interested in one of two spatial scales determined by the relative importance of rotation to inertia in the associated processes. The separation between these scales is generally taken to be the Rossby Radius of deformation, R o , which, at mid-latitude is ∼ O(30) km. Processes with scales larger than R o are referred to as mesoscale processes for which the importance of rotation dominates. At smaller scales the processes are referred to as sub-mesoscale. For this study, we chose to focus on the former and extracted 128 × 128 pixel images, which we refer to as cutouts, from the MODIS granules. Cutouts are approximately 128 km on a side. We are confident, supported by limited experimentation, that the techniques described here will apply to other scales as well.
The analysis was further restricted to data within 480 pixels of nadir. This constraint was added to reduce the influence of pixel size on the selection process for outliers; the along-scan size of pixels increases away from nadir as does the rate of this increase. To distances of ∼480 km the change in along-scan pixel size is less than a factor of two; at the edge of the swath the along-scan pixel size is approximately ∼ 5 times that at nadir.
The L2 MODIS product includes a quality flag -a measure of confidence of the retrieved SST -with values from 0 (best) to 4 (value not retrieved). The primary reason for assigning a poor quality to a pixel is due to cloud contamination although there are other issues that result in a poor quality rating Kilpatrick et al. (2019). A quality threshold of 2 was used for this study. Because the incidence, sizes, and shapes of clouds are highly variable (both temporally and spatially), an out-of-distribution (OOD) algorithm trained on images with some cloud contamination may become more sensitive to cloud patterns than unusual SST patterns. Indeed, our initial experiments were stymied by clouds with the majority of outlier cutouts showing unusual cloud patterns, suggesting an application of this approach to the study of clouds as well. To mitigate this effect, we further restricted the dataset to cutouts with very low cloud cover (CC), defined as the fraction of the cutout image masked for clouds or other image defects. After experimenting with model performance for various choices of CC, we settled on a conservative limit of CC ≤ 5% as a compromise between dataset size and our ability to further mitigate clouds (and other masked pixels) with an inpainting algorithm (see next section).
From each granule, we extracted a set of 128x128 cutouts satisfying CC ≤ 5% and distance to nadir of the central pixel ≤480 km. To well-sample the granule while limiting the number of highly overlapping cutouts, we drew at most one cutout from a pre-defined 32x32 pixel grid on the granule. This procedure yields ≈ 700, 000 cutouts per year and 12,358,049 cutouts for the full analysis.
Of course, by requiring regions largely free of clouds (CC < 5%), we are significantly restricting the dataset and undoubtedly biasing the regions of ocean analyzed both in time and space. Figure 1 shows the spatial distribution of the full dataset across the ocean. The coastal regions show the highest incidence of multiple observations, but nearly all of the ocean was covered by one or more cutouts. Given this spatial distribution, one might naively expect the results to be biased against coastal regions because these were sampled at higher frequency and comprises a greater fraction of the full distribution. This is mitigated, in part, by the fact that the non-coastal regions cover a much larger area of the ocean but, in practice, we find that a majority of the outlier cutouts are in fact located near land. Figure 1. Mollweide projection depicting the log10 of the spatial distribution of all cutouts analyzed in this manuscript. Note the higher incidence of data closer to land driven by the lower CC in those areas.

METHODOLOGY
In this section, we describe the preprocessing of the SST cutouts and the architecture of our ulmo algorithm designed to discover outliers within the dataset.

Preprocessing
While modern machine learning algorithms are designed with sufficient flexibility to learn underlying patterns, gradients, etc. of images (Szegedy et al. 2016), standard practice is to apply initial "preprocessing" to each image to boost the performance by accentuating features of interest, or suppressing uninteresting attributes. For this project, we adopted the following pre-processing steps prior to the training and evaluation of the cutouts.
First, we mitigated the presence of clouds. As described in § 2, this was done primarily by restricting the cutout dataset to regions with CC < 5%. We found, however, that even a few percent cloud contamination can significantly affect results of the OOD algorithm. Therefore, we considered several inpainting algorithms to replace the flagged pixels with estimated values from nearby, unmasked SST values. After experimentation, we selected the Navier-Stokes method Bertalmio et al. (2001) based on its superior performance at preserving gradients within the cutout. Figure 2 presents an example, which shows masking along a strong SST gradient (the white pixels between the red (∼ 22 • C) and yellow (∼ 19 • C) regions). We see that the adopted algorithm has replaced the masked data with values that preserve the sharp, underlying gradient without producing any obviously spurious patterns. Because inpainting directly modifies the data, however, there is risk that the process will generate cutouts that are preferentially OOD. However, we have examined the set of outlier cutouts to find that these do not have preferentially higher CC.
Second we applied a 3x1 pixel median filter in the along-track direction, which reduces the presence of striping that is manifest in the MODIS L2 data product. Third, we resized the cutout to 64x64 pixels using the local mean, in anticipation of a future study on ocean models, which have a spatial resolution of ≈ 2 km Qiu et al. (2019). Last, we subtracted the mean temperature from each cutout to focus the analysis on SST differences and avoid absolute temperature being a determining characteristic. We refer to the mean-subtracted SST values as SSTa. Figure 2. (left) Cutout which shows masking (white pixels) due primarily to sharp temperature gradients in this case, which tend to be flagged as low quality by the standard MODIS processing algorithm. (right) Same image but with masked pixels replaced by estimated values using the Navier-Stokes in-painting algorithm.

Architecture
ulmo is a probabilistic autoencoder (PAE), a likelihood-based generative model which combines an autoencoder with a normalizing flow. In our model, a deep convolutional autoencoder reduces an input cutout to a latent representation with N latent dimensions which is then transformed via the flow.
Flows Durkan et al. (2019) are invertible neural networks which map samples from a data distribution to samples from a simple base distribution, solving the density estimation problem by learning to represent complicated data as samples from a familiar distribution. The likelihood of the data can then be computed using the probability of its transformed representation under the base distribution and the determinant of the Jacobian of the transformation.
Though a flow could be applied directly to image cutouts in our use case, recent research Nalisnick et al. (2018) in the use of normalizing flows for OOD has revealed their sensitivity to uninformative background features which skew their estimation of the likelihood. To circumvent this issue, the PAE proposes to first reduce the input to a set of the most pertinent features via the non-linear compression of an autoencoder. The flow is then fit to the compressed representations of the image cutouts where its estimates of the likelihood are robust to the noisy or otherwise uninformative background features of the input image.
An alternative approach is the variational autoencoder (VAE) Kingma & Welling (2013) which provides a lower bound on the likelihood, though empirically we find PAEs boast faster and more stable training, and are less sensitive to the user's choice of hyperparameters.
Therefore, to summarize the advantages of our approach: (1) explicit parameterization of the likelihood function; (2) robustness of likelihood estimates to noisy and/or uninformative pixels in the input; and (3) speed and stability in training for a broad array of hyperparameter choices.
The key hyperparameters for the results that follow are presented in Table 1. Regarding N latent , we were guided by a Principal Components Analysis (PCA) decomposition of the imaging dataset which showed that 512 components captured > 95% of the variance. The full model with 4096 input values per cutout, is comprised of ∼ 6, 000, 000 parameters for the auto-encoder and ∼ 22, 000, 000 parameters for the normalizing flow. It was built with PyTorch and the source code is available on GitHub -https://github.com/AI-for-Ocean-Science/ulmo.

Training
Training of the complete model consists of two, independent phases: one to develop an autoencoder that maps input cutouts into 512-dimensional latent vectors, and the other to transform the latent vectors into samples from a 512-dimensional Gaussian probability distribution function (PDF) to estimate their probability. For the autoencoder, the loss function is the standard mean squared error reconstruction loss between all pixels in the input and output cutouts. In practice, the model converged to a small loss in ∼ 10 epochs of training.
The flow is trained by directly maximizing the likelihood of the autoencoder latent vectors. This equates to minimizing the Kullback-Leibler divergence between the data distribution and flow's approximate distribution. Minimizing this divergence encourages the flow to fit the data distribution and thereby produce meaningful estimates of probability.
Throughout training, we used a random subset of ≈ 20% of the data from 2010 (135,680 cutouts). These cutouts were only used for training and are not evaluated in any of the following results. Figure 3 shows an example of a preprocessed input SSTa cutout and the resultant reconstruction cutout from the autoencoder. As designed, the output is a good reconstruction albeit at a lower resolution that does not capture all of the finer features due to the information bottleneck in the autoencoder's latent space but it does capture the mesoscale structure of the field. For the normalizing flow, we used a cutout batch size of 64 and a learning rate of 0.00025. Similarly, we found ≈ 10 epochs were sufficient to achieve convergence. We performed training on the Nautilus distributed computing system with a single GPU. In this training setup, a single epoch for the auto-encoder requires 100 s while a single epoch for the flow requires ≈ 900 s.

RESULTS AND DISCUSSION
In this section, we report on the main results of our analysis with primary emphasis on outlier detection. We also begin an exploration of the ulmo model to better understand the implications of deep learning for analyzing remote-sensing imaging; these will be expanded upon in future works.  Figure 4 shows the LL distribution for all extracted cutouts modulo the set of training cutouts from 2010. The distribution peaks at LL ≈ 240 with a tail to very low values. The latter is presented in the inset which shows the lowest 0.1% of the distribution; these define the outlier cutouts of the full sample (or outliers for short).

The outlier cutouts sample
The striping apparent in the inset of Figure 4 indicates a non-uniform, temporal dependence in the outlier cutouts. Figure 5 examines this further, plotting the occurrence of outliers as a function of year and month. The only significant trend apparent is seasonal, i.e., a higher incidence of outliers during the boreal winter. We speculate this is due to Figure 5. Incidence (counts) of outlier cutouts broken down by month and year. The primary feature is seasonal, i.e., a higher number of outlier cutouts during boreal winter months than summer. There is also a weak, but possible increase in the incidence of outliers over the past ∼ 10 years.
the predominance of northern hemisphere cutouts/outliers -approximately 60%/64% of the total -and the reduced thermal contrast of northern hemisphere surface waters in the boreal summer. As will be shown, the range of SSTa in a cutout is correlated with the probability of the cutout being identified as an outlier; the larger the range the more likely the cutout will be so flagged. This is especially true in the vicinity of strong currents such as western boundary currents, which separate relatively warm, poleward moving equatorial and subtropical waters from cooler water poleward of the currents. In summer months the cooler water warms substantially faster than the surface water of the current dramatically reducing the contrast between the two water bodies, often masking the dynamical nature of the field in these regions rendering them less atypical. We also see variations during the ∼ 20 years of the full dataset, including a possible increase over the past ∼ 10 years. These modest trends aside, ulmo identifies outliers in all months and years of the dataset.
A question that naturally arises is whether there is any structure to the geographic distribution of outliers. Figure 6 shows the count distribution of the outliers across the entire ocean. Remarkably, the ulmo algorithm has rediscovered that the rarest phenomena occur primarily in western boundary currents -following the continental boundary and/or shortly after separation. These regions of the ocean have been studied extensively because of their highly dynamical nature. In short, the ulmo algorithm identified (or even rediscovered!) without any predisposition a consistent set of dynamically important oceanographic regions.
To a lesser extent, one also finds outliers in the vicinity of the connection between large gulfs or seas and the open ocean -the Gulf of California, the Red Sea and the Mediterranean. Also of interest are the outliers in the Gulf of Tehuantepec. These result from very strong winds blowing from the Gulf of Mexico to the Pacific Ocean through the Chivela Pass, resulting in significant mixing of the near-shore waters.
There are two ways to view the results in Figure 6: (1) as the contrarian, i.e., the ulmo algorithm has simply reproduced decades-old, basic knowledge in physical oceanography on where the most dynamical regions of the ocean lie; or (2) as the optimist, i.e., the ulmo algorithm -without any direction from its developers -has rederived one of the most fundamental aspects of physical oceanography. It has learned central features of the ocean from the patterns of SSTa alone. In this regard, ulmo may hold greater potential/promise to learn and derive other, not-yet-identified behaviors in the ocean. Figure 6. Depiction of the spatial distribution of the outliers discovered by ulmo. These are primarily in the well-known western boundary currents off Japan, North and South America, and South Africa. Note that the scaling is logarithmic. Figure 7 shows a gallery of 9 outliers selected to uniformly span time and location in the ocean. These exhibit extreme SSTa variations and/or complexity and (presumably) mark significant mesoscale activity. A common characteristic of these cutouts is the presence of a strong and sharp gradient in SSTa which separates two regions exhibiting a large temperature difference. Typically, such gradients are associated with strong ocean currents, often at mid-latitudes on the western edge of ocean basins. We define a simple statistic of the temperature distribution ∆T ≡ T 90 − T 10 where T X is the temperature at the Xth percentile of a given cutout. All of the outliers in Figure 7 exhibit ∆T > 7K, a point we return to in the following sub-section.

Scrutinizing examples of the outliers
As an example of the anomalous behaviour associated with outlier cutouts, we examine the evolution of the SST field in the vicinity of the 19 June 2010 cutout (Fig. 7) located in the Gulf Stream region; Fig. 8a shows the cutout and (b) its location in the 5-minute granule. We selected this cutout because it is in a region with which we have significant experience. Fig. 9 shows an expanded version of the SST field in the vicinity of the cutout. The main feature in Fig. 9 is the Gulf Stream, the bright red, fading to orange, band meandering from the bottom left hand corner of the image to the middle of the right hand side. A portion of the Gulf Stream loops through the lower half of the cutout and a streamer extends to the north (Fig. 8) from the northernmost excursion of the stream. To aid in the interpretation of this cutout, we make use of the mesoscale eddy dataset produced by Chelton et al. (2011). It shows an eddy, most probably a Warm Core Ring (WCR), moving to the west at approximately 5 km/day to the north of the stream from 17 May (very light gray circle) to 14 June (red circle) when it began to interact with the Gulf Stream, drawing warm Gulf Stream Water on its western side to the north and cold Slope Water on its eastern side to the south. The eddy disappears from the altimeter record two weeks later and is replaced by a very large anticyclone (the dotted black circle) to the west southwest of the eddy's last position. This is likely a detaching meander of the Gulf Stream resulting from the absorption of the eddy into an already chaotic configuration.
Of particular interest is that the Gulf Stream appears to have lost its coherence between approximately 63 • and 59 • W. Specifically, note the very thin band of cooler water (∼ 21 • C) in the middle of the warm band (∼ 24 • C) of, presumably, Gulf Stream Water between 63.5 • and 62 • W and a second similar band (but moving in the opposite direction) between 61.5 • and 60.5 • . The western cool band appears to separate one branch of Gulf Stream Water that has been advected from the southwestern edge of the large meander centered at 65.5 • W, 38.5 • N, and a second branch advected from its southeastern edge. These two branches may result from a general instability of the Gulf Stream associated with formation, or in this case the likely aborted formation, of a WCR. In the normal formation process, the initial state is a large meander of the Gulf Stream and the final state is a relatively straight Gulf Stream with a WCR to the north. In this case the process appears to have begun but inspection of subsequent images suggests that a ring was not formed; the meander reformed after initially beginning the detachment process. However, this is all quite speculative, the important point is that the stream appears to have lost its coherence immediately upstream of the cutout, which we believe to be a very unusual process. Admittedly the cutout only 'sees' a very small portion of this but we have found the suggestion of convoluted dynamics in the immediate vicinity of a large fraction of other outliers as well. Bottom line: Cornillon, who has been looking at SST fields derived from satellite-borne sensors for over 40 years, found that more than one-in-ten of the anomalous fields discovered by ULMO suggested intriguing dynamics that he has not previously encountered; recall that this is one-in-ten of one-in-a-thousand (the definition of an outlier) or approximately one field in ten thousand.

Digging Deeper
It is evident from the preceding sub-sections (e.g., Figure 7) that ulmo has discovered a set of highly unusual and dynamic regions of the ocean. Scientifically, this is extremely useful -irrespective of the underlying processes -as it can launch future, deeper inquiry into the physical processes generating such patterns. On the other hand, as scientists we are inherently driven to understand -as best as possible -what/how/why ulmo triggered upon. We begin that process here and defer further exploration to future work.
In Section 4.1, we emphasized that the entire gallery of outliers ( Figure 7) exhibits a large temperature variation ∆T > 7K. Exploring this further, Figure 10 plots LL vs. ∆T for the full set of cutouts analyzed. Indeed, the two are  anti-correlated with the lowest LL values corresponding to the largest ∆T . This suggests that a simple rules-based algorithm of selecting all cutouts with ∆T > 10 K would select the most extreme outliers discovered by ulmo. One may question, therefore, whether a complex and hard-to-penetrate AI model was even necessary to reproduce our results.
Further analysis suggests that there may be more to the distribution of LLs. Specifically, note that there is substantial scatter about the mean relation between LL and ∆T ; for example, at LL = −1000 one finds ∆T values ranging from 1 − 10 K. Similarly, any cutout with ∆T < 8 K includes a non-negligible set of images with LL 0. Figure 10 indicates that the patterns that ulmo flags as outliers are not solely determined by ∆T . Figure 10. Distribution of LL values as a function of ∆T . While there is a strong anti-correlation apparent, the relationship exhibits substantial scatter such that ∆T is not a precise predictor of LL nor the underlying SSTa patterns characterized by ulmo. The horizontal line at LL ≈ −1123 corresponds to the 0.1% threshold; cutouts with log-likelihood values beneath this line are considered to be outliers.
This becomes especially clear in the following exercise. Consider the full set of cutouts within the small range ∆T = [2 − 2.1] K. From Figure 10, we see these exhibit LL ≈ [−2600, 590] and find that the LL distribution is well described by a Gaussian (not shown) with < LL >≈ 10 and σ(LL) ≈ 150. Now consider the cutouts with the lowest/highest 10/90% of the distribution, i.e., the 'outlier'/'inlier' sub-samples within this small range of ∆T . We refer to these as LL 10 and LL 90 cutouts, respectively. Figure 11 shows the spatial distribution of these cutouts. Remarkably, there are multiple areas dominated by only one of the sub-samples (e.g., LL 90 cutouts along the Pacific equator). It is evident that ulmo finds large spatial structures in the log-likelihood distribution of cutouts that are independent of ∆T .
Furthermore, there are several locations in the ocean where LL 10 and LL 90 cutouts are adjacent to one another but still separate. One clear example is within the Brazil-Malvinas Confluence, off the coast of Argentina. Figure 12a shows a zoom-in of that region with the colors corresponding to the LL values (not strictly the LL 10 or LL 90 distributions shown in Fig. 11). Figure 12a highlights the clear and striking separation of the LL values in this region as do the histograms (Figure 12b) for the LL values of cutouts in the two rectangles shown in panel a. The dynamics of the ocean in this region is well-studied Piola et al. (2018). Higher LL regions tend to be found on the Patagonian Shelf where the dynamics are dominated by tides, buoyancy and wind-forcing the circulation at the local level-and off-shore currentsforcing the circulation remotely. In contrast, the lower LL regions track more dynamic, current-driven motions of the main Brazil-Malvinas Confluence. Of particular interest is the rather abrupt switch at ∼ 40 • S from higher LL values to the south to lower values to the north. This is consistent with the observation of Combes & Matano (2018), based on numerical simulations, that "[t]here is an abrupt change of the dynamical characteristics of the shelf circulation at 40 • S". They attribute this change in dynamics to this region being a sink for Patagonian Shelf waters, which are being advected offshore by the confluence of the Brazil and Malvinas Currents. Again, ulmo has captured striking detail in regional dynamics with no directed input. Further analysis of the region (not shown) suggests that ulmo has also captured seasonal differences in the dynamics, with a region of lower LL cutouts in waters approximately 100 m deep between 38 • and 45 • S in austral winter but not austral summer. Figure 11. Spatial distribution of the LL10/LL90 cutouts (red/blue) defined as the upper/lower tenth percentile of the LL distribution for the set of cutouts with ∆T = [2, 2.1] K. It is evident that these cutouts occupy distinct regions of the ocean, i.e., the ulmo algorithm has identified patterns with significant spatial coherence. More remarkable, note the several areas (e.g. in the Brazil-Malvinas current) where one identifies adjacent but separate patches of LL90 and LL10 cutouts.
Intrigued by ulmo's ability to spatially separate these regions based on SSTa patterns alone, we inspected a set of 25 randomly selected samples from R1, the eastern rectangle in Figure 12a, and 25 randomly selected samples from R2 to further explore its inner-workings (see lower panels of Figure 12). The comparison is striking and we easily identify qualitative differences in the observed patterns despite their nearly identical ∆T values. The higher LL cutouts show large-scale gradients and features with significant coherence whereas the lower LL cutouts exhibit gradients and features with a broader range of scales and a suggested richer distribution of relative vorticity.
Another area, which stands out in Figure 11, is that in the Northwest Atlantic where a region of LL 90 cutouts (red) are surrounded by LL 10 cutouts (blue). The structure (not shown) of the LL 90 cutouts in this region, which are on the Grand Banks of Newfoundland, resemble the structure of the LL 90 cutouts shown in Figure 12 and the structure of the LL 10 cutouts in this region is much closer to that of the LL 10 cutouts shown in Figure 12 than to the LL 90 cutouts in either region. In fact, a gallery of randomly selected LL 90 cutouts from the world ocean are similar to those off of Argentina and Newfoundland and a gallery of randomly selected LL 10 cutouts from the world ocean are more similar to the LL 10 cutouts off of Argentina and Newfoundland than to the LL 90 cutouts. Simply put, the structure of the SST cutouts shown in blue in Figure 11 tend to be similar to one another and quite different from those shown in red although the cutouts in both cases have virtually the same dynamic ranges in SST. This observation raises intriguing questions about the similarities and the differences in upper ocean processes in these regions -questions to be addressed in further analyses of the fields. (a) Distribution of the LL values for cutouts near the Brazil-Malvinas Confluence, restricted to those with temperature difference ∆T ≈ 2 K. One identifies a clear separation, where the lower LL values lie within the current and the higher values lie close to the Argentinian coast. Marked are two rectangles (R1, R2), one in each region, referred to in the other panels. Also marked is the 200m bathymetry described in the text. (b) Histograms of the LL values for the cutouts from two regions (R1/R2) chosen to show lower/higher LL values near the confluence. (lower panels) Representative cutouts from each subset -the left set are drawn from the R2 rectangle and therefore exhibit higher LL values. The right set are from R1. These galleries reveal qualitative differences in the SST patterns, i.e., unique ocean dynamics. The palette for the lower panels ranges linearly from −1 to 1 K.
We can capture some of the differences between the higher/lower LL sub-samples of Figure 12 with another simple statistic -the RMS in SSTa, σ T . On average, the lower LL cutouts exhibit ≈ 9% higher σ T than those with a higher LL. Furthermore, we find LL correlates with σ T in a fashion similar to ∆T . On the other hand, it is evident from Figure 12 that there is significant structure apparent in the cutouts that is not described solely by σ T . The correlations of LL with ∆T and σ T manifest from the underpinnings of ulmo: the distribution of ocean SSTa patterns reflect the distribution of simple statistics like ∆T or σ T , which exhibit large and non-uniform variations across the ocean. The complexity of these patterns, however, belie the information provided by simple statistics alone.

CONCLUSIONS AND FUTURE WORK
With the design and application of a machine learning algorithm ulmo, we set out to identify the rarest sea surface temperature patterns in the ocean through an out-of-distribution analysis yielding a unique log-likelihood (LL) value for every cutout. On this goal, we believe we were successful (e.g., Figures 6,7). In examining the nature of the outliers we found that these exhibited extrema of two simple metrics: the temperature difference ∆T and standard deviation σ T . With the full privilege of hindsight, we expect that any metric introduced to describe the cutouts which exhibits a broad and non-uniform distribution would correlate with LL. However, no single metric can capture the inherent pattern complexity and therefore none correlates tightly with LL ( Figure 10).
Looking to the future, the greatest potential of algorithms like ulmo may be that the patterns it learns are more fundamental than measures traditionally implemented in the scientific community (e.g., Fast Fourier Transform (FFT), Empirical Orthogonal Function EOF). We hypothesize that the mathematical nature of convolutional neural networks CNN -convolutional features and max-pooling, which synthesizes data across the scene while remaining invariant to translation -captures aspects of the data that EOF analysis could not (nor any other simple linear approach). Indeed, referring back to Figure 12, while as humans we trivially distinguish between the two sets of cutouts marking the ocean dynamics in the Brazil-Malvinas Confluence and can identify metrics on which they differ, these metrics offer incomplete descriptions. Going forward, we will determine the extent (e.g., via analysis of ocean model outputs) to which the patterns mark fundamental, dynamical processes within the ocean. Potentially, the patterns learned by ulmo (or its successors) hold the optimal description of any such phenomena.
As emphasized at the onset, this manuscript offers only a first glimpse at the potential for applying advanced artificial intelligence techniques to the tremendous ocean datasets obtained from satellite-borne sensors. The techniques introduced here will translate seamlessly to sea surface height or ocean color imaging to identify extrema/complexity in geostrophic currents and biogeochemical processes. These too will be the focus of future works. Warm Core Ring CONTRIBUTIONS Prochaska led the writing of the manuscript, including figure generation. He also ran the majority of models presented. Cornillon proposed the original idea of searching for extremes in the MODIS L2 SST dataset, undertook a significant fraction of the analysis of the resulting LL fields, guided the oceanographic interpretation of the results and contributed to the writing of the manuscript. Reiman prototyped and developed ulmo's deep learning components.