Unsupervised Learning Reveals Geography of Global Ocean Dynamical Regions

Abstract

Dynamically similar regions of the global ocean are identified using a barotropic vorticity (BV) framework from a 20‐year mean of the Estimating the Circulation and Climate of the Ocean state estimate at 1° resolution. An unsupervised machine learning algorithm, K‐means, objectively clusters the standardized BV equation, identifying five unambiguous regimes. Cluster 1 covers 43 ± 3.3% of the ocean area. Surface and bottom stress torque are balanced by the bottom pressure torque and the nonlinear torque. Cluster 2 covers 24.8 ± 1.2%, where the beta effect balances the bottom pressure torque. Cluster 3 covers 14.6 ± 1.0%, characterized by a “Quasi‐Sverdrupian” regime where the beta effect is balanced by the wind and bottom stress term. The small region of Cluster 4 has baroclinic dynamics covering 6.9 ± 2.9% of the ocean. Cluster 5 occurs primarily in the Southern Ocean. Residual “dominantly nonlinear” regions highlight where the BV approach is inadequate, found in areas of rough topography in the Southern Ocean and along western boundaries.

Motivation

Before the advent of modern observational and modeling techniques, understanding of the physical/dynamical state of the ocean focused on large‐scale quasi‐laminar descriptions such as Sverdrup balance, abyssal recipes, or Stommel‐Arons flows (Munk, ; Munk & Palmén, ; Stommel, ). Recent advances in instrumentation and modeling capability have revealed that ocean physics is characterized by complex spatial and temporal variability. It is possible that every location in the ocean has a unique physical state depending upon many factors, including local topography, meteorology, proximity to eastern and western boundaries, or latitude, rendering any global scale interpretation lacking in general applicability.

The ocean is spatially and temporally diverse, but delineating spatial and temporal commonalities and continuities are central to understanding emergent patterns that lead to a global geography of dynamical regimes. Dominant physics underlying the emergent patterns become evident when common features are identified. This note's purpose is to explore unsupervised machine learning as a method for depicting and understanding the gross features of global oceanic physics. The study is restricted to the barotropic vorticity (BV) balance of a time mean global circulation as calculated from a noneddy resolving state estimate (Fukumori et al., ). Our approach appears to be both interesting and useful and is readily generalized to far more complex oceanic states. In unsupervised machine learning approaches, data are not prelabeled or pregrouped (Bishop, ). The noneddy resolving case is explored in this initial work, analogous to Coupled Model Intercomparison Project Phases 5 and 6 (CMIP5 and CMIP6) efforts (Church & Miles, 2013; Eyring et al., ; Stouffer et al., ). The presented analysis is intended to provide a description relevant to the CMIP climate models.

The presence of differing dynamical regimes is already suggested by the known structures of the wind‐stress forcing and the geometry of the ocean basins, including the underlying topography. Classifying and identifying regions in the world ocean is done here using the variations in the dominant terms of the BV budget, following the demonstration that the global budget can be closed. Schoonover et al. () and Yeager () have assessed the dynamics of the BV budget focusing on the North Atlantic Ocean variability and sensitivity to resolution, respectively. Here the procedures are global.

Distinct geographical physics were demonstrated by, for example, Xu and Fu (), who used altimeter data to show differing spatial regimes of geostrophic turbulence. Their global patterns are presumably connected to circulation structures, planetary waves, and topography. Similarly, Hughes and Williams () and Sonnewald et al. () used linear statistical models to infer patterns suggestive of global regimes. The present work extends the pattern determination methodology, as working only with data from the surface limits the application to a comprehensive assessment of global dynamical regimes. Use of the Estimating the Circulation and Climate of the Ocean (ECCO) state estimate extends the previous surface data analyses to the full three dimensions of the ocean.

