Multiscale fractal dimension analysis of a reduced order model of coupled ocean–atmosphere dynamics

Atmosphere and ocean dynamics display many complex features and are characterized by a wide variety of processes and couplings across different timescales. Here we demonstrate the application of multivariate empirical mode decomposition (MEMD) to investigate the multivariate and multiscale properties of a reduced order model of the ocean–atmosphere coupled dynamics. MEMD provides a decomposition of the original multivariate time series into a series of oscillating patterns with time-dependent amplitude and phase by exploiting the local features of the data and without any a priori assumptions on the decomposition basis. Moreover, each oscillating pattern, usually named multivariate intrinsic mode function (MIMF), represents a local source of information that can be used to explore the behavior of fractal features at different scales by defining a sort of multiscale and multivariate generalized fractal dimensions. With these two complementary approaches, we show that the ocean–atmosphere dynamics presents a rich variety of features, with different multifractal properties for the ocean and the atmosphere at different timescales. For weak ocean–atmosphere coupling, the resulting dimensions of the two model components are very different, while for strong coupling for which coupled modes develop, the scaling properties are more similar especially at longer timescales. The latter result reflects the presence of a coherent coupled dynamics. Finally, we also compare our model results with those obtained from reanalysis data demonstrating that the latter exhibit a similar qualitative behavior in terms of multiscale dimensions and the existence of a scale dependency of the statistics of the phase-space density of points for different regions, which is related to the different drivers and processes occurring at different timescales in the coupled atmosphere–ocean system. Our approach can therefore be used to diagnose the strength of coupling in real applications.

1 Introduction

The atmosphere and the ocean form a complex system whose dynamical variability extends over a wide range of spatial and temporal scales (Liu, 2012; Xue et al., 2020). As an example, the tropical regions are markedly characterized by inter-/multi-annual processes like the El Niño–Southern Oscillation (ENSO) (Neelin et al., 1994; Meehl et al., 2003), while the North Atlantic Oscillation (NAO) affects extratropical northern hemispheric regions at seasonal and decadal timescales (Ambaum et al., 2001). The sources of these processes have been widely investigated by means of multiple data analysis methods and various types of modeling (e.g., Philander, 1990; Czaja and Frankignoul, 2002; Van der Avoird et al., 2002; Mosedale et al., 2006; Kravtsov et al., 2007; Feliks et al., 2011; Liu, 2012; L'Hévéder et al., 2014; Farneti, 2017; Vannitsem and Ghil, 2017; Wang, 2019; Xue et al., 2020, and references therein), highlighting how the atmospheric low-frequency variability (LFV) is related to the ocean. The latter develops thanks to the interaction with the ocean mixed layer (OML) driven by a mixing process due to the development of an instability within the water column (Czaja and Frankignoul, 2002; D'Andrea et al., 2005; Wunsch and Ferrari, 2004; Gastineau et al., 2012) that also shows a strong seasonal variability. The relation between the OML and the LFV can be investigated from a dynamical system point of view by developing suitable reduced order ocean–atmosphere models dealing with the modeling of the coupling between the atmosphere and the underlying surface layer of the ocean. Recently, by means of a 36-variable model displaying marked LFV Vannitsem et al. (2015) demonstrated that the LFV in the atmosphere could be a natural outcome of the ocean–atmosphere coupling. Other sources could be invoked to explain and to contribute to the development of LFV in the atmosphere, such as the long-range system memory as a consequence of the heat storage mechanism of the land–ocean–atmosphere system (e.g., Lovejoy, 2021; Lovejoy et al., 2021), the internal dynamics of the atmosphere itself (e.g., Legras and Ghil, 1985), or even the interaction between the tropical and extratropical regions (e.g., Alexander et al., 2002; Vannitsem et al., 2021), just to quote a few.

The current work presents an investigation on how a recently introduced concept of multiscale generalized fractal dimensions can be used to analyze the statistics of attractors in coupled ocean–atmosphere systems (Alberti et al., 2020a). This demonstration is done by means of the reduced order model developed in Vannitsem et al. (2015). Indeed, the dynamical properties of physical systems can be related to their support fractal dimension as well as its singularities by means of different established concepts like the box-counting dimension (e.g., Steinhaus, 1954; Mandelbrot, 1967; Ott, 2002), generalized correlation integrals (Grassberger, 1983; Hentschel and Procaccia, 1983; Pawelzik and Schuster, 1987), the pointwise dimension method (Farmer et al., 1983; Donner et al., 2011), and related characteristics (Badii and Politi, 1984; Primavera and Florio, 2020). These methods are based on partitioning the phase space into hypercubes of size ℓ to define a suitable invariant measure through the filling probability of the ith hypercube by N_k points as $p_k=N_k/N$, with N being the total number of points. With M(ℓ) denoting the number of filled hypercubes, we can define some useful dynamical invariants such as the box-counting (or capacity or simply fractal) dimension

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{D}_{\mathrm{0}}\doteq -\underset{\mathrm{\ell }\to \mathrm{0}}{lim}\underset{N\to \mathrm{\infty }}{lim}{\displaystyle \frac{\log M\left(\mathrm{\ell }\right)}{\log \mathrm{\ell }}},\end{array}$$

the information dimension

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{D}_{\mathrm{1}}\doteq \underset{\mathrm{\ell }\to \mathrm{0}}{lim}\underset{N\to \mathrm{\infty }}{lim}{\displaystyle \frac{{\sum }_{k=\mathrm{1}}^{M\left(\mathrm{\ell }\right)}{p}_k\log p_k}{\log \mathrm{\ell }}},\end{array}$$

and the correlation dimension

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{D}_{\mathrm{2}}\doteq \underset{\mathrm{\ell }\to \mathrm{0}}{lim}\underset{N\to \mathrm{\infty }}{lim}{\displaystyle \frac{\frac{\mathrm{1}}{N^{\mathrm{2}}}{\sum }_{i\ne j}\mathrm{\Theta }\left(\mathrm{\ell }-|x_i-x_j|\right)}{\log \mathrm{\ell }}},\end{array}$$

with Θ(⋯) being the Heaviside function. More specifically, D₀ is a measure of the sparseness of the phase space by the studied system's dynamics, D₁ is a measure of the information gained on the phase space with a given accuracy, while D₂ is a measure of correlations, i.e., mutual dependence, between phase-space points. All these fractal dimension measures, as well as their higher order extensions D_q measuring qth order correlations between points in the phase space, have been used to characterize the statistics of the phase-space scaling of a given system (Hentschel and Procaccia, 1983). More details on the estimation of generalized fractal dimensions are provided in the Supplement. However, the above concepts only give us a global view on the phase-space system's properties, without exploring how these evolve at different scales in the real space (Alberti et al., 2020a). More recently, by means of a suitable combination between a state-of-the-art time series decomposition method (the empirical mode decomposition) and the concept of generalized fractal dimensions, Alberti et al. (2020a) introduced a multiscale approach to deal with the investigation of the evolution of the statistics of the phase-space scaling in dynamical systems.

Here, we extend for the first time the concept of multiscale generalized fractal dimensions in a multivariate framework by means of the multivariate empirical mode decomposition (MEMD), allowing us to investigate the multiscale and multivariate properties of a reduced order model of the ocean–atmosphere coupled dynamics. By using the oscillating patterns forming the decomposition basis of the MEMD algorithm, usually named multivariate intrinsic mode functions (MIMFs), we define the new concept of multiscale and multivariate generalized fractal dimensions. The MEMD results allow us to capture the essential dynamics of the phase-space trajectory that can be used for reconstructing the skeleton of the phase-space dynamics, while the evaluation of the fractal dimensions at different timescales provides a quantitative characterization of the intrinsic complexity of oscillating patterns that can be related to the attractor properties. Our results also allow for associating the statistics of the phase-space scaling to the dynamical regimes at different timescales of the coupled ocean–atmosphere system. Finally, our findings for the reduced order model well reconcile with corresponding results for reanalysis data, thus supporting and encouraging the use of reduced order models for investigating the essential aspects of the coupled ocean–atmosphere system in terms of the statistics of the phase-space scaling.

