Multi-scale systems, such as the climate system, the atmosphere, and the ocean, are hard to understand and predict due to their intrinsic nonlinearities and chaotic behavior. Here, we apply a physics-consistent machine learning method, the multi-resolution dynamic mode decomposition (mrDMD), to oceanographic data. mrDMD allows a systematic decomposition of high-dimensional data sets into time-scale dependent modes of variability. We find that mrDMD is able to systematically decompose sea surface temperature and sea surface height fields into dynamically meaningful patterns on different time scales. In particular, we find that mrDMD is able to identify varying annual cycle modes and is able to extract El Nino–Southern Oscillation events as transient phenomena. mrDMD is also able to extract propagating meanders related to the intensity and position of the Gulf Stream and Kuroshio currents. While mrDMD systematically identifies mean state changes similarly well compared to other methods, such as empirical orthogonal function decomposition, it also provides information about the dynamically propagating eddy component of the flow. Furthermore, these dynamical modes can also become progressively less important as time progresses in a specific time period, making them also state dependent.

The climate system exhibits variability on a multitude of temporal and spatial scales. Due to the nonlinearity of the equation of motions, all these scales interact with each other, thereby hampering the understanding and predictability of the climate system. Here, we use multi-resolution dynamic mode decomposition (mrDMD), a physics-consistent machine learning approach, to systematically examine ocean variability of different time and spatial scales. We show that this method is able to systematically extract dynamically meaningful patterns of ocean variability.

## I. INTRODUCTION

The ocean covers about 72% of Earth’s surface and is an integral part of the global climate system. It provides essential environmental services, such as food and transportation, and affects atmospheric predictability and extreme events. Therefore, the ocean has a strong impact on vital aspects of our society, and it is of paramount importance to study and understand its underlying physical phenomena, which can vary on a multitude of time and space scales. One of the most important modes of ocean variability is the El Nino–Southern Oscillation (ENSO) (e.g., Timmermann *et al.*, 2018). ENSO describes variations in winds and sea surface temperature over the tropical eastern Pacific ocean, which have widespread effects on surface weather and climate conditions due to teleconnection patterns (Feldstein and Franzke, 2017). ENSO appears irregularly with enhanced frequency power in the range of 3–7 years.

Two other important modes of ocean variability are the western boundary currents: the Gulf Stream and the Kuroshio (Kang and Curchitser, 2013). The Gulf Stream varies on monthly through decadal time scales (Seidov *et al.*, 2019) and can shift between a northerly and southerly location (e.g., Pérez-Hernández and Joyce, 2014). The Kuroshio current tends to vary between a stable and an unstable state. For the former, the current is relatively strong and zonal, while in the latter case, the Kuroshio Extension tends to meander on large scales with higher eddy kinetic energy levels (e.g. Qiu and Chen, 2005 and Oka *et al.*, 2015).

A fascinating aspect of climate and ocean variability is that it occurs on all time scales and that due to the underlying nonlinear equations of motion, the different time scales interact with each other (Franzke *et al.*, 2020). However, this property also makes it difficult to understand climate and ocean variability because it is not straightforward to disentangle variability on different time scales. Moreover, another important aspect of understanding ocean variability is the identification of coherent structures with dynamical relevance, such as ENSO. A widely used method for the identification of modes of variability is empirical orthogonal functions (EOFs) (von Storch and Zwiers, 2003 and Hannachi, 2021), also known as a principal component analysis or a proper orthogonal decomposition. Global modes of sea surface temperature have been computed by Messié and Chavez, (2011). Global EOFs identify the well-known modes of ocean variability, such as the ENSO, the Pacific decadal oscillation (PDO), and the Atlantic multidecadal oscillation (AMO). However, the dynamical relevance of some of these modes has been questioned (Clement *et al.*, 2015 and Mann *et al.*, 2020). While EOFs are a powerful tool for multivariate data analysis, they have the drawback that the EOF patterns are mutually orthogonal. Thus, the EOF patterns lose physical interpretability since the ocean modes, or any physical modes, need not to be mutually orthogonal (North, 1984).