Objective classification via K‐means clustering allows unbiased identification of patterns in data. This form of unsupervised machine learning is common in many fields, ranging from pharmaceuticals to engineering (Breuhl et al., ; Hauser & Rybakowski, ; Kulis & Jordan, ). Ardyna et al. () applied a similar method to identify regions with distinct biological activity, and K‐means have been used to identify key regions for data collection to build maps of nitrate in the Southern Ocean (Y.‐C. Liang, personal communication, February 15, 2018). Applications in the earth sciences have been explored both in the prognostic and diagnostic sense (Krasnopolsky et al., ; Schneider et al., ).

Methods

The ECCO Version 4 State Estimate

The BV equation is applied to version 4, release 2, of the ECCO (ECCOv4) global state estimate, described by Forget et al. () and Wunsch and Heimbach (), see also ECCO Consortium (, ). The state estimate has a nominally 1° resolution, available at ecco.jpl.nasa.gov. A least squares with Lagrange multipliers (4DVAR) approach is used to obtain the state estimate. The result is a free‐running version of the MIT General Circulation Model (MITgcm; Adcroft et al., ), with adjusted initial and boundary conditions and internal model parameters. The ECCO state satisfies basic conservation laws for enthalpy, salt, volume, and momentum while remaining largely within error estimates of a diverse set of global data (Stammer et al., ; Wunsch & Heimbach, , ). Regions without data are brought into consistency in a dynamically consistent way using the dynamics, still relying on parameterizations but avoiding the use of untested statistical hypotheses or infilling e.g., Reynolds et al. ().

Barotropic Vorticity

The momentum and continuity equations of an ocean in a thin shell on a rotating sphere are $${\partial }_t\mathbf{u}+f\mathbf{k}\times \mathbf{u}=-\frac{1}{{\rho }_0}{\nabla }_{h}p+\frac{1}{{\rho }_0}{\partial }_z\tau +\mathbf{a}+\mathbf{b},{\partial }_{z}p=-g\rho ,$$ $${\nabla }_h·\mathbf{u}+{\partial }_{z}w=0$$

Pressure, gravity, density, and vertical shear stress are denoted p, g, ρ, and τ, respectively, with ρ₀ the reference density; the three‐dimensional velocity field v = (u, v, w ) = (u, w); the gradient ∇ = (∇_h, ∂_z); the unit vector is denoted k; planetary vorticity is a function of latitude ϕ in $f\mathbf{k}=(0,0,2\mathrm{\Omega }\sin \varphi )$; the viscous forcing by vertical shear is denoted ∂_zτ; the nonlinear torque is a, and the horizontal viscous forcing b includes subgrid‐scale parameterizations. Assuming a steady state, the vertical integral from the surface z = η(x,y,t) to the water depth below the surface z = H(x,y) is $$\beta V=\frac{1}{{\rho }_0}\nabla p_b\times \nabla H+\frac{1}{{\rho }_0}\nabla \times \mathrm{𝛕}+\nabla \times \mathbf{A}+\nabla \times \mathbf{B},$$ where ∇·U = 0,U·∇f = βV, the bottom pressure is denoted p_b, $\mathbf{A}={\int }_H^{\eta }\mathbf{a}\mathrm{d}z$, and $\mathbf{B}={\int }_H^{\eta }\mathbf{b}\mathrm{d}z$. The curl operator ∇× yields a scalar, representing the vertical component of the operator. The left‐hand side of equation  is the planetary vorticity advection term, while the right‐hand side of equation  is the bottom pressure torque (BPT), the wind and bottom stress curl, the nonlinear torque, and the viscous torque, respectively. The subgrid‐scale parameterization introduces a torque, which is included in the viscous torque term. Nonlinear torque is composed of three terms: $$\nabla \times \mathbf{A}=\nabla \times \left[{\int }_{-H}^{\eta }\nabla ·(\mathbf{u}\mathbf{u})\mathrm{d}z\right]+{\left[w\zeta \right]}_{z=H}^{z=\eta }+{\left[\nabla w\times \mathbf{u}\right]}_{z=H}^{z=\eta },$$ where uu is a second‐order tensor. The right‐hand side of equation  represents the curl of the vertically integrated momentum flux divergence, the nonlinear contribution to vortex tube stretching, and the conversion of vertical shear to barotropic vorticity. Horizontal viscous forcing includes that induced by subgrid‐scale parameterizations. Twenty‐year averaged fields of the BV equation are used after a Laplacian smoother is applied, with an effective averaging range of three grid cells.