2 The reduced order ocean–atmosphere model

Reduced order coupled ocean–atmosphere models are key tools in the hierarchy of climate models, allowing for an extensive analysis of the features of the coupled dynamics that would otherwise be impossible to evaluate (Lorenz, 1984; Nese and Dutton, 1993; Roebber, 1995; Jin, 1996; Timmermann et al., 2003; Van Veen, 2003; De Cruz et al., 2016; Vannitsem, 2017). These models allow for obtaining key insights into the role of coupling for the development of LFV in the atmosphere associated with the presence of the ocean.

Recently, dynamical analysis has been conducted by means of the development of a suitable reduced order model of the coupled ocean–atmosphere system. This model has been developed starting from the quasi-geostrophic equations describing the interaction between a two-layer atmosphere and a one-layer ocean over an infinitely deep quiescent ocean layer (Vannitsem et al., 2015; Vannitsem, 2015, 2017; De Cruz et al., 2016, 2018). The ocean flow passively advects the temperature within the ocean, while momentum, radiative, and heat transfer mechanisms realize the coupling between the atmosphere and the ocean. By expanding the solutions of these equations into Fourier series, by truncating them at low wavenumbers, and by projecting onto the Fourier modes retained, a set of ordinary differential equations is derived. The fields are defined over a rectangular domain with $\mathrm{0}\le x\le \mathrm{2}\mathit{\pi }L/n$ and $\mathrm{0}\le y\le \mathit{\pi }L$, with n denoting the aspect ratio between the meridional and the zonal extents of the domain and L the characteristic spatial scale. Moreover, periodic boundaries along the zonal direction and free slip along the meridional direction are chosen for the atmosphere, while a closed basin with no flux through the boundaries is imposed for the ocean.

In the reduced order coupled model version proposed in Vannitsem et al. (2015), a long-periodic attracting orbit combining atmospheric and oceanic variables emerges from a Hopf bifurcation for large values of the meridional gradient of radiative input and frictional coupling. Beyond a certain value of the meridional gradient for the radiative input, a chaotic behavior appears, which is still dominated by LFV on decadal and multi-decadal timescales.

Here we use the original version of the model (Vannitsem et al., 2015) where the four relevant fields, i.e., the barotropic and baroclinic atmospheric streamfunctions, the ocean streamfunction and the ocean temperature, are given by ${\mathit{\psi }}_a={\sum }_{i=\mathrm{1}}^{\mathrm{10}}{\mathit{\psi }}_{\mathrm{a},i}{F}_i$, ${\mathit{\theta }}_a={\sum }_{i=\mathrm{1}}^{\mathrm{10}}{\mathit{\theta }}_{\mathrm{a},i}{F}_i$, ${\mathrm{\Psi }}_o={\sum }_{i=\mathrm{1}}^{\mathrm{8}}{\mathrm{\Psi }}_{\mathrm{o},i}{\mathit{\varphi }}_i$ and $T_o={\sum }_{i=\mathrm{1}}^{\mathrm{8}}{T}_{\mathrm{o},i}{\mathit{\varphi }}_i$, where F_i and ϕ_i are simplified notations for the sets of modes used, compatible with the boundary conditions of both the atmosphere and the ocean. The parameter values used are the ones given in Figs. 8 and 9 of Vannitsem (2017). Depending on the choice of the surface friction coefficient C, different solutions are found with a highly chaotic dynamics without marked LFV in the atmosphere for small values of C, but a more moderately chaotic dynamics with stronger LFV in both the ocean and the atmosphere (related to the development of a coupled mode) for larger values of C.

3 Methods

Traditional multivariate and/or spatiotemporal data analysis methods are commonly based on fixing an orthogonal decomposition basis, satisfying certain mathematical properties such as linearity and/or stationarity (Chatfield, 2016). However, these conditions are not usually met when real-world geophysical data are analyzed, which calls for more adaptive methods (Huang et al., 1998). Indeed, adaptive methods can be helpful for overcoming some limitations of fixed-basis methods, which implicitly assume that a given natural phenomenon or a superposition of physical processes can be represented in terms of a priori defined mathematical functions like sine and/or cosine or some other kinds of wave functions (Chatfield, 2016). Since this cannot be assured, adaptive methods (as the MEMD) could be more suitable for reducing some mathematical assumptions and a priori constraints (Huang et al., 1998; Huang and Wu, 2008; Rehman and Mandic, 2010). Moreover, geophysical data are usually also characterized by scale-invariant features over a wide range of scales with different complexity and show a scale-dependent behavior due to several factors like forcings, coupling, intrinsic variability, and so on (e.g., Lovejoy and Schertzer, 2013; Franzke et al., 2020). For the above reasons, in this work we put forward a novel approach based on combining two different data analysis methods for investigating the multiscale fractal behavior of the coupled ocean–atmosphere system: multivariate empirical mode decomposition (MEMD; Rehman and Mandic, 2010) and generalized fractal dimensions (Hentschel and Procaccia, 1983). One of the main advantages of combining the MEMD with generalized fractal dimensions instead of classical approaches deals with the limited number of intrinsic components that can be also visually inspected. Indeed, if we, for example, use Fourier decomposition we will have a large number of (harmonic) oscillating components at different fixed frequencies that should be summed up for exploiting our proposed procedure. Furthermore, if we, for example, use wavelets we will deal with some a priori assumptions on the decomposition basis onto which we are projecting our data that could produce misleading results in our procedure of evaluating fractal measures on a priori fixed scales. Another advantage is that MEMD allows to preserve some intrinsic properties of signals related to the nonlinear and/or non-stationary nature of processes they are associated with, since the decomposition is based on the local characteristic scale of the data in deriving intrinsic components with time-dependent amplitudes and phases (Huang et al., 1998; Huang and Wu, 2008; Rehman and Mandic, 2010). However, we do not question the appropriateness of conventional analysis techniques but rather acknowledge that other approaches can provide a new perspective on what we can learn from the respective system under study (Alberti et al., 2020a).

3.1 Multivariate empirical mode decomposition (MEMD)

The multivariate empirical mode decomposition (MEMD) is the “natural” multivariate extension of the univariate empirical mode decomposition (EMD) (Huang et al., 1998; Rehman and Mandic, 2010). MEMD directly works on the data domain, instead of defining a conjugate space as for Fourier or wavelet transforms, with the aim of being as adaptive as possible to minimize mathematical assumptions and definitions (Huang et al., 1998) in extracting embedded structures in the form of so-called multivariate intrinsic mode functions (MIMFs) (Rehman and Mandic, 2010). Each MIMF is an oscillatory pattern of the multivariate coordinates having the same number (or differing at most by one) of local extremes and zero crossings, and whose upper and lower envelopes are symmetric (Huang et al., 1998; Rehman and Mandic, 2010). MIMFs are derived through the sifting process (Huang et al., 1998). This process is easily realized for univariate signals (Huang et al., 1998), while it needs to be carefully implemented for multivariate processes (Rehman and Mandic, 2010), since it is based on the cubic spline interpolation of local extremes that cannot easily be defined on a k-dimensional space (Rehman and Mandic, 2010). Rehman and Mandic (2010) proposed an alternative definition of local extremes for multivariate signals by considering the k-variate data as composed by k-dimensional signals projected onto appropriate directions in this k-dimensional space. This allows us to perform cubic spline interpolation in each direction, with the suitable directions chosen by means of a combination of a quasi-Monte Carlo-based low-discrepancy sequences and a uniform angular sampling method (Rehman and Mandic, 2010). These allow providing a more uniform set of direction vectors over which to compute the local mean of envelopes, without introducing any smoother dynamics in the data, via the following procedure:

1. Given a k-dimensional space we need to find the direction vectors by considering that these reduce to points in a (k−1)-dimensional space.

2. The simplest choice is to employ uniform angular sampling on a k-dimensional hypersphere, but this will lead to a non-uniform filling of the k-dimensional space (a higher density of points would be observed near the poles).

3. A quasi-Monte Carlo method is used to provide a more uniform distribution of direction vectors.

4. Once the direction vectors are chosen, the signal is projected onto these vectors, the extrema of the resulting projected signals are evaluated and interpolated component-wise to yield multidimensional envelopes that are then averaged to obtain the multivariate mean.

This means that the quasi-Monte Carlo method is needed only for selecting a uniform sampling of direction vectors, thus avoiding implicitly preferred directions that could be more dominant with respect to the others, which could introduce a source of errors in evaluating signal projections (Rehman and Mandic, 2010).

Having now defined the procedure needed to compute envelopes over each direction, the main steps of the sifting process acting on a k-variate signal $\mathit{s}\left(t\right)=\left[s_{\mathrm{1}}\left(t\right),s_{\mathrm{2}}\left(t\right),\mathrm{\dots },s_k\left(t\right)\right]$ can be summarized as below:

1. identify local extremes (i.e., data points where abrupt changes in the local tendency of the series under study are observed);

2. interpolate local extremes separately by cubic splines (i.e., produce continuous functions with smaller error than other polynomial interpolation);

3. derive the upper and lower envelopes u(t) and l(t), respectively;

4. derive the mean envelope m(t) as $\mathit{m}\left(t\right)=\frac{\mathit{u}\left(t\right)+\mathit{l}\left(t\right)}{\mathrm{2}}$;

5. evaluate the resulting candidate MIMF as $\mathit{h}\left(t\right)=\mathit{s}\left(t\right)-\mathit{m}\left(t\right)$.

The previous steps are iteratively repeated until the obtained candidate MIMF h(t) can be identified as a multivariate intrinsic mode function (also called multivariate empirical mode) (Huang et al., 1998; Rehman and Mandic, 2010), while the full sifting process ends when no more MIMFs c_j(t) can be filtered out from the data. Hence, we can write

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathit{s}\left(t\right)=\sum _{j=\mathrm{1}}^{N_j}{\mathit{c}}_j\left(t\right)+\mathit{r}\left(t\right).\end{array}$$

In this way a multivariate signal is decomposed into N_j k-dimensional functions, each containing the same frequency distribution, e.g., into a set of k-dimensional embedded oscillating patterns c_j(t) which form the multivariate decomposition basis, plus a multivariate residue r(t).

For each MIMF we can define a k^⋆-variate mean timescale as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{\mathit{\tau }}_{j,k^{\star }}={\displaystyle \frac{\mathrm{1}}{T}}\underset{\mathrm{0}}{\overset{T}{\int }}{t}^{\prime }{c}_{j,k^{\star }}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime },\end{array}$$

representing the typical oscillation scale of the jth mode for the k^⋆th univariate component $c_{j,k^{\star }}$ extracted from the multivariate signal $s_{k^{\star }}\left(t\right)$ for $k^{\star }\in [\mathrm{1},k]$. Similarly, by ensemble averaging over the k-dimensional space we can introduce the concept of a multivariate mean timescale as

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{\mathit{\tau }}_j={\displaystyle \frac{\mathrm{1}}{T}}\underset{\mathrm{0}}{\overset{T}{\int }}{t}^{\prime }\langle {\mathit{c}}_j\left(t^{\prime }\right){\rangle }_k\mathrm{d}{t}^{\prime },\end{array}$$

with 〈⋯〉_k denoting an ensemble average over the k-dimensional space. Thus, the k^⋆-variate timescale ${\mathit{\tau }}_{j,k^{\star }}$ is evaluated for each mode and for each k^⋆-dimensional data, while the multivariate mean timescale τ_j is the mean over all $k^{\star }\in (\mathrm{1},k]$. Moreover, as for univariate EMD (Huang et al., 1998), we can introduce the concepts of instantaneous amplitudes a_j(t) and phases ϕ_j(t) of each MEMD mode via the Hilbert transform along the different directions of the k-dimensional space. The instantaneous energy content is then derived as E_j(t)=a_j(t)². Thereby, we can characterize the spectral content by introducing an alternative yet equivalent definition of the power spectral density (PSD) as

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }S\left(\mathit{\tau }\right)={\displaystyle \frac{\mathrm{1}}{T^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{T}{\int }}\langle {\mathbf{E}}_j\left(t^{\prime }\right){\rangle }_k\mathrm{d}{t}^{\prime }\cdot \underset{\mathrm{0}}{\overset{T}{\int }}{t}^{\prime }\langle {\mathit{c}}_j\left(t^{\prime }\right){\rangle }_k\mathrm{d}{t}^{\prime }\doteq {\mathit{\sigma }}^{\mathrm{2}}\left(\mathit{\tau }\right)\cdot \mathit{\tau },\end{array}$$

with σ²(τ) being the k-variate variance of MIMFs and τ the mean timescale defined as in Eq. (6). Moreover, from the instantaneous energy content E_j(t) the relative contribution e_j can be derived as

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{e}_j={\displaystyle \frac{\frac{\mathrm{1}}{T}{\int }_{\mathrm{0}}^T\langle {\mathbf{E}}_j\left(t^{\prime }\right){\rangle }_k\mathrm{d}{t}^{\prime }}{{\sum }_{j=\mathrm{1}}^{N_j}\frac{\mathrm{1}}{T}{\int }_{\mathrm{0}}^T\langle {\mathbf{E}}_j\left(t^{\prime }\right){\rangle }_k\mathrm{d}{t}^{\prime }}}.\end{array}$$

Finally, as for the univariate decomposition (Huang et al., 1998), also the MIMFs are empirically and locally orthogonal with respect to each other, the decomposition basis is a complete set (Rehman and Mandic, 2010), and partial sums of Eq. (4) can be obtained (Alberti, 2018; Alberti et al., 2020b).

3.2 Multivariate and multiscale generalized fractal dimensions

The dynamics of complex systems is usually characterized by a multitude of scales whose dynamical features determine their collective behavior. Nevertheless, vast efforts have been made to determine collective properties of systems (e.g., Hentschel and Procaccia, 1983), instead of considering to measure scale-dependent features. Recently, Alberti et al. (2020a) introduced a new formalism allowing measuring information at different scales by combining a data-adaptive decomposition method and the classical concept of generalized fractal dimensions. The starting point is that a multivariate signal manifesting a multiscale behavior can be written as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }\mathit{s}\left(t\right)=\langle \mathit{s}\rangle +\sum _{\mathit{\tau }}\mathit{\delta }{\mathit{s}}_{\mathit{\tau }}\left(t\right)={\mathit{s}}_{\mathrm{0}}+{\mathit{s}}_{\mathrm{1}}\left(t\right),\end{array}$$

with 〈⋯〉 representing a steady-state average operation and δ indicating a fluctuation at scale τ. For any given τ we can introduce a local natural probability measure dμ_τ such that the probability p_i of visiting the ith hypercube $B_{{\mathit{s}}^{*},\mathit{\tau }}\left(\mathrm{\ell }\right)$ of size ℓ centered at the point s^* on the considered (d-dimensional) phase space of s₁(t) can be defined as

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{p}_i\doteq \underset{{\mathit{s}}_{\mathrm{1}}\in B_{{\mathit{s}}^{*},\mathit{\tau }}\left(\mathrm{\ell }\right)}{\int }\mathrm{d}{\mathit{\mu }}_{\mathit{\tau }}.\end{array}$$