Hence, there is a need for better methods, which are able to systematically identify dynamically relevant patterns. A promising method is dynamic mode decomposition (DMD) (Tu *et al.*, 2014; Kutz *et al.*, 2016a; Rowley *et al.*, 2009; Mezić, 2005; 2013; and Brunton *et al.*, 2016), a machine learning method. DMD decomposes high-dimensional fields into complex patterns whose eigenvalues describe the growth rates and oscillation frequencies of the modes. DMD is widely used in many areas (Kutz *et al.*, 2016a; Tu *et al.*, 2014; Rowley *et al.*, 2009; Brunton *et al.*, 2016; and Mezić, 2005; 2013) and recently also in geophysical and climate research (Kutz *et al.*, 2016b; Gottwald and Gugole, 2019; and Gugole and Franzke, 2019). DMD is a dimension reduction method, which has a strong theoretical and dynamical underpinning. For a given high-dimensional time series, DMD computes a set of complex modes; each of these modes represents an oscillation with a fixed frequency and a growth rate. For linear systems, these modes are analogous to normal modes. Furthermore, DMD is closely connected to principal oscillation patterns and linear inverse models (Hasselmann, 1988; Penland and Magorian, 1993; and Tu *et al.*, 2014). However, DMD is more general. DMD approximates the modes and eigenvalues of the Koopman operator and, thus, can represent nonlinear dynamics (Tu *et al.*, 2014 and Kutz *et al.*, 2016a). DMD is different from other popular dimension reduction methods, such as EOFs. EOFs are not directly associated with a temporal behavior, while DMDs are. However, in contrast to EOFs, DMD modes are not orthogonal. Hence, DMDs might provide a less parsimonious description of the full data set than EOFs, but on the other hand, DMDs are dynamically more meaningful, which we will show below.

The multi-scale space-time structure of ocean variability calls also for multi-scale methods. The multi-resolution DMD is an attractive option for this problem since it provides a systematic multi-scale decomposition into dynamical modes (Kutz *et al.*, 2016a; 2016b). Here, we will demonstrate that multi-resolution DMD is able to identify dynamically meaningful patterns of ocean variability, which are of practical concern. While our study does not advance theory, our aim is to demonstrate the ability of DMD to systematically extract multi-scale dynamics from a complex real-world system, the ocean. Furthermore, our study shows how DMD can be used to deepen our understanding of ocean dynamics, and we also show how DMD can systematically extract multi-scale dynamics of a component of the climate system, which not many methods can do. We also demonstrate by extracting the changing annual cycle of SST how DMD can potentially lead to better predictions of the climate system. In Sec. II, we describe the ocean data sets we are using and the multi-resolution DMD method. In Sec. III, we present our results for global SST and sea surface height dynamics in the Kuroshio and Gulf stream. In Sec. IV, we summarize our study results.

## II. DATA AND METHODS

### A. Data