Unsupervised Learning: K‐Means Clustering

Our goal is to determine the spatial patterns that correspond to various balances between the dominant terms of the BV equation. Various combinations of terms dominate the BV equation in different regions, and spatial patterns emerge, which we seek to determine. Clustering identifies groups in data based on how the data are distributed in parameter space, the dimensions of which are defined by a set of chosen variables (or “features” as they are also called in the machine learning literature such as ; Kubat, ). In this application, the terms of the BV equation are used to define the dimensions of the parameter space in which the clusters are identified. If groups of dominant terms are present that differ significantly, a robust separation into distinct “clusters” is feasible. For each cluster, the area‐weighted histogram is calculated of the values of each term in the BV equation. Next, the clusters are sorted by the amount of geographical area that they cover using a sorting algorithm, associating the name of the cluster as listed in Table . This focuses on those clusters that featured clearly different balances between terms in the BV equation.

Percentage of Area Covered by the Area‐Specific Balance of the BV Equation  and the Corresponding Map Figure

Cluster	Area	Leading terms
1	43 ± 3.3%, Depth coherent (Figure a)	∇ × τsb + ∇ × A≈∇pb × ∇H (Figure b)
2	24.8 ± 1.2%, Interior flow (Figure c)	∇ × τsb≈∇pb × ∇H + ∇·(fU) (Figure d)
3	14.6 ± 1%, Quasi‐Sverdrupian (Figure e)	∇ × τsb≈∇·(fU) (Figure f)
4	6.9 ± 2.9%, Interior flow, vertical (Figure a)	∇ × τsb≈∇·(fU) + ∇pb × ∇H (Figure b)
5	1.9 ± 1%, Interior flow, Southern Ocean (Figure c)	∇ × τsb≈∇·(fU) + ∇pb × ∇H (Figure d)
6–50	8.9 ± 0.3%, Dominantly nonlinear (Figure e)	∇·(fU)≈∇ × A + ∇ × τsb (Figure f)

Note. Leading order terms are sorted by magnitude, colors indicating if barotropic vorticity is added (red in font) or removed (blue in font) by the leading order term, the corresponding bar chart figure shows the full breakdown. The quoted percentage coverage and StD is the mean of 100 runs of the algorithm.

The terms in the BV equation were normalized so that each term individually has a zero mean and unit variance globally. This scaling of the variables ensures that the patterns in the variance of the variables form gridpoint to gridpoint are what is highlighted, rather than the relative magnitudes (Not applying normalization and feature scaling produced a poor result where clusters are harder to robustly identify.). Using the normalized and scaled fields, the clusters are more easily identified. The K‐means algorithm (MacQueen, 1967) involves an iterative minimization of the sum of squares of the Euclidean distance partitioning of the hyperspace given by the terms in the BV equation: $$J={\mathrm{\Sigma }}_{j=1}^K{\mathrm{\Sigma }}_{i=1}^n||{\mathbf{x}}_i^j-{\mathbf{c}}_j|{|}^2$$ the number of K clusters is a free parameter that is chosen a priori; the cluster centers have random initial values scattered throughout the parameter space. The parameter x is a vector field that is defined at every grid cell on a discretized sphere, with each element x_i representing a five‐dimensional vector on the model's horizontal grid, such that index i uniquely identifies a grid point on the sphere, with (lon,lat) = (ϕ_i,θ_i). The components (features) of each vector x_i correspond to the five terms in the BV budget. Each cluster j = 1,…, K is represented by the five‐dimensional characterizing vector c_j, and the K‐means classification attributes each vector x_i to a unique cluster c_j, thus ${\mathbf{x}}_i={\mathbf{x}}_i^j$. The distance between a data point ${\mathbf{x}}_i^j$ and c_j is given by $||{\mathbf{x}}_i^j-{\mathbf{c}}_j|{|}^2$. Each data point is associated with the closest K‐cluster, the position of c_j is recalculated, and the association reassessed until the solution converges. Assumptions regarding the covariance of the data are discussed in the appendix.