By defining a qth order partition function

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{\mathrm{\Gamma }}_q({\mathit{\mu }}_{\mathit{\tau }},B_{{\mathit{s}}^{*},\mathit{\tau }}(\mathrm{\ell }\left)\right)=\sum _i{p}_i^q=\int \mathrm{d}{\mathit{\mu }}_{\mathit{\tau }}\left(s\right){\mathit{\mu }}_{\mathit{\tau }}\left(B_{{\mathit{s}}^{*},\mathit{\tau }}\right(\mathrm{\ell }){)}^q\end{array}$$

and taking the limit ℓ→0, the multiscale generalized fractal dimensions are derived as

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{D}_{q,\mathit{\tau }}={\displaystyle \frac{\mathrm{1}}{q-\mathrm{1}}}\underset{\mathrm{\ell }\to \mathrm{0}}{lim}{\displaystyle \frac{\log {\mathrm{\Gamma }}_q({\mathit{\mu }}_{\mathit{\tau }},B_{{\mathit{s}}^{*},\mathit{\tau }}(\mathrm{\ell }\left)\right)}{\log \mathrm{\ell }}}.\end{array}$$

Here we identify the intrinsic oscillations by using the MEMD, and then we investigate the phase-space properties at different scales by deriving the generalized dimensions (Alberti et al., 2020a). We summarize this process as follows:

1. We extract multiscale components from s(t) by using the MEMD.

2. We evaluate the intrinsic scale τ_j of each MIMF.

3. We evaluate reconstructions of modes by means of Eq. (4):

   $$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }\sum _{\mathit{\tau }}\mathit{\delta }{\mathit{s}}_{\mathit{\tau }}\left(t\right)\to F_{j^{\star }}\left(t\right)=\sum _{j=\mathrm{1}}^{j^{\star }}{\mathit{c}}_j\left(t\right),\end{array}$$

   with $j^{\star }=\mathrm{1},\mathrm{\dots },N_j$ (by construction, MIMFs are ordered from short to long scales, i.e., ${\mathit{\tau }}_j<{\mathit{\tau }}_{j^{\prime }}$ if $j<j^{\prime }$).

4. We evaluate the generalized dimensions D_q,τ from $F_{j^{\star }}\left(t\right)$ for each j^⋆ (i.e., for each scale ${\mathit{\tau }}_{j^{\star }}$).

5. We evaluate the singularities and singularity spectrum:

   $$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathit{\tau }}={\displaystyle \frac{\mathrm{d}}{\mathrm{d}q}}\left[(q-\mathrm{1})D_{q,\mathit{\tau }}\right]\text{(15)}& {\displaystyle }& {\displaystyle }{f}_{\mathit{\tau }}=f\left({\mathit{\alpha }}_{\mathit{\tau }}\right)=q{\mathit{\alpha }}_{\mathit{\tau }}-\left[(q-\mathrm{1})D_{q,\mathit{\tau }}\right].\end{array}$$

From Eq. (13) we can inspect the local properties of fluctuations in terms of the geometry of the phase space, thus providing a characterization of dynamical features of different regimes and disentangling the different dynamical components of (possibly) different origin.

Figure 1 3-D projection of the full system attractor in the subspace ($T_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}},{\mathit{\psi }}_{\mathrm{a},\mathrm{1}}$) for C=0.008 (red) and C=0.015 (black), respectively.

Our proposed formalism provides a novel way to investigate how phase-space properties (geometry, correlations) change when dynamical components at different mean scales with different dynamics are considered. In other words, we can highlight the role of scale-dependent phenomena in defining the global properties of a system. Indeed, global measures proposed in the past (e.g., Grassberger, 1983; Hentschel and Procaccia, 1983) only allow us to investigate the statistics of the phase-space scaling properties of the whole system; conversely, our proposed approach allows us to investigate how the different scales contribute to the global properties of a system (Alberti et al., 2020a). Moreover, our framework also provides consistency with established measures for characterizing time series from an integral (not scale-resolved) perspective, since the scale-dependent measures we evaluate converge to the associated global measures as all scales are considered, i.e., when the full system dynamics, composed by all accessible scales, is reached (Hentschel and Procaccia, 1983). Within this framework, our approach is promising for investigating scale-dependent properties, as measured by fractal dimensions, of the system. Furthermore, since we are indeed interested in nonlinear variability characteristics at different timescales, employing perfectly linear and/or stationary (harmonic) functions as components would leave out any information on nonlinear dynamics. Moreover, simply looking at the behavior of spectral densities would leave out any higher-order statistical properties, only focusing on the autocorrelation function (i.e., the second-order moment). By looking at the behavior of fractal dimensions we can explore how the different scales contribute to change the phase-space properties for higher-order statistics (i.e., for different values of q).

Figure 2 Spectral energy content S(τ) of the different MIMFs as a function of their mean timescales τ as in Eq. (7) for the full system (atmosphere+ocean, blue circles), only for the atmosphere (orange asterisks), and only for the ocean (yellow diamonds). Panel (a, b) refer to the two values of the friction coefficient, C=0.008 and C=0.015, respectively.

4 Results

4.1 Multivariate empirical mode decomposition

Figure 1 reports the 3-D projection of the full system attractor in the subspace ($T_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}},{\mathit{\psi }}_{\mathrm{a},\mathrm{1}}$) for two representative values of the friction coefficient C (0.008 and 0.015 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ as indicated by red and black points, respectively). In the following, we will omit the physical units of this parameter for the sake of brevity. The considered subspace characterizes the dynamics of the system as represented by the dominant mode of the meridional temperature gradient in the ocean (T_o,2), by the double-gyre transport within the ocean (Ψ_o,2), and by the vertically averaged zonal flow within the atmosphere (ψ_a,1), respectively.

The behavior of the system is clearly dependent on the friction coefficient, with both the location and the topology of the attractor changing as C is increased from 0.008 (red points in Fig. 1) to 0.015 (black points in Fig. 1). This behavior has also been previously reported by Vannitsem et al. (2015) and Vannitsem (2015), indicating a drastic qualitative change of the nature of the dynamics at about C=0.011 above which substantial LFV emerges (Vannitsem et al., 2015; Vannitsem, 2015, 2017). However, all model components are clearly characterized by multiscale variability, spanning a wide range of timescales that can contribute to the dynamics in different ways, depending on the values of the friction coefficient and the intrinsic variability of the coupled ocean–atmosphere system.

Figure 2 displays the behavior of the spectral energy content S(τ) of the different MIMFs as a function of their mean timescales τ as in Eq. (7) for the full system (atmosphere+ocean) and for the two subsystems separately (i.e., the atmosphere and the ocean, respectively).

First of all, it is important to underline that a different number of MIMFs has been identified for the two different cases: N_j=17 for C=0.008 and N_j=22 for C=0.015. This underlines that the respective dynamical behavior of the system is different, being characterized by different sets of empirical modes and consequently by a different number of relevant timescales. Moreover, by keeping in mind that for pure noise the expected number of MIMFs is log ₂N with N being the number of data points, both situations cannot be related to a purely stochastic dynamics. Indeed, in both cases we have used N=10⁵ data points; thus the expected number of MIMFs is $N_j^{\text{noise}}=\mathrm{16}$ (Flandrin et al., 2004). However, an interesting feature is that for the lower C value a number of MIMFs closer to that expected for noisy data are found, possibly related to the more irregular dynamics in this low-friction coefficient case. Conversely, a marked departure from N_j=16 is found for the higher C case, corresponding to a more regular dynamics characterized by significant LFV.