To demonstrate the abilities of DMD and to examine ocean variability, we use two datasets. The first one is the extended reconstructed sea surface temperature (ERSST) version 5 data set. This is a global monthly mean sea surface temperature (SST) data set on a $2\xb0\xd72\xb0$ regular horizontal grid (https://doi.org/10.7289/V5T72FNM) (Huang *et al.*, 2017). The data set covers the period January 1854–December 2020. This data set allows us to examine variability on medium to long time scales.

The second data set is the Aviso satellite altimetry sea surface height (SSH). This is a global daily sea surface height data set on a $0.25\xb0\xd70.25\xb0$ regular horizontal grid (Saraceno *et al.*, 2008) (https://resources.marine.copernicus.eu/?option=com_cswtask=results?option=com_cswview=detailsproduct_id=SEALEVEL_GLO_PHY_L4_REP_OBSERVATIONS_008_047 and https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/mss.html). The data set covers the period 1993–2018. This data set allows us to examine ocean eddies and larger-scale flow structures on short to medium time scales. Because of the higher temporal and spatial resolution for this data set, it is computationally challenging to consider the whole globe as a domain. Furthermore, scientifically, a better understanding of the dynamics of smaller scale features, such as eddies and meanders, are also needed. Hence, we focus on two important ocean currents, the Gulf Stream and the Kuroshio, and apply mrDMD to these two areas.

We define the Gulf Stream region as the area covering $280\xb0$E–$340\xb0$E, $30\xb0$N–$60\xb0$N and the Kuroshio region as the area covering $120\xb0$E–$170\xb0$E, $25\xb0$N–$50\xb0$N. For spatial pattern correlations, we use bandpass filtering of the SSH data using a Fourier transformation based approach where the cut-off frequencies correspond to the respective mrDMD frequency bands. The bandpass filtering is necessary for the SSH data because the eddy scale changes considerably for different time scales. Thus, pattern correlations between fast time-scale mrDMD patterns and the full flow fields would lead to small correlation values. On the other hand, no filtering is necessary for SST data because on monthly time scales, the anomalies are still relatively large scale. For the SST data, we also tested the impact of detrending the data. Detrending the data leads to qualitatively similar results to using non-detrended data. Hence, our results are robust.

### B. Dynamic mode decomposition

We consider the following dynamical system (Kutz *et al.*, 2016a):

where $x$ denotes the state vector, $t$ time, $\mu $ the parameters of the system, and $f$ is the possible nonlinear function representing the dynamics. Equation (1) can also induce a discrete-time representation for time step $\Delta t$,

In general, it is impossible to derive a solution to the nonlinear system equation (1). DMD takes an equation-free, machine learning view where we have no knowledge of the dynamics of the system. DMD only uses observed data from the system to approximate and forecast the system. Hence, DMD computes an approximate, locally linear, dynamical system,

with discrete-time representation,

with

where the subscript $k$ denotes discrete time. The solution of this system can be represented in terms of its eigenvalues $\lambda j$ and corresponding eigenvectors $\varphi j$ of the discrete-time matrix $A$,

where $b$ is the matrix of the initial conditions $bj$ and $j$ is an index. $\Phi $ is the matrix consisting of the eigenvectors $\varphi j$. DMD now derives a low-rank eigen-decomposition of $A$ that optimally captures the trajectory of the system in a least-squares sense so that the object

is minimized across all grid points, and this is achieved by an eigen-decomposition of $A$.

The DMD algorithm is as follows: The data can be described by two parameters.

- $\u2219$
n: number of spatial grid points per time step and

- $\u2219$
m: number of time steps.

We now have the following two sets of data:

so that $x\u2032k=F(xk)$ for time step $\Delta t$. The DMD modes now correspond to the eigen-decomposition of $A$. $A$ relates the data $X\u2032\u2248AX$, and thus, $A=X\u2032X\u2020$, where $\u2020$ is the Moore–Penrose pseudo-inverse (Kutz *et al.*, 2016a).

The DMD method has a strong theoretical underpinning, as it is connected to the Koopman operator (Tu *et al.*, 2014 and Kutz *et al.*, 2016a). DMD is a finite dimensional approximation of the modes of the Koopman operator. The Koopman operator is an infinite dimensional linear operator describing the dynamics of nonlinear systems.

The Koopman operator is defined as follows (Kutz *et al.*, 2016a):

Consider a continuous–time dynamical system,

where $x\u2208M$ is a state on a smooth n-dimensional manifold $M$. The Koopman operator $K$ is an infinite-dimensional linear operator that acts on all observable functions $g:M\u2192C$ so that

The Koopman operator propagates states along with the flow $F$.

### C. Multiresolution dynamic mode decomposition

Multiresolution dynamic mode decomposition (mrDMD) is an advanced DMD method for analyzing multi-scale systems, such as the ocean and the atmosphere (Kutz *et al.*, 2016a; 2016b). Basically, it performs DMD on different time scales, similar to a wavelet analysis (Lau and Weng, 1995 and Kutz *et al.*, 2016a). Figure 1 shows a schematic of the mrDMD approach. mrDMD starts with analyzing the full time series and by identifying the slowest modes of variability; then, this window is divided into two equally long windows as displayed in Fig. 1, and the DMD analysis is repeated. This is recursively repeated until the fastest dynamics are reached in the data set. The frequency of the modes is given by the eigenvalues of the mrDMD modes. As a cut-off frequency, we choose that slow modes can only perform a maximum of two oscillations in a window in order to eliminate faster modes from this level. See Kutz *et al.* (2016a; 2016b) for more details. The eigenvalues of the mrDMD modes are related to frequencies as follows:

where $s=1\rho 8\pi \Delta t$ with $\rho =2/T$, where $T$ is the window length. Only mrDMD modes with frequencies $\omega $ smaller than $\rho $ are considered at a given level. The power of the mrDMDs is computed as in Jovanović *et al.*, (2014) and Kutz *et al.*, (2016a) and is based on the dynamics. The mrDMD power is computed by separating the DMD amplitude into a product of the normalized DMD modes, a diagonal matrix of mode amplitudes, and the Vandermonde eigenvalue matrix (Jovanović *et al.*, 2014 and Kutz *et al.*, 2016a). The Vandermonde matrix captures the exponentiation of the DMD eigenvalues. The exponentiated eigenvalues determine the power of the DMD modes. Since DMDs are not normalized, the amplitude of the power spectra does not correspond to a physical unit.

We have the following notation: mrDMD(i,j,k) denotes the DMD from the ith level and the jth segment, while k denotes the number of the corresponding DMD mode. For instance, mrDMD starts at the first level, i.e., the full time series. The second level denotes the two halves of the full time series. The first half is the first segment, while the second half is the second segment. The first DMD mode of the third level and the second segment is denoted: mrDMD(3,2,1). As for EOFs, the sign of the DMDs is arbitrary.

## III. RESULTS

### A. Global sea surface temperature

We start with the monthly global SST data set. The mrDMD power spectrum (Fig. 2) reveals that the maximum power is contained in the first level, while the second largest power is in the range containing the annual cycle. The first mrDMD mode of the first level, the mode whose real component of the eigenvalue is closest to one and, thus, corresponds to an almost neutral mode, has a geographical structure, which is very similar to the climatological mean state [Figs. 3(a) and 3(b)]. This level of the mrDMD decomposition has two more mrDMD modes, which have purely real eigenvalues [Figs. 3(d) and 3(f)]. These are almost neutral modes, though their eigenvalues are somewhat smaller than 1 with values of 0.9712 and 0.8911 and, thus, are damped modes. mrDMD(1,1,2) is likely representing low-frequency behavior of the sea ice edge since in comparison with mrDMD(1,1,1), it represents an equatorward extension in the polar regions of negative anomalies and to a reduction of the meridional temperature gradient.

Our mrDMD modes of the first level are also different from the global SST EOF modes of the study by Messié and Chavez, (2011) in that they do not correspond to the well-known modes of ocean variability, such as the ENSO, PDO, or AMO. The comparison with EOFs is not straightforward since with mrDMD, we focus on certain time scales at each level, while EOFs do not systematically distinguish between different time scales, though they tend to be ordered by an integrated auto-correlation time scale (Franzke *et al.*, 2005 and Franzke and Majda, 2006).

We also computed pattern correlations by projecting the instantaneous SST fields onto the mrDMD modes [Figs. 3(f)–3(h)]. The pattern correlations also confirm the almost neutral behavior of these three mrDMD modes, though mrDMD(1,1,2) and mrDMD(1,1,3) show a slight reduction of correlation strength, which could be an imprint of their damped nature.

A powerful feature of mrDMD is that it can identify the modes of the annual cycle, which are encoded at level 8. Figure 4 shows the periods of the mrDMD modes associated with the annual cycle, which are all around 12 months, but also vary, indicating that mrDMD has the ability to capture year to year variations in the annual cycle. Note that our segments do not correspond to the calendar annual cycle; the annual cycle dynamics are determined by the dynamics of the climate system. As an example of an annual cycle mrDMD mode, we choose the tenth segment of the eighth mrDMD level, i.e., mrDMD(8,10,2) (Fig. 5). The other annual cycle related mrDMDs of these other segments look similar, suggesting that this is a robust feature (not shown). The corresponding first mrDMD mode corresponds to the mean state over this segment and is an almost neutral mode. The second mrDMD mode, consisting of a real and imaginary component, corresponds to the annual cycle. Both components together make up a propagating mode with temperature anomalies of opposite sign in both hemispheres where the imaginary part corresponds to the transition seasons, while the real part corresponds to the peak seasons [Fig. 5(e)]. Our results are consistent with the analysis by Pezzulli *et al.*, (2005), which found substantial interannual variability in the seasonal cycle of the univariate NINO3.4 index.

The length of the annual cycle seems to be determined by other large-scale processes. In Figs. 5(f) and 5(g), we display composites of SST over long [Fig. 5(f)] and short [Fig. 5(g)] annual cycle events, respectively. The composites are averaged over those segments when the mrDMD annual cycle period is either one standard deviation above or below its long term mean period of about 12 months, respectively. The composites indicate a hemispheric seesaw behavior of SST between both hemispheres.

The mrDMD method is also able to identify ENSO events as is demonstrated in Fig. 6. For instance, mrDMD(7,52,3) corresponds to the El Nino of 1987. This mrDMD mode has a period of about 14 months, which is in the range of the typical El Nino duration of between 7 and 24 months. Both parts of mrDMD(7,52,3) show the typical El Nino anomaly in the tropical Pacific. This shows that mrDMD is able to extract physically meaningful patterns from ocean data sets. This is consistent with the results of Kutz *et al.*, (2016a). Since ENSO is a transient phenomenon, mrDMD represents this as real and imaginary components of a DMD mode. This is in contrast to EOFs, in which ENSO would be represented by just one EOF pattern (Messié and Chavez, 2011). This illustrates that mrDMD provides patterns, which are dynamically directly interpretable since DMD analysis provides eigenvalues determining the pattern’s oscillation frequency and growth rate.

Two widely recognized modes of SST variability are the Atlantic multidecadal oscillation (Knight *et al.*, 2006 and Ting *et al.*, 2011) and the Pacific decadal oscillation (Mantua and Hare, 2002). The mrDMD power spectrum (Fig. 2) does not show enhanced power at decadal time scales; the power at those scales is actually rather low. This is consistent with recent studies, which questioned the physical relevance of these modes, which are identified by a global EOF analysis (Messié and Chavez, 2011). Mann *et al.*, (2020) provide evidence that both modes are not distinguishable from the noise background. Also, Clement *et al.*, (2015) provide evidence that the AMO is not a dynamical oceanographic mode of variability since their model experiments do not contain ocean dynamics but still show AMO type variability. These studies are consistent with our mrDMD results that those modes are potentially not dynamically meaningful.

To further demonstrate the ability of mrDMD to identify physically meaningful patterns, we now examine variability at the fourth and fifth levels in more detail. As Fig. 2 shows, the fourth level has one big amplitude event and the fifth level has five large amplitude events. First, we focus on the large amplitude event of the fourth level (Fig. 7). As can be seen from Fig. 2, this event occurred between November 1874 and October 1895. In Fig. 7(a), we display the average over this period. This period is characterized by warm SST anomalies in the North Pacific along 40$\xb0$N, in the Labrador sea and the Fram strait, and in the Southern Ocean between South America and Antarctica. Most of the remaining ocean is anomalously cold, especially the Arctic Ocean. In contrast, the anomalies averaged over all other times are much weaker [Fig. 7(b)]; whether this is just due to averaging over a longer period or whether this suggests that DMD picks systematically dynamically relevant and active states needs further research, ideally with very long climate model data. mrDMD(4,2,2) has a period of about 20 years. Its real component is similar to the Pacific decadal oscillation (Mantua and Hare, 2002), but this mrDMD describes more complex dynamics than just a standing pattern.

We now turn to the fifth level where we have five large amplitude events (Fig. 2). We now average over the periods of these five events [Fig. 8(a)] and average over all other times [Fig. 8(b)]. The high amplitude composite shows increased SST over most of the ocean areas with cold anomalies only in the northern North Pacific and the south of Greenland. The composite of all other times displays mainly cold anomalies. The modes mrDMD(5,3,2) and mrDMD(5,6,2) again demonstrate that those states are the result of dynamic processes. Both mrDMD modes also have similarities to the PDO. This suggests that the PDO is an important mode of ocean variability on decadal time scales. Moreover, the fact that we identify multiple DMD modes resembling the PDO is consistent with the finding that the PDO is not a single physical mode of variability, but rather is an aggregation of multiple processes, such as ENSO teleconnections, reemergence of SST, and stochastic atmospheric forcing (Newman *et al.*, 2003; 2016; Qiu *et al.*, 2007; Schneider and Cornuelle, 2005; and Vimont, 2005).

### B. Kuroshio sea surface height

In the following, we examine the Kuroshio current using mrDMD. For this purpose, we use daily Aviso SSH data. The mrDMD power frequency–time plot (Fig. 9) shows that again, the maximum power is contained in the first level. However, the third, fifth, and sixth levels contain a sizable amount of power. The power of the first level is associated with the mean state (Fig. 10) as mrDMD(1,1,1) corresponds to the mean state [compare Figs. 10(a) and 10(c)]. The mrDMD(1,1,2) mode projects onto the linear trend [Fig. 10(b)] for most of the area with the exception of the southern area of our chosen box, though with the opposite sign. The pattern correlation is positive and has a trend toward zero [Fig. 10(f)]. This suggests that mrDMD(1,1,2) is a damped mode, which is also indicated by its positive real eigenvalue with modulus smaller than 1. The pattern itself represents weakening of the central SSH gradient of the Kuroshio and has a positive anomaly at the location of the Kuroshio large meander, southeast of Japan.

To focus on specific high amplitude events in the frequency plot of Fig. 9, we display the relevant mrDMDs of the third level in Fig. 11. The time scales of this level correspond to about 7.5 years. The by far largest and dominating segment at this level is the second segment (04.07.1999–03.01.2006), while all other segments are rather inconspicuous when it comes to the power spectrum. mrDMD(3,2,1) corresponds to the local mean state for that time segment. While mrDMD(3,2,1) of this time period is very similar to the overall mrDMD(1,1,1) in Fig. 10, one noteworthy aspect is the variations in the pattern correlations [Fig. 11(b)]. For most of the time, the correlations fluctuate between 0.955 and 0.96. There are two excursions to values around 0.94, corresponding to the years 1999 and 2001. These years were characterized by an exceptionally meandering current with a pair of large persistent eddies off the coast of Japan (see Fig. 2 in Qiu and Chen, 2005), i.e., a northern warm core eddy and a southern cold core eddy. If we now take higher modes into consideration, these exceptional dynamics are confirmed [Figs. 11(c)–11(h)]. The modes mrDMD(3,2,2) and mrDMD(3,2,4) have complex eigenvalues and are, thus, propagating patterns with periods of about 986 and 874 days. Correspondingly, mrDMD(3,2,3) and mrDMD(3,2,5) are the complex conjugates of mrDMD(3,2,2) and mrDMD(3,2,4), respectively. All these modes highlight a propagating large meander around 30–35$\xb0$N, which is characteristic for this time period. Furthermore, modes 2 and 3 propagate the signal of the eddy pair between 32–37$\xb0$N and 140–145$\xb0$E, which is such a dominating feature for both 1999 and 2001 (see Fig. 2 of Qiu and Chen, 2005). A part of this signal is also visible in modes 4 and 5, although the dominant role of these modes seems to lie in the general meandering of the current starting from this dipole eddy perturbation. To support this interpretation, the corresponding pattern correlations have been computed by projecting the mrDMD patterns onto bandpass filtered SSH fields where the bandpass filter frequencies correspond to the mrDMD frequencies associated with the respective third level modes, which correspond here to periods between 874 and 986 days. One important observation here is that the real and complex modes remain at a 90$\xb0$ angle with respect to the correlations for the first 1.5 to two cycles between positive and negative correlations, which corresponds to the period up until the end of 2001. After that, they become increasingly mixed and also damped with respect to correlation amplitude. This suggests that these specific dynamical modes only play an important propagating role during the first years of this time period, coinciding with the years of exceptional large-scale and persistent eddy and meandering activity of the Kuroshio.

To conclude the mrDMD analysis of the Kuroshio SSH, we would like to point out that this diagnostic does not necessarily distinguish between the stable and unstable years of the Kuroshio (as discussed by Qiu and Chen, 2005). Instead, our results suggest that the mrDMD analysis discriminates between years, which are dominated by propagating large-scale anomalies and those years where the conditions are either more persistent (very long time scales) or shaped by (short lived) chaotic behavior on temporal and spatial scales that do not correspond to the respective mrDMD level. A potential drawback of mrDMD here is the strict, non-overlapping decomposition of the total time period, which can lead to specific events being split between different segments. The consequence may be that the mrDMD does not pick up those events at certain levels but may hint at them at higher levels when the data are further subdivided. One solution for this potential drawback would be to use an overlapping windowing approach. However, this would be computationally much more expensive.

### C. Gulf Stream sea surface height

Next, we examine the Gulf Stream SSH. The multiresolution DMD power spectrum is displayed in Fig. 12. Both the frequency-time and time averaged frequency plots show that the lowest frequencies dominate the spectrum. This is mainly due to the mean state, which is captured by mrDMD(1,1,1) [compare Figs. 13(a) and 13(b)]. The climatological mean state and mrDMD(1,1,1) are very similar so that most of the low-frequency information is associated with the mean. As the eigenvalue of mrDMD(1,1,1) is 1.0, it represents a temporally neutral mode.

In order to have a closer look at the structure of the associated mrDMD patterns, we also compute pattern correlations between the respective mrDMD pattern and the SSH fields. Figure 13 shows these time lag correlations for mrDMD(1,1,1) of level 1. The pattern shows for the whole time pattern correlation values between 0.94 and 0.98. For mrDMD(1,1,2), the pattern correlation shows, in absolute terms, a decreasing trend for the first 20 years before it is increasing to large absolute values again. mrDMD(1,1,2) may show an imprint of long time-scale changes of the Gulf Stream due to climate change signals and low-frequency changes in the Gulf Stream intensity and stability. Even though it has a robust correlation of around 0.3, the features in the pattern are rather small scale with some slight imprint of an emphasized north–south gradient along the Gulf Stream path. It is, therefore, a mixture of eddy driven and large-scale changes in the structure, intensity, and position of the Gulf Stream.

Similarly to the mrDMD of the Kuroshio region in Fig. 9, the Gulf Stream SSH also exhibits isolated large amplitude events in frequency space for low-frequencies (0.0008; about 3.5 years), although these occur at higher frequencies than the level 3 mrDMDs described in Sec. III B for the Kuroshio, which occurred at periods of about 7 years. These events are associated with a strong meandering of the Gulf Stream (not shown). The next interesting mrDMD occurs at the sixth level during the period 20.04.2000–08.02.2001 (Fig. 14). mrDMD(6,9,1) again corresponds to the mean state for this time window. mrDMD(6,9,2) is again a standing pattern since it has a purely real eigenvalue. The pattern correlation of this mode reveals that this mode shifts the Gulf Stream from south to north and back to south over the period 20.04.2000–08.02.2001 with an emphasis between around 48$\xb0$W and 72$\xb0$W. The pattern has relatively close resemblance to the first EOF mode discussed by Pérez-Hernández and Joyce, (2014) (their Fig. 3). The strengthening and weakening of the correlations with mrDMD(6,9,2) coincide with a noticeable northward shift of the Gulf Stream in October 2000, also discussed by Pérez-Hernández and Joyce, (2014) (their Fig. 4). mrDMD(6,9,3) and mrDMD(6,9,5) correspond to eddy modes with periods of about 200 and 162 days, respectively. The pattern correlations confirm that these mrDMDs are propagating eddy fields as the correlations have a periodic structure and are shifted by about 90$\xb0$ between real and imaginary components (Fig. 14). The mrDMD of the sixth level, therefore, highlights this dynamic shifting event by a peak in the DMD power (Fig. 12) for the ninth segment and then splits the event into a mean component, a Gulf Stream shift, and further eddy dynamics with a propagating signature.

According to Pérez-Hernández and Joyce, (2014), other extreme northward shifts of the Gulf Stream occurred in July 1995 and April 2012. These events also show high DMD power. As already noted above, the ability of DMD to capture such events depends on the length of the segments (i.e., the frequency), the dynamical nature of these shifts, which may be picked up by low DMD modes and whether these events are fully captured within one segment or whether they may be split between two consecutive segments. One, therefore, needs to look at individual segments to investigate why the power of that segment is comparatively large or low. Also, pattern correlations is a useful diagnostic for the interpretation of the DMDs.

## IV. SUMMARY

In this study, we have demonstrated that the physics-consistent machine learning method multi-resolution dynamic mode decomposition is able to extract dynamically relevant patterns of ocean variability. We applied mrDMD to sea surface temperatures and sea surface height fields. We find that mrDMD is able to systematically decompose SST and SSH fields into meaningful patterns on different time scales. This allows for a systematic analysis of multi-scale systems and the climate system, in particular. We show that mrDMD is able to identify annual cycle modes, which can vary from year to year, without supervision. This is an important aspect in the analysis of climate dynamics since a time-varying annual cycle can provide an alternative basic state for the study of climate anomalies (Wu *et al.*, 2008 and Pezzulli *et al.*, 2005). When using a time fixed annual cycle, all changes, e.g., due to global warming, will be part of the anomalies. However, changes in the annual cycle can have pronounced impacts on the climate system and, therefore, our understanding of it.

We also show that mrDMD is able to extract actual ENSO events from the SST data set without supervision. ENSO is one of the most important modes of climate variability and occurs on a broad range of time scales (Timmermann *et al.*, 2018). This makes mrDMD a very attractive method for the analysis of multi-scale systems since no *a priori* filtering is necessary. mrDMD also seamlessly provides a decomposition into a local basic state and eddy fields allowing state-dependent eddy-mean flow interaction studies. Here, we used non-overlapping time windows for our mrDMD analysis, putting some constraints on identifying specific events (such as ENSO events) if they are spread across two windows. However, this can be relaxed to an overlapping windows analysis. mrDMD is also relatively computationally inexpensive, making it an attractive analysis method and potentially a prediction method. For instance, Gottwald and Gugole, (2019) use DMD to identify regime transitions in the North Atlantic region and the Southern Hemisphere. DMD also has the potential for subgrid-scale modeling as shown by Gugole and Franzke, (2019).

Considering sea surface height fields for the Kuroshio and Gulf Stream, mrDMD is capable of identifying dynamically interesting and complex events related to changes in the position and intensity of the currents. While it can highlight these mean state changes similarly well compared to other methods, such as EOF decomposition, it also provides information about the dynamically propagating eddy component of the flow. Such dynamically evolving components consist of a real and an imaginary DMD mode, whose correlations with the original SSH data are cyclic and are shifted by 90$\xb0$ with respect to each other, signifying a propagating signal. As highlighted by the weakening of the correlations, these dynamical DMD modes can also become less important as time progresses throughout a specific time window. The detailed mrDMD decomposition of the flow allows us to investigate isolated events and associate relevant drivers to the respective modes.

## ACKNOWLEDGMENTS

C.L.E.F. was supported by the Institute for Basic Science (IBS), Republic of Korea (No. IBS-R028-D1) and the Pusan National University research grant 2021. S.J. was supported by subprojects M3 and L4 of the Collaborative Research Center TRR181 Energy Transfer in Atmosphere and Ocean funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) (Project No. 274762653).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Christian L. E. Franzke:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Federica Gugole:** Conceptualization (equal); Writing – review & editing (equal). **Stephan Juricke:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## REFERENCES

*Nonlinear and Stochastic Climate Dynamics*, edited by C. L. E. Franzke and T. O’Kane (Cambridge University Press, 2017), pp. 54–104.

*et al.*, “