The solution is sensitive to the initialization and choice of K, and the algoritwhm partitions the parameter subspace using linear hyperplanes. This linearity constraint means that higher numbers of K can both assist in partitioning the subspace more appropriately, and isolate noise. The appendix demonstrates the small sensitivity of the result to the algorithm's initial random seed, and the impact of varying K. An optimal value of K is determined as K > 35 using the Akaike and Bayesian Information Criteria (AIC and BIC; Akaike, 1973). Information criteria provide a measure of the quality of a statistical model, rewarding increased likelihood across a data set and penalizing overfitting. AIC and BIC indicate robust regimes as they both asymptote in the bottom left panel of Figure , suggesting that no information is gained by further increasing K. K = 50 is used for the remaining analysis, where five clusters are individually analyzed as they are considered to represent somewhat classical dynamical regimes, and the remaining 45 clusters are taken together as a single “nonlinear regime.”

Results

The closure in ECCOv4 for the 20‐year average of the BV terms in equation  is illustrated in Figure . Individual terms are of order ±10⁻⁹ m s⁻¹ and the residual has magnitude of less than ±10⁻¹² m/s. For 36% of the ocean the residual is ≪1% (Figure ). This small residual permits going forward with confidence. Some numerical issues do exist on the continental shelf and in shallow water generally, but these regions only amount to 3% of the area of the global ocean and will be ignored.

The breakdown of the barotropic vorticity budget (m/s) over 1992–2013 in the ECCOv4 State Estimate. From equation (3), Figure 1a is the planetary vorticity advection term, Figure 1b is the bottom pressure torque (BPT), Figure 1c is the wind and bottom stress curl, Figure 1d is the nonlinear torque, and the viscous torque is Figure 1e.

Figure b illustrates where the beta term is important from equation . This term is balanced by the BPT term shown in Figure c and the wind and bottom stress BV terms shown in Figure d. The remainder is largely found in the nonlinear BV contributions seen in Figure e, with the lateral viscous dissipation largely being an order of magnitude smaller, apart from localized regions in the Southern Ocean. Wind and bottom stress BV terms in Figure d are largely zonally symmetric, with large patterns of negative BV in the Southern Ocean, and large gyre patterns visible in the Pacific and Atlantic basins. BPT in Figure c is associated with interactions with steep bathymetry. For example, in the Southern Ocean a large positive patch leads toward the Antarctic‐Pacific ridge, with a negative patch beyond. This structure is consistent with vortex stretching as circulation crosses the ridge. Along Western Boundaries, BPT is positive to the west and negative just adjacent to the east, consistent with studies such as Myers et al. (). The BV of the nonlinear torque shown in Figure e is concentrated along the western edge of basins where WBCs are found, but it is less spatially coherent than the BPT term. Large activity stands out in the Southern Ocean region, particularly in the Atlantic sector. Lateral viscous dissipation is small.

Picking out globally coherent dynamical regimes, the K‐means algorithm results are presented in Figure a, where the numbering on the color bar is arbitrary. The structure is mainly found in five named regimes, summarized in Table , which are briefly examined. Each is numbered, and a partially descriptive label is attached. Figure c shows the area and 2σ uncertainty covered by the named clusters.

Top figure (a) illustrates the area selected by the clusters. Colors represent clusters in arbitrary order. The depth coherent ocean region (43%, Figure a) in dark blue, Interior flow (24.8%, Figure c) in light brown, Quasi‐Sverdrupian (14.6%, Figure e) in light green, Interior flow, vertical, (6.9%, Figure a) in dark green, Interior flow, Southern Ocean, (1.9%, Figure c) in lighter blue, and the dominantly nonlinear torque cover remaining 8.9% (Figure e). Panel (b) illustrates that the Akaike Information Criteria (AIC) and Bayesian Information Criteria (BIC) asymptoting and a K of 50 is chosen for the analysis. Error bars of 2σ capturing the stochastic seed are shown. Panel (c) demonstrates the robustness of the algorithm with the ocean area, with 100 runs of the classification algorithm finding nearly identical areas (2σ error bars).