Furthermore, from Fig. 2 it is easy to note that the behavior of S(τ) depends on both the friction coefficient C and the different components of the model. For the full system (i.e., atmosphere+ocean) S(τ) decreases as τ increases for both values of C, while it is characterized by increasing spectral energy content at larger scales (i.e., at lower frequencies). By discriminating between the atmospheric and the oceanic contribution we are able to see that (as expected), the short-term variability of the full system can be attributed to the atmosphere, while the long-term one is a reflection of the ocean dynamics. Moreover, when C increases we note an increase of the spectral energy content at all timescales, together with a flattening of the atmospheric spectral behavior, while the ocean dynamics seems to preserve its spectral features. These behaviors can be related to the existence of multiscale variability of the full system that can be linked to the different components operating at different timescales and to the different dynamics of the system as the friction coefficient C is changed.

Figure 3 Relative contribution (in percentage) E_χ,τ of each variable $\mathit{\chi }=\mathit{\{}{\mathit{\psi }}_{\mathrm{a},i},{\mathit{\theta }}_{\mathrm{a},i},{\mathrm{\Psi }}_{\mathrm{o},i},T_{\mathrm{o},i}\mathit{\}}$ in dependence on the mean timescale τ. Panels (a) and (b) refer to the two values of the friction coefficient C=0.008 and C=0.015, respectively. The white line separates the atmospheric variables from the oceanic ones.

Figure 4 3-D projection of the full system attractor in the subspace ($T_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}},{\mathit{\psi }}_{\mathrm{a},\mathrm{1}}$) for C=0.008 (red) and C=0.015 (black), respectively, as obtained from reconstructions based on the multivariate empirical modes R_χ,95 %(t) accounting for 95 % of the total variance of the model dynamics.

To further clarify the latter aspect, we evaluate the relative contribution (in percentage) E_χ,τ of the different MIMFs (i.e., at different timescales τ) for each variable $\mathit{\chi }=\mathit{\{}{\mathit{\psi }}_{\mathrm{a},i},{\mathit{\theta }}_{\mathrm{a},i},{\mathrm{\Psi }}_{\mathrm{o},i},T_{\mathrm{o},i}\mathit{\}}$ as reported in Fig. 3. It can be clearly noted that the oceanic variability mainly contributes to the low-frequency dynamics ($E_{\mathit{\chi },\mathit{\tau }}>\mathrm{95}\mathrm{\%}$ for $\mathit{\chi }=\mathit{\{}{\mathrm{\Psi }}_{\mathrm{o},i}$, T_o,i} and τ≳10⁴ d), while the atmosphere is mainly characterized by short-term variability for C=0.008 ($E_{\mathit{\chi },\mathit{\tau }}>\mathrm{95}\mathrm{\%}$ for $\mathit{\chi }=\mathit{\{}{\mathit{\psi }}_{\mathrm{a},i}$, θ_a,i} and τ≲10 d) and by both short- and long-term dynamics for C=0.015. This points towards the C-dependent behavior of the atmospheric dynamics, with the ocean multiscale variability being less affected by changes in the values of the friction coefficient, and to the role of the ocean in developing LFV in the atmosphere as C increases.

Thanks to the completeness property of the MEMD we can explore the dynamics of the system as reproduced by the most energetic empirical modes via partial sums of Eq. (4). By using the information coming from the energy percentage distribution across the different timescales for each variable χ, we can provide MIMF reconstructions accounting for a certain percentage of energy with respect to the total spectral energy content. By ordering the empirical modes with decreasing relative contribution e_j and summing up those contributing at least 95 % of the total spectral content, we are able to investigate the 3-D projection of the full system attractor onto the subspace ($T_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}},{\mathit{\psi }}_{\mathrm{a},\mathrm{1}}$) and compare it with the projection obtained by considering all timescales (as in Fig. 1). Thus, for each variable $\mathit{\chi }=\mathit{\{}{\mathit{\psi }}_{\mathrm{a},i},{\mathit{\theta }}_{\mathrm{a},i},{\mathrm{\Psi }}_{\mathrm{o},i},T_{\mathrm{o},i}\mathit{\}}$, we can define a reconstruction based on empirical modes, R_χ,95 %, as

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{R}_{\mathit{\chi },\mathrm{95}\mathrm{\%}}\left(t\right)\doteq \sum _{j^{\prime }|e_{j^{\prime }}\ge \mathrm{95}\mathrm{\%}}{\mathit{c}}_{\mathit{\chi },j^{\prime }}\left(t\right),\end{array}$$

with ${\mathit{c}}_{\mathit{\chi },j^{\prime }}\left(t\right)$ being the j^′th multivariate empirical mode extracted via the MEMD of the variables χ. The 3-D projection onto the subspace ($T_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}},{\mathit{\psi }}_{\mathrm{a},\mathrm{1}}$) of R_χ,95 % is shown in Fig. 4, while Table 1 summarizes the mode indices j^′ and corresponding k^⋆-variate timescales ${\mathit{\tau }}_{j^{\prime },k^{\star }}$ (see Eq. 5) used for the reconstruction.

Figure 5 Multiscale correlation dimension D_2,τ for C=0.008 at different timescales τ_j for different cases: (a) for each MIMF individually ($D_{\mathrm{2}}^j$), (b) for reconstructions of MIMFs as in Eq. (12) (D_2,τ), and (c) for reconstructions of MIMFs separately for each variable (barotropic modes – blue circles, baroclinic modes – orange asterisks, transport modes – yellow diamonds, and temperature modes – violet symbol). Each panel also shows the 95 % confidence intervals as error bars.

Table 1 Mode indices j^′ and corresponding k^⋆-variate timescales ${\mathit{\tau }}_{j^{\prime },k^{\star }}$ (see Eq. 5) used for the reconstruction based on empirical modes R_χ,95 %.

By comparing Figs. 1 and 4 it can be easily noted that the underlying structure of the 3-D projection of the full attractor is essentially the same, thus suggesting that the subspace statistics of the phase-space scaling information can be recovered by a subset of multivariate empirical modes. This underlines that the dynamics of the full system can be reproduced by only few relevant timescales without too much loss of information, thus reducing the complexity of the low-order model itself. These results appear relevant if put into the wider context of coupled ocean–atmosphere dynamics, allowing us to recover the main features by only considering the most relevant (in terms of energy) timescale dynamical components.

4.2 Multiscale generalized fractal dimensions

Under general conditions, the complexity of a dynamical system can be conveniently investigated by means of the nonlinear properties of its phase-space trajectory (e.g., its attractor or repeller in case of dissipative dynamics) (Ott, 2002). One of the most common ways to characterize the topology of an attractor is to compute its spectrum of generalized fractal dimensions, allowing us to statistically characterize important properties of the dynamics as reflected by its phase-space geometry, including its information content, complexity, and underlying fractal structure (Grassberger, 1983; Hentschel and Procaccia, 1983; Donner et al., 2011). However, classical approaches can only provide global information on the phase-space topology (Hentschel and Procaccia, 1983; Ott, 2002), while multiscale dynamical systems can be characterized by the statistics of the phase-space scaling changing as different real-space scales are considered (Alberti et al., 2020a). For this purpose, we investigate the statistics of the phase-space scaling of the coupled ocean–atmosphere model by evaluating the multiscale generalized fractal dimensions described in Section 3.2. Figures 5 and 6 report the behavior of the correlation dimension D₂ for both values of the friction coefficient and for three different cases: (a) for each MIMF individually ($D_{\mathrm{2}}^j$), (b) for reconstructions of MIMFs (D_2,τ), and (c) for reconstructions of MIMFs performed separately for each variable $\mathit{\chi }=\mathit{\{}{\mathit{\psi }}_{\mathrm{a},i},{\mathit{\theta }}_{\mathrm{a},i},{\mathrm{\Psi }}_{\mathrm{o},i},T_{\mathrm{o},i}\mathit{\}}$.

Figure 6 Same as in Fig. 5, but for C=0.015.

As expected, the multiscale correlation dimension for each MIMF decreases with increasing timescale, being representative of a more regular, less stochastic/chaotic, behavior of large-scale MIMFs as compared with the short-term ones (Alberti et al., 2020a). Particularly, when approaching the largest timescales, $D_{\mathrm{2},\mathit{\tau }}\to \mathrm{1}$ suggesting the existence of fixed-scale MIMFs, i.e., with the instantaneous frequencies being almost constant (as expected, e.g., Rehman and Mandic, 2010). Conversely, when the multiscale correlation dimensions are evaluated by summing up the different MIMFs, starting from the shortest up to the largest scale, a clearly scale-independent behavior of D_2,τ is highlighted for both values of the friction coefficient C. This suggests that the short-term variability mostly defines the correlations between pairs of points in the phase space, thus setting the minimum number of variables needed to describe the dynamics of the system, i.e., its degrees of freedom. However, the role of C clearly emerges in determining the values of D_2,τ, being lower for the larger C value. Indeed, $D_{\mathrm{2},\mathit{\tau }}\sim \mathrm{8}$ for C=0.008, while $D_{\mathrm{2},\mathit{\tau }}\sim \mathrm{1.5}$ for C=0.015. This reflects the different statistics of the attractor scaling of the full system associated with a different dynamical behavior of the model variables (Faranda et al., 2019), also suggesting a less chaotic nature of the system as C increases, together with a reduced number of degrees of freedom. This points towards the possibility of recovering the main features of the model with a reduced number of variables and scales. However, the most interesting features emerge when the different variables of both atmosphere and ocean are separately investigated by means of the multiscale generalized fractal dimensions. It is indeed evident that a scale-independent behavior is found for the atmosphere for both values of C, while a scale-dependent behavior is observed for the ocean. The former can be easily related to the dominant role of the short-term variability for the atmosphere, while the latter is a reflection of the long-term dynamics of the ocean. Moreover, it is also particularly interesting to note that higher (lower) D_2,τ values are found for the atmosphere with respect to the ocean for C=0.008 (C=0.015). This reflects the role of the ocean in developing LFV in the atmosphere as C increases, although the complexity of the full system seems to be determined by the atmosphere for both C values, being indeed characterized by a scale-independent behavior of D_2,τ.

Figure 7 Multiscale capacity dimension D_0,τ, multiscale information dimension D_1,τ, and multiscale correlation dimension D_2,τ for C=0.008 at different timescales τ_j for different cases: (a) for each MIMF individually ($D_q^j$), (b) for reconstructions of MIMFs as in Eq. (12) (D_q,τ), and (c–f) for reconstructions of MIMFs separately for each variable (barotropic modes – c, baroclinic modes – d, transport modes – e, and temperature modes – f).

The described findings are not only valid for the multiscale correlation dimension D_2,τ but are also observed for both the multiscale capacity dimension D_0,τ and the multiscale information dimension D_1,τ as reported in Figs. 7 and 8, together with the multiscale correlation dimension D_2,τ, for both values of C.

Our formalism reveals the expected property that for $q<q^{\prime }$, $D_{q,\mathit{\tau }}>D_{q^{\prime },\mathit{\tau }}\forall \mathit{\tau }$ (Hentschel and Procaccia, 1983; Alberti et al., 2020a). Moreover, when evaluating the multiscale generalized fractal dimensions for each MIMF separately (e.g., Figs. 7a and 8a) a decreasing value for $D_q^j$ is found as τ increases, with all $D_q^j$ converging towards the same value of 1 at large timescales. As for $D_{\mathrm{2}}^j$ this behavior can be easily interpreted in terms of more chaotic vs. more regular MIMFs when moving from short to large scales. This indeed reflects the existence of large-scale MIMFs that are characterized by a linear phase, i.e., a constant timescale (e.g., Rehman and Mandic, 2010). Thus, this is a trivial result. Conversely, when the D_q,τ are evaluated for reconstructions based on MIMFs a scale-independent behavior is found for the full system for both values of C (e.g., Figs. 7b and 8b). However, the key role of the friction coefficient clearly emerges by looking at the larger values of D_q,τ for C=0.008 with respect to the lower values found for C=0.015. This clearly indicates the existence of a completely different dynamics between the two values of C, where the coupled ocean–atmosphere dynamics can be interpreted as a higher-dimensional chaotic system for reduced ocean–atmosphere coupling (i.e., C=0.008) as opposed to a lower-dimensional one for a strong ocean–atmosphere coupling (i.e., C=0.015). Although C acts as a control parameter for the dimensionality of the system, it is not able to change the underlying fractal nature of the full system. Indeed, for both C values we clearly observe different D_q,τ for different q, thus suggesting the existence of a multifractal nature of the ocean–atmosphere dynamics at all timescales. Furthermore, by separately looking at the two subsystems (i.e., the ocean and the atmosphere) a completely different behavior emerges (e.g., Figs. 7c–f and 8c–f). In this case, the atmospheric variables are characterized by scale-independent D_q,τ, being representative of a high-dimensional system whose prime dynamics occurs at short timescales and with little effects of large-scale processes on the collective dynamics of the atmosphere. By contrast, a clearly scale-dependent behavior is found for the oceanic variables, with the multiscale generalized dimensions decreasing at larger timescales, reflecting the effects of large-scale dynamics dominating with respect to the short-term one for the ocean variability. Again the friction coefficient C controls the values of D_q,τ, decreasing as C increases, while both the atmosphere and the ocean are clearly characterized by multifractal features at all timescales.

Figure 8 Same as in Fig. 7, but for C=0.015.

By estimating the Lyapunov spectra (cf. Fig. S11 in the Supplement) separately for the ocean and the atmosphere, we obtained that for C=0.008 the instability is large for the atmosphere with a Lyapunov dimension D_L∼10, while for C=0.015 the instability is weaker for the atmosphere, and the Lyapunov dimension is slightly larger than 4. Following the Kaplan–Yorke conjecture (Kaplan and Yorke, 1979), the Lyapunov dimension can be used as a proxy of D₀. Hence, our results are clearly consistent with the dimension estimates for the atmosphere. For the ocean, however, there seems to be a less good agreement, with D_L≈2 while we found that $D_{\mathrm{0},\mathit{\tau }}\approx \mathrm{4}$. This quantitative disagreement could be related to the fact that the ocean can be viewed as a relatively stable system perturbed by high-frequency “noise” due to the atmosphere. Deeper investigations will be devoted to clarify this point in future research.

Figure 9 D_q,τ spectra for the coupled ocean–atmosphere dynamics at different timescales τ_j (indicated by different line colors) for reconstructions of MIMFs as in Eq. (12) (D_q,τ) for (a) C=0.008 and (b) C=0.015.

Figure 10 Multifractal width Δ at different timescales τ_j for reconstructions of MIMFs as in Eq. (12) (D_q,τ) for (a) C=0.008 and (b) C=0.015. The different colors refer to the full system (atmosphere+ocean, black circles), only the atmosphere (red circles), and only the ocean (green circles), respectively.