Figures and isolate the geographical area (left column) for each cluster determined as distinct by the K‐means algorithm. From this geographical area, the area‐weighted histogram is computed (right column) across all terms of the BV equation. These histograms relate the clusters to different dynamical regimes, illustrated in Figure . Area averaging was done for comparison.

Maps of the selected locations (a, c, e) and corresponding area averaged histogram (b, d, f) of the terms in the barotropic vorticity equation. The color bar is kept, but the color/ordering of the map are arbitrary. Colors in the barchart indicate if barotropic vorticity is added (red) or removed (blue).

Maps of the selected locations (a, c, e) and corresponding area averaged histogram (b, d, f) of the terms in the barotropic vorticity equation. The colorbar is kept, but the color/ordering of the map are arbitrary. Colors in the barchart indicate if barotropic vorticity is added (red) or removed (blue).

The largest cluster, accounting for 43% of the global ocean (Cluster 1 in Figure a), is determined by a balance between wind stress curl and bottom torque. This depth‐coherent “negative wind curl/bottom torque” region is found primarily in zonal streaks in the tropics, and in a thin ribbon in the Southern Ocean mainly in the Pacific sector. In the Northern Hemisphere, Cluster 1 areas surround the subtropical and subpolar gyres. Large areas of the Arctic Seas are also in this Cluster. Figure b demonstrates that the balance of terms is dominantly between the input of negative vorticity by the wind‐stress curl largely balanced by the positive input by the bottom interaction terms.

The next largest dynamical region covers 25% of the ocean area (Cluster 2, Figure c: Interior flow “positive wind curl/beta and bottom torque”), where the wind stress curl inputs positive vorticity, nearly balanced by the beta and bottom interaction terms. In the Northern Hemisphere, this cluster covers the southern extent of the subpolar gyres. A zonal streak crosses the equator in both the Atlantic and Pacific, but is absent in the Indian Ocean. The Southern Hemisphere has large Cluster 2 expanses in both the Pacific and Atlantic, but again not in the Indian Ocean.

The 15% of the ocean area selected by Cluster 3 are illustrated in Figure e (Quasi‐Sverdrupian “negative wind torque/beta effect”). In Sverdrup balance, the wind torque and the beta effect are the only important terms. This theoretical relation is seen in Cluster 3 in the subtropical gyres where they are expected, together with regions in the Southern Ocean that the classical theory did not consider. Dominant areas in the subtropical gyres in the Northern Hemisphere Atlantic and Pacific stand out, together with thin streaks on the equator. Isolated streaks are seen in the Southern Ocean, and in a large area of the Southern Hemisphere tropical Indian Ocean. This region might be considered also as corresponding to quasi‐Sverdrup balance.

Cluster 4 (Figure a: Interior flow, vertical, “positive wind torque, beta and bottom stress”), covers 7% of the ocean and is a complement to Cluster 1. In the Northern Hemisphere, the Cluster largely represents the northern edge of the subpolar gyre. In the Southern Hemisphere, it is found on the eastern edge of the Pacific and Atlantic basins, just to the south, and flaring out westward of the continental barrier. In the Indian Ocean, this barrier can be seen to be New Zealand. The area of this dynamical regime fills the subtropical Indian Ocean down to the border with the Southern Ocean, where this regime is absent. Figure b illustrates that it is an amplified version of the dominant terms seen in Figure d, being an order of magnitude larger, but still having the wind as the major source of barotropic vorticity, with sinks in the Coriolis term and BPT. A small source exists in the nonlinear torque.

The “Southern Ocean gyre” is Cluster 5, covering only 2% of the global world ocean seen in Figure c. This Cluster is mainly seen in a series of zonal streaks in the Southern Ocean with negative wind torque and complements Cluster 4. Again, nonlinear torque is a small sink.