As a further step, we evaluate the full spectrum of generalized fractal dimensions for each MIMF by considering a wide range of statistical moments q. As suggested in Lovejoy and Schertzer (2013) the range of significant moments can be evaluated by means of the tail of the cumulative distribution function of the data. Indeed, the effect of sample size and its implications for spurious scaling may be due either to first- or second-order multifractal phase transitions (Lovejoy and Schertzer, 2013). To mitigate these effects (e.g. Supplement) since we deal with the investigation of scale-dependent fractal dimensions, we evaluate the cumulative statistics at different scales and we observe that extreme fluctuations follow a power law decay leading to the divergence of the sixth-order and the fourth-order moment for C=0.008 and C=0.015, respectively. Thus we fix our range of moments $-\mathrm{6}<q<\mathrm{6}$ and $-\mathrm{4}<q<\mathrm{4}$ for C=0.008 and C=0.015, respectively. This analysis allows characterizing how the (multi)fractal properties of the system evolve with the timescale τ. Indeed, there are ongoing discussions on the fractal structure of both the atmosphere and the ocean, especially dealing with the short-term variability and in terms of scaling-law behavior and statistics of increments (e.g., Lovejoy and Schertzer, 2013; Franzke et al., 2020).

The D_q,τ spectrum is reported in Fig. 9, where colored lines correspond to different timescales. It can be observed that for both values of the friction coefficient C, different values of D_q,τ are obtained for different q, with being D_q,τ a nonlinear decreasing function of q. This means that the full system exhibits signatures of multifractality at all timescales, especially at very short and very long timescales. A simple and direct measure of the degree of multifractality¹ is the so-called multifractal width $\mathrm{\Delta }\doteq D_{q_{\text{min}},\mathit{\tau }}-D_{q_{\text{max}},\mathit{\tau }}$. We observe (see Fig. 10a, b, black circles) that Δ≈2 for $\mathit{\tau }\in [{\mathit{\tau }}_S,{\mathit{\tau }}_L]$ days, while Δ>2 for both τ<τ_S and τ>τ_L, with τ_S∼20 d and τ_L∼1 year. This behavior could be the reflection of processes operating at different timescales for both the atmosphere (at short timescales) and the ocean (at long timescales). In order to further disentangle those processes, we also evaluated the full spectra of the generalized multifractal dimensions for each subsystem (i.e., atmosphere and ocean) individually. For both values of C, the corresponding results are shown in Fig. 11.

Figure 11 D_q,τ spectra for the dynamics of atmosphere and ocean individually at different timescales τ_j (indicated by different line colors) for reconstructions of MIMFs as in Eq. (12) (D_q,τ) for (a, c) C=0.008 and (b, d) C=0.015. Panels (a, b) refer to the atmosphere, (c, d) to the ocean.

We clearly see that for the atmosphere, there is a scale-independent behavior of D_q,τ for all q, rendering the different curves almost invariant with respect to the scale. By contrast, a scale-dependent behavior emerges for the ocean for the lower value of C. Indeed, it is evident that as the timescale increases the multiscale generalized dimensions tend to decrease for all values of q, moving from $D_{q,{\mathit{\tau }}_{\mathrm{1}}}\in [\mathrm{5},\mathrm{8}]$ to $D_{q,{\mathit{\tau }}_{\mathrm{17}}}\in [\mathrm{2},\mathrm{3}]$ for C=0.008. Conversely, although there is an overall reduction in the D_q,τ values for C=0.0015 with respect to those evaluated for C=0.008, the decrease with the timescale is less evident for this higher C value, although it is still present for τ>1 year (see orange and red curves in comparison with the blue ones in Fig. 11d). This clearly suggests that the presence of strong multifractality in the full system can be essentially attributed to the atmosphere, with only a marginal role of the ocean variability in determining the fractal structure of the full system. By evaluating the difference between $D_{q_{\text{min}},\mathit{\tau }}$ and $D_{q_{\text{max}},\mathit{\tau }}$, we can clearly see that larger values, of the order of 3, are found for the atmosphere, at almost all timescales (and especially at shorter timescales), for both values of C. Conversely, larger values are found at shorter timescales for both values of C for the ocean. As the timescale increases, this difference tends to be reduced to values close to 1, suggesting a reduced multifractality of the ocean with respect to the atmosphere, especially for the lower value of C at larger timescales when the role of the ocean becomes dominant as compared to the atmosphere (see Fig. 2).

4.3 Comparison with regional averages from reanalysis data

As a final step we compare our previous results for the reduced order coupled ocean–atmosphere model with those obtained from reanalysis data (Poli, 2015). More specifically, we use three different sets of regional time series based on the European Centre for Medium-range Weather Forecasts (ECMWF) ORA-20C project (De Boisséson and Balmaseda, 2016; De Boisséson et al., 2017) that is a 10-member ensemble of ocean reanalyses covering the complete 20th century using atmospheric forcing from the ERA-20C reanalysis (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-20c). Here, we focus on data from January 1958 to December 2009 at monthly resolution in terms of different monthly averaged time series, the set of data also used previously in Vannitsem and Ekelmans (2018). This period has been chosen in the latter study because of the ocean reanalysis dataset showing here smaller uncertainties than during the first half of the 20th century (De Boisséson and Balmaseda, 2016).

Three different representative regions are chosen: the North Atlantic region, corresponding to the domain defined by λ ∈ [55, 15^∘ W] and ϕ ∈ [25, 60^∘ N], the North Pacific region, i.e., a spherical–rectangle domain with λ ∈ [165, 225^∘ E] and ϕ ∈ [25, 60^∘ N], and the tropical Pacific region, corresponding to λ ∈ [165, 225^∘ E] and ϕ ∈ [25^∘ S, 25^∘ N] (Vannitsem and Ekelmans, 2018). The individual series for the two extratropical regions have been derived by projecting the reanalysis fields on two dominant Fourier modes: (i) $F_{\mathrm{1}}=\sqrt{\mathrm{2}}\cos \left(\mathit{\pi }y/L_y\right)$ and (ii) ${\mathit{\varphi }}_{\mathrm{2}}=\mathrm{2}\sin \left(\mathit{\pi }x/L_x\right)\sin \left(\mathrm{2}\mathit{\pi }y/L_y\right)$ (Vannitsem and Ekelmans, 2018). For the tropical Pacific region, the series are formed by spatial averages. In this way, we obtain two sets of three time series each for both the North Atlantic and the North Pacific (i.e., one for the atmosphere and two for the ocean), as well as a third set of three time series for the tropical Pacific (two for the atmosphere at two different pressure levels and one for the ocean). This allows us to build a 3-D projection of the local atmosphere–ocean coupled dynamics for each region (see Vannitsem and Ekelmans, 2018, for more details).

Figure 12 Spectral energy content S(τ) of the different MIMFs as a function of their mean timescales τ as in Eq. (7) for the North Atlantic (blue circles), the North Pacific (orange asterisks), and the tropical Pacific (yellow diamonds).

By using the MEMD analysis to investigate the multivariate patterns of reanalysis data, we found the same number of N_j=9 MIMFs for each region, whose mean timescales range from ∼2 months up to ∼20 years, suggesting the existence of multiscale variability over a wide range of scales. As for the reduced order model, we first investigate the behavior of the spectral energy content S(τ) of the different MIMFs as a function of their mean timescales τ as in Eq. (7) for the three different regions as shown in Fig. 12. We clearly observe an increase of the spectral energy content up to a timescale τ∼1 year for all regions, then declining for both the North Atlantic and the North Pacific. Conversely, the tropical Pacific is characterized by larger spectral content also for timescales larger than 1 year, up to τ∼5 years, which coincide with the typical timescales of the El Niño–Southern Oscillation (ENSO). Furthermore, for all regions a decreasing spectral energy content is found at the largest timescales (i.e., τ>5 years).

To further compare our above model results with those obtained for the reanalysis data, we evaluate the multiscale generalized fractal dimensions for the three different regions. For each region, we derive both the multifractal width $\mathrm{\Delta }\doteq D_{q_{\text{min}},\mathit{\tau }}-D_{q_{\text{max}},\mathit{\tau }}$ and the full multiscale multifractal spectrum at different timescales τ_j for reconstructions of MIMFs as in Eq. (12) (D_q,τ). Figure 13 shows the corresponding results for the North Atlantic region, the North Pacific region, and the tropical Pacific region.

Figure 13 (a–c) Multifractal width Δ and (d–f) D_q,τ spectra at different timescales τ_j for reconstructions of MIMFs as in Eq. (12) (D_q,τ) for (a, d) the North Atlantic, (b, e) the North Pacific, and (c, f) the tropical Pacific, respectively.

First of all, it is important to underline that the multiscale generalized fractal dimensions are clearly different with respect to those obtained from the ocean–atmosphere model. This directly follows from the different numbers of variables (time series) in the model, being a 36-dimensional dynamical system, with respect to the reanalysis data, being a 3-dimensional projection of the regional ocean–atmosphere dynamics. Nevertheless, although different in terms of absolute values, both the model and the reanalysis data show a similar qualitative behavior with varying scale τ, although some differences are found between the different regions.

On one hand, both the North Atlantic and the North Pacific regions (see Fig. 13d, e) are characterized by a scale-dependent behavior, with decreasing D_q,τ as τ increases. Moreover, by looking at the multifractal width as a function of the scale (Fig. 13a, b) we find evidence for a decreasing Δ as τ increases, being representative of a transition from a short-term multifractal nature to a long-term monofractal one. These features can be interpreted in terms of the different multiscale dynamical processes affecting the atmosphere on short scales and the ocean on larger scales.

On the other hand, by looking at the tropical Pacific region we clearly see an enhancement of Δ, i.e., the emergence of multifractal features (see Fig. 13c), at annual/multi-annual timescales (i.e., τ∼1–8 years), being also characterized by the largest values of the multiscale generalized fractal dimensions (see Fig. 13f). This could be related to the role of the El Niño–Southern Oscillation (ENSO) cycle manifesting at these timescales (between 2 and 7 years), which is likely responsible for the different scale-dependent behavior of D_q,τ as compared to the two other extratropical regions.

In summary, by means of the reanalysis data, we have been able to demonstrate that (i) the reduced order coupled ocean–atmosphere model and the reanalysis data show some qualitatively similar behavior of the multiscale generalized fractal dimensions, although they are characterized by different absolute values due to the different numbers of variables considered in the model and the projections on a few modes of the reanalysis data, and that (ii) interesting features emerge when looking at the scale dependency of the statistics of the phase-space scaling for different regions, being the reflection of different driving mechanisms and processes operating at different timescales in the coupled ocean–atmosphere system. However, further investigations are needed to characterize the role of the different processes as well as their intrinsic dimensionality, occurrence, and spatial dependency in more detail. Such an in-depth investigation is outlined as a part of our future work.

5 Conclusions

We have provided a first time systematic investigation of the multiscale dynamics of a reduced order coupled ocean–atmosphere model (Vannitsem et al., 2015) as described by means of the statistics of the phase-space scaling (Alberti et al., 2020a).

First, by means of the multivariate empirical mode decomposition (MEMD) we have been able to detect oscillating patterns with time-dependent amplitude and frequency that are directly linked to a rich variety of features of the coupled ocean–atmosphere system. We have found that the underlying structure of the 3-D projection of the full attractor is essentially reproduced by a subset of multivariate intrinsic mode functions (MIMFs) corresponding to the most relevant timescales without too much loss of information, thus further reducing the complexity of the reduced order model itself. These results appear relevant if put into the wider context of coupled ocean–atmosphere dynamics, allowing us to recover the main features by only considering the most relevant (in terms of energy) timescale dynamical components.

Second, by exploiting the novel concept of multiscale and multivariate generalized fractal dimensions we have investigated the different multifractal properties for the ocean and the atmosphere at different timescales. We have demonstrated that for weak ocean–atmosphere coupling (i.e., for low values of the friction coefficient C), the resulting dimensions of the two model components are very different, while for strong coupling (larger C) at which coupled modes develop at low frequencies, the scaling properties are more similar especially at longer timescales. These results suggest that as C increases, we observe the development of a coherent coupled dynamics, primarily at large timescales. In terms of the underlying fractal structure, we have found that for both considered values of the friction coefficient C, the full system exhibits signatures of multifractality at all timescales, especially pronounced at short and long as compared to intermediate timescales. By means of the full spectrum of generalized fractal dimensions, we have clearly evidenced that for the atmosphere, there is a scale-independent behavior of D_q,τ for all q, rendering the multifractal spectra almost invariant with respect to the timescale. By contrast, a scale-dependent behavior emerges for the ocean for the lower value of C. This clearly suggests that the presence of strong multifractality in the full system can be attributed to the atmosphere, with only a marginal role of the ocean variability in determining the fractal structure of the full system.

Finally, we have compared our results for the reduced order coupled ocean–atmosphere model with those derived from reanalysis data (Poli, 2015) by using three sets of different regional time series from the ORA-20C project (De Boisséson and Balmaseda, 2016; De Boisséson et al., 2017). Although the resulting multiscale generalized fractal dimensions clearly differ quantitatively from those obtained from the ocean–atmosphere model – which can be easily understood by considering the different dimensions of the model (a 36-dimensional dynamical system) and the reanalysis data (3-dimensional projections of the local ocean–atmosphere dynamics) – we observed a similar qualitative behavior with changing scale τ. Interestingly, the multiscale multifractal features of different regions show different scale-dependent behaviors. Specifically, both the North Atlantic and the North Pacific regions are characterized by a scale-dependent behavior, with decreasing D_q,τ as τ increases, with a transition from a short-term multifractal nature to long-term monofractal one. These features can be interpreted in terms of the different multiscale dynamical processes affecting the atmosphere at short timescales and the ocean at longer timescales. Conversely, the tropical Pacific region is characterized by the emergence of multifractal features at annual/multi-annual timescales (i.e., τ∼1–8 years), being also characterized by the largest values of the multiscale generalized fractal dimensions. This behavior can be seen as a manifestation of the El Niño–Southern Oscillation (ENSO) cycle that typically acts at these timescales and can be considered the key driving factor of a different scale-dependent behavior of D_q,τ as compared to the two extratropical regions.

Our findings for both the model and the reanalysis data suggest that our approach can be used to diagnose the strength of coupling in the ocean–atmosphere system and to investigate the statistics of the phase-space scaling of the system. We have demonstrated that the model and the reanalysis data show a qualitatively similar behavior of the multiscale generalized fractal dimensions. However, the different scale dependency of the statistics of the phase-space scaling for different regions can contribute to a better understanding of the different driving mechanisms and processes operating at different timescales in the coupled ocean–atmosphere system. Indeed, our results highlight that the complexity of the coupled ocean–atmosphere system significantly depends not only on model parameters, which can be helpful for reproducing different features of the dynamics, but also on the particular scale we are looking at that can be related to different phenomena and source mechanisms, of both intrinsic and external origin to the ocean–atmosphere system. This means that our results could also be helpful for understanding the dimensionality of the system at different timescales, thus being useful for forecasting the dynamics at different scales and for building empirical models based on dynamical system approaches in a similar fashion to models developed considering real space scaling behavior (e.g., Del Rio Amador and Lovejoy, 2019, 2021). These observations suggest that further investigations are needed to better characterize the role of the different processes as well as their intrinsic dimensionality, occurrence, and spatial dependency, which shall be further addressed in our future work.