A summary of the area of the remaining clusters that account for 9% of the world ocean is shown in Figure e: “Dominantly nonlinear.” Separate clusters have different colors. Areas of rough bathymetry stand out, such as the Pacific‐Antarctic Ridge and the Drake Passage area. Figure f is the overall average, illustrating that the nonlinear contribution to the barotropic vorticity dominates, together with the Coriolis term. The different constituents are quite varied, but strong contributions from the nonlinear torque are consistently present. Their detailed discussion is the subject of a subsequent study.

Discussion and Conclusions

In the ECCOv4 state estimate, the barotropic vorticity equation closes very accurately. The 20‐year time average ECCOv4 state estimate, was analyzed globally using K‐means clustering to find regions dominated by groups of balanced terms in the BV equation. Large regions of the global ocean exhibit structures of consistent term balances as shown in Figure . Those balances vary among the wind stress, Coriolis, and BPT terms. Areas where the nonlinear torque are small suggest that the linearized BV is a good approximation. Areas where the nonlinear torque are important are found in western boundary regions, as well as the Southern Ocean where the Antarctic Circumpolar Current interacts with bathymetric obstacles. The momentum‐dominated area (Cluster 1) implies a coherent vertical structure. Cluster 3 mainly in the subtropical gyre is unique in lacking significant contributions by BPT, implying it is shielded from topography. Transition zones have a stronger momentum‐driven portion of the BPT and topographic interactions to become important. Cluster 4 has a stronger baroclinic component to the BPT, feeling topography. The Southern Ocean Cluster is like Cluster 2, but with contributions of opposite sign. Important nonlinear contributions are present in the remaining ocean, and the linearized barotropic interpretation is not appropriate.

Schematic of identified regions with names and cluster numbers. The depth‐coherent area implies a coherent vertical structure in Cluster 1. The quasi‐Sverdrupian gyre in Cluster 3 is unique due to lack of bottom pressure torque (BPT). The interior flow, vertical, in Cluster 4 has a stronger momentum driven portion of the BPT, and topographic interactions begin to become important. The interior flow in Cluster 2 has a stronger baroclinic component to the BPT and feels topography. The interior flow, Southern Ocean, in Cluster 4 is like the interior flow in Cluster 2, but with contributions of opposite sign. The remainder is dominated by nonlinear contributions, and the barotropic interpretation is not appropriate.

In the North Atlantic Ocean, results are generally consistent with the inferences of the relative magnitude of the terms of the BV equation (Schoonover et al., ; Yeager, ). Cluster analysis reveals a shift from a barotropic flow in Clusters 1 and 3, to a strong interior flow (baroclinic meridional, North Atlantic Current, and North Atlantic Deep Water, and flow over the Mid Atlantic Ridge) in Clusters 2 and 4. Globally, the clustering illustrates that strong interior flow is present in vast expanses of the Southern Hemisphere, as well as in the North Pacific. Cluster 4 coincides with regions identified by Lumpkin and Speer (), Perez et al. (), and Speich et al. () as areas of water mass transformation and intermediate pathways in the overturning circulation between surface and deep water. The quasi‐Sverdrupian regime in Cluster 3 is not present in the South Atlantic and Pacific. In the Southern Ocean Clusters 3 and 5 dominate, with significant nonlinear contributions.

Five regions cover 91% of the world ocean. Residual areas collected here as “dominantly nonlinear” have a small spatial extent, but are dynamically important in the overall circulation. These regions include the Drake Passage region as well as the Antarctic‐Pacific Ridge regions where the circulation interacts with topography and cross frontal transport likely takes place. In the Northern Hemisphere, areas in the Labrador Sea and on the continental shelf stand out as nonlinear. These nonlinear regions will be the subject of a separate study at higher resolution.

The sign and spatial distribution of the wind stress term suggests the importance of Ekman pumping (negative) or suction (positive). Equatorial and Southern Ocean regions show Ekman pumping, whereas the subpolar gyre areas have Ekman suction where mode waters are created. The BPT term mirrors the wind stress term, suggesting it acts as either a source or a sink in opposite complement to the wind stress. A lack of symmetry in wind driven gyres in the Southern and Northern Hemisphere show that the gyres are not driven solely by the sign of the Ekman pumping. This complex relationship among the terms is the subject of future study.

The use of vertical integrals to describe the circulation is a simplification of what is a three‐dimensional problem, as is the use of a model in which the important eddy field is only included through parameterizations such as in CMIP5 and CMIP6 (Church & Miles, 2013; Eyring et al., ; Stouffer et al., ). Future work is intended to apply this and related machinery to fully eddy‐resolving ocean states.

A Appendix K‐Means and influence of Information Criteria

The K‐means algorithm is related to methods such as Principal Component Analysis (PCA), more traditionally applied to oceanography. Where PCA attempts to represent all data vectors using a low‐order combination of eigenvectors, minimizing the mean squared reconstruction error, the K‐means algorithm represents the data vectors via a small number of clusters. This is also done to minimize the mean squared reconstruction error. In this manner, the K‐means algorithm can be interpreted as a very sparse PCA.

Robustness of the regions in terms of the stochastic initialization is highlighted in Figure c, where the K‐means clustering was run 100 times and mean and 2σ used as metrics in Table . The regimes identified are robust, with the extent of the subpolar gyre being the main area where the algorithm shows appreciable variance.

The K‐means algorithm is initiated by scattering K first guesses of where the parameters/clusters could be. This initial guess introduces a stochastic element. The success of the algorithm is sensitive to K, as this determines how the hyperspace given by the dimensions is partitioned. As with regression analysis, adding parameters can increase the accuracy, but overfitting should be avoided. Determining the appropriate value of K, information criteria (AIC and BIC) are used to assess the quality of the statistical model. These measures weight the added accuracy with the cost of adding additional parameters, minimizing the expectation of the prediction error, are used: $$\text{AIC}=2K-2\ln (\mathcal{L}),$$ $$\text{BIC}=K\ln (n)-2\ln (\mathcal{L}),$$ where n is the number of data points and $\mathcal{L}$ is the likelihood: $$\mathcal{L}={\mathrm{\Pi }}_{i=1}^N\frac{1}{\sqrt{2\mathrm{\pi }{\sigma }^2}}\exp \left(-\frac{{({\zeta }_i-{\widehat{\zeta }}_i)}^2}{2{\sigma }^2}\right).$$

The parameter ζ_i is the observed, and ${\widehat{\zeta }}_i$ is the prediction, so ${({\zeta }_i-{\widehat{\zeta }}_i)}^2$ are the prediction residuals. In the estimate, the AIC value is minimized, which determines the smallest appropriate order to represent the time series. As discussed by Priestley () and Yang (), the AIC can overestimate the order. Figure b demonstrates that both the AIC and BIC stabilize at >35 K and the asymptotic nature of the regime.

The Euclidian distance is used, meaning the variance is assumed to be isotropic (meaning round). This leads to the standard practice of normalizing and standardizing data. To elucidate the impact of assumptions the algorithm makes for the classification, a more generalized form of clustering was also tested: Gaussian Mixture Models (GMM). GMM are used to assess the impact of assumptions relating to the covariance structure; spherical, diagonal, tied, or full covariance. Using the BIC, the results building on the BV data were not seen to be sensitive to this. However, this could be important at higher resolution as the K‐means clustering problem is NP‐hard and GMM could perform better.

Acknowledgments

This work was funded by the U.S. National Aeronautics and Space Administration Sea Level Change Team (contract NNX14AJ51G) and through the ECCO Consortium funding via the Jet Propulsion Laboratory. M. S. acknowledges Anne Reinarz, Roosa Tikkanen, and Katherine Rosenfeld as well as the python scikit‐learn toolbox. ECCOv4, release 2, model output is available in this website (ftp://mit.ecco-group.org/ecco_for_las/version_4/release2/). ECCOv4, release 2, model output is available in this website (ftp://mit.ecco-group.org/ecco_for_las/version_4/release2/). M. S. developed the concept, designed the method, and performed the numerical simulations. M. S. wrote the paper under the guidance of C. W. C. W. and P. H. contributed equally to the final version of the manuscript.