the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Multiscale fractal dimension analysis of a reduced order model of coupled ocean–atmosphere dynamics
Tommaso Alberti
Reik V. Donner
Stéphane Vannitsem
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 timedependent 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 phasespace 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.
 Article
(3551 KB) 
Supplement
(1711 KB)  BibTeX
 EndNote
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/multiannual 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 lowfrequency 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 36variable 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 longrange 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 boxcounting 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 boxcounting (or capacity or simply fractal) dimension
the information dimension
and the correlation dimension
with Θ(⋯) being the Heaviside function. More specifically, D_{0} is a measure of the sparseness of the phase space by the studied system's dynamics, D_{1} is a measure of the information gained on the phase space with a given accuracy, while D_{2} is a measure of correlations, i.e., mutual dependence, between phasespace 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 phasespace 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 phasespace 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 stateoftheart 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 phasespace 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 phasespace trajectory that can be used for reconstructing the skeleton of the phasespace 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 phasespace 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 phasespace scaling.
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 quasigeostrophic equations describing the interaction between a twolayer atmosphere and a onelayer 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 longperiodic 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 multidecadal 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.
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 realworld 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 fixedbasis 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 scaleinvariant features over a wide range of scales with different complexity and show a scaledependent 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 nonstationary 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 timedependent 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 socalled 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 kdimensional space (Rehman and Mandic, 2010). Rehman and Mandic (2010) proposed an alternative definition of local extremes for multivariate signals by considering the kvariate data as composed by kdimensional signals projected onto appropriate directions in this kdimensional space. This allows us to perform cubic spline interpolation in each direction, with the suitable directions chosen by means of a combination of a quasiMonte Carlobased lowdiscrepancy 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:

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

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

A quasiMonte Carlo method is used to provide a more uniform distribution of direction vectors.

Once the direction vectors are chosen, the signal is projected onto these vectors, the extrema of the resulting projected signals are evaluated and interpolated componentwise to yield multidimensional envelopes that are then averaged to obtain the multivariate mean.
This means that the quasiMonte 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 kvariate 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:

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

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

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

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}}$;

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
In this way a multivariate signal is decomposed into N_{j} kdimensional functions, each containing the same frequency distribution, e.g., into a set of kdimensional 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
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 kdimensional space we can introduce the concept of a multivariate mean timescale as
with 〈⋯〉_{k} denoting an ensemble average over the kdimensional 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 kdimensional space. The instantaneous energy content is then derived as E_{j}(t)=a_{j}(t)^{2}. Thereby, we can characterize the spectral content by introducing an alternative yet equivalent definition of the power spectral density (PSD) as
with σ^{2}(τ) being the kvariate 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
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 scaledependent features. Recently, Alberti et al. (2020a) introduced a new formalism allowing measuring information at different scales by combining a dataadaptive 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
with 〈⋯〉 representing a steadystate 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 (ddimensional) phase space of s_{1}(t) can be defined as
By defining a qth order partition function
and taking the limit ℓ→0, the multiscale generalized fractal dimensions are derived as
Here we identify the intrinsic oscillations by using the MEMD, and then we investigate the phasespace properties at different scales by deriving the generalized dimensions (Alberti et al., 2020a). We summarize this process as follows:

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

We evaluate the intrinsic scale τ_{j} of each MIMF.

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}$).

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}}$).

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.
Our proposed formalism provides a novel way to investigate how phasespace 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 scaledependent 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 phasespace 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 scaleresolved) perspective, since the scaledependent 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 scaledependent 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 higherorder statistical properties, only focusing on the autocorrelation function (i.e., the secondorder moment). By looking at the behavior of fractal dimensions we can explore how the different scales contribute to change the phasespace properties for higherorder statistics (i.e., for different values of q).
4.1 Multivariate empirical mode decomposition
Figure 1 reports the 3D 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}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\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 doublegyre 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 _{2}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^{5} 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 lowfriction 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 shortterm variability of the full system can be attributed to the atmosphere, while the longterm 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.
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 lowfrequency dynamics (${E}_{\mathit{\chi},\mathit{\tau}}>\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathrm{\%}$ for $\mathit{\chi}=\mathit{\{}{\mathrm{\Psi}}_{\mathrm{o},i}$, T_{o,i}} and τ≳10^{4} d), while the atmosphere is mainly characterized by shortterm variability for C=0.008 (${E}_{\mathit{\chi},\mathit{\tau}}>\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathrm{\%}$ for $\mathit{\chi}=\mathit{\{}{\mathit{\psi}}_{\mathrm{a},i}$, θ_{a,i}} and τ≲10 d) and by both short and longterm dynamics for C=0.015. This points towards the Cdependent 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 3D 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
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 3D 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.
By comparing Figs. 1 and 4 it can be easily noted that the underlying structure of the 3D projection of the full attractor is essentially the same, thus suggesting that the subspace statistics of the phasespace 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 loworder 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 phasespace 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 phasespace 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 phasespace topology (Hentschel and Procaccia, 1983; Ott, 2002), while multiscale dynamical systems can be characterized by the statistics of the phasespace scaling changing as different realspace scales are considered (Alberti et al., 2020a). For this purpose, we investigate the statistics of the phasespace 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_{2} 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{\}}$.
As expected, the multiscale correlation dimension for each MIMF decreases with increasing timescale, being representative of a more regular, less stochastic/chaotic, behavior of largescale MIMFs as compared with the shortterm ones (Alberti et al., 2020a). Particularly, when approaching the largest timescales, ${D}_{\mathrm{2},\mathit{\tau}}\to \mathrm{1}$ suggesting the existence of fixedscale 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 scaleindependent behavior of D_{2,τ} is highlighted for both values of the friction coefficient C. This suggests that the shortterm 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 scaleindependent behavior is found for the atmosphere for both values of C, while a scaledependent behavior is observed for the ocean. The former can be easily related to the dominant role of the shortterm variability for the atmosphere, while the latter is a reflection of the longterm 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 scaleindependent behavior of D_{2,τ}.
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 largescale 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 scaleindependent 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 higherdimensional chaotic system for reduced ocean–atmosphere coupling (i.e., C=0.008) as opposed to a lowerdimensional 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 scaleindependent D_{q,τ}, being representative of a highdimensional system whose prime dynamics occurs at short timescales and with little effects of largescale processes on the collective dynamics of the atmosphere. By contrast, a clearly scaledependent behavior is found for the oceanic variables, with the multiscale generalized dimensions decreasing at larger timescales, reflecting the effects of largescale dynamics dominating with respect to the shortterm 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.
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_{0}. 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 highfrequency “noise” due to the atmosphere. Deeper investigations will be devoted to clarify this point in future research.
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 secondorder multifractal phase transitions (Lovejoy and Schertzer, 2013). To mitigate these effects (e.g. Supplement) since we deal with the investigation of scaledependent 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 sixthorder and the fourthorder 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 shortterm variability and in terms of scalinglaw 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^{1} is the socalled 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.
We clearly see that for the atmosphere, there is a scaleindependent behavior of D_{q,τ} for all q, rendering the different curves almost invariant with respect to the scale. By contrast, a scaledependent 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 Mediumrange Weather Forecasts (ECMWF) ORA20C project (De Boisséson and Balmaseda, 2016; De Boisséson et al., 2017) that is a 10member ensemble of ocean reanalyses covering the complete 20th century using atmospheric forcing from the ERA20C reanalysis (https://www.ecmwf.int/en/forecasts/datasets/reanalysisdatasets/era20c). 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}}\mathrm{cos}\left(\mathit{\pi}y/{L}_{y}\right)$ and (ii) ${\mathit{\varphi}}_{\mathrm{2}}=\mathrm{2}\mathrm{sin}\left(\mathit{\pi}x/{L}_{x}\right)\mathrm{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 3D projection of the local atmosphere–ocean coupled dynamics for each region (see Vannitsem and Ekelmans, 2018, for more details).
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.
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 36dimensional dynamical system, with respect to the reanalysis data, being a 3dimensional 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 scaledependent 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 shortterm multifractal nature to a longterm 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/multiannual 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 scaledependent 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 phasespace 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 indepth investigation is outlined as a part of our future work.
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 phasespace scaling (Alberti et al., 2020a).
First, by means of the multivariate empirical mode decomposition (MEMD) we have been able to detect oscillating patterns with timedependent 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 3D 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 scaleindependent behavior of D_{q,τ} for all q, rendering the multifractal spectra almost invariant with respect to the timescale. By contrast, a scaledependent 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 ORA20C 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 36dimensional dynamical system) and the reanalysis data (3dimensional 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 scaledependent behaviors. Specifically, both the North Atlantic and the North Pacific regions are characterized by a scaledependent behavior, with decreasing D_{q,τ} as τ increases, with a transition from a shortterm multifractal nature to longterm 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/multiannual 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 scaledependent 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 phasespace 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 phasespace 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.
All codes used for the analysis and generating the figures can be obtained from the authors upon request.
The time series of the model used in the present article are available from the authors upon request. The reanalysis dataset is available on Zenodo: https://doi.org/10.5281/zenodo.1135134 (Vannitsem, 2018).
The supplement related to this article is available online at: https://doi.org/10.5194/esd128372021supplement.
TA, RVD, and SV designed the study. TA conducted the analysis and drafted the manuscript. RVD and SV contributed to the interpretation of the results. All authors contributed to the writing of the manuscript.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work has been partially supported by the German Federal Ministry for Education and Research (BMBF, grant no. 01LP2002B) and the Belgian Science Policy Office (Belspo) through the JPI Climate/JPI Oceans project ROADMAP.
This research has been supported by the Bundesministerium für Bildung und Forschung (grant no. 01LP2002B) and the Belgian Federal Science Policy Office (project ROADMAP).
This paper was edited by Rui A. P. Perdigão and reviewed by two anonymous referees.
Alberti, T.: Multivariate empirical mode decomposition analysis of Swarm data, Nuovo Cimento C, 41, 1–10, 2018. a
Alberti, T., Consolini, G., Ditlevsen, P. D., Donner, R. V., and Quattrociocchi, V.: Multiscale measures of phasespace trajectories, Chaos, 30, 123116, https://doi.org/10.1063/5.0008916, 2020a. a, b, c, d, e, f, g, h, i, j, k
Alberti, T., Giannattasio, F., De Michelis, P., and Consolini, G.: Linear Versus Nonlinear Methods for Detecting Magnetospheric and Ionospheric Current Systems Patterns, Earth Space Sci., 7, e2019EA000559, https://doi.org/10.1029/2019EA000559, 2020b. a
Alexander, M. A., Bladé, I., Newman, M., Lanzante, J. R., Lau, N.C., and Scott, J. D.: The atmospheric bridge: The influence of ENSO teleconnections on air–sea interaction over the global oceans. J. Climate, 15, 2205–2231, 2002. a
Ambaum, M. H. P., Hoskins, B. J., and Stephenson, D. B.: Arctic Oscillation or North Atlantic Oscillation?, J. Climate, 14, 3495–3507, 2001. a
Badii, R. and Politi, A.: Hausdorff dimension and uniformity factor of strange attractors, Phys. Rev. Lett., 52, 1661–1664, 1984. a
Chatfield, C.: The analysis of time series – an introduction, 6th edn., Chapman and Hall/CRC, London, 2016. a, b
Czaja, A. and Frankignoul, C.: Observed impact of Atlantic SST anomalies on the North Atlantic Oscillations, J. Climate, 15, 606–623, 2002. a, b
D'Andrea, F., Czaja, A., and Marshall, J.: Impact of anomalous ocean heat transport on the North Atlantic Oscillation, J. Climate, 18, 4955–4969, 2005. a
De Boisséson, E. and Balmaseda, M.: An ensemble of 20th century ocean reanalyses for providing ocean initial conditions for CERA20C coupled streams, ERA report series, 24, European Centre for Mediumrange Weather Forecasts, Reading, UK, 2016. a, b, c
De Boisséson, E., Balmaseda, M., and Mayer, M.: Ocean heat content variability in an ensemble of twentieth century ocean reanalyses, Clim. Dynam., 50, 3783–3798, https://doi.org/10.1007/s0038201738450, 2018. a, b
De Cruz, L., Demaeyer, J., and Vannitsem, S.: The Modular ArbitraryOrder OceanAtmosphere Model: MAOOAM v1.0, Geosci. Model Dev., 9, 2793–2808, https://doi.org/10.5194/gmd927932016, 2016. a, b
De Cruz, L., Schubert, S., Demaeyer, J., Lucarini, V., and Vannitsem, S.: Exploring the Lyapunov instability properties of highdimensional atmospheric and climate models, Nonlin. Processes Geophys., 25, 387–412, https://doi.org/10.5194/npg253872018, 2018. a
Del Rio Amador, L., and Lovejoy, S.: Predicting the global temperature with the Stochastic Seasonal to Interannual Prediction System (StocSIPS), Clim. Dynam., 53, 4373–4411, https://doi.org/10.1007/s00382019047914, 2019. a
Del Rio Amador, L. and Lovejoy, S.: Longrange forecasting as a past value problem: Untangling correlations and causality with scaling, Geophys. Res. Lett., 48, e2020GL092147, https://doi.org/10.1029/2020GL092147, 2021. a
Donner, R. V., Heitzig, J., Donges, J. F., Zou, Y., Marwan, N., and Kurths, J.: The geometry of chaotic dynamics – a complex network perspective, Eur. Phys. J. B, 84, 653–672, 2011. a, b
Faranda, D., Messori, G., and Vannitsem, S.: Attractor dimension of timeaveraged climate observables: insights from a loworder oceanatmosphere model, Tellus A, 71:1, https://doi.org/10.1080/16000870.2018.1554413, 2019. a
Farmer, J. D., Ott, E., and Yorke, J. A.: The dimension of chaotic attractors, Physica D, 7, 153–180, 1983. a
Farneti, R.: Modelling interdecadal climate variability and the role of the ocean, WIREs Clim. Change, 8, e441, https://doi.org/10.1002/wcc.441, 2017. a
Feliks, Y., Ghil, M., and Robertson, A. W.: The atmospheric circulation over the North Atlantic as induced by the SST field, J. Climate, 24, 522–542, https://doi.org/10.1175/2010JCLI3859.1, 2011. a
Flandrin, P., Rilling, G., and Goncalves, P.: Empirical Mode Decomposition as a Filter Bank, IEEE Sign. Proc. Lett., 11, 112, https://doi.org/10.1109/LSP.2003.821662, 2004. a
Franzke, C. L. E., Barbosa, S., Blender, R., Fredriksen, H.B., Laepple, T., Lambert, F., Nilsen, T., Rypdal, K., Rypdal, M., Scotto, M. G., Vannitsem, S., Watkins, N. W., Yang, L., and Yuan, N.: The Structure of Climate Variability Across Scales, Rev. Geophys., 58, e00657, https://doi.org/10.1029/2019RG000657, 2020. a, b
Gastineau, G., D'Andrea, F., and Frankignoul, C.: Atmospheric response to the North Atlantic Ocean variability on seasonal to decadal time scales, Clim. Dynam., 40, 2311–2330, https://doi.org/10.1007/s0038201213330, 2013. a
Grassberger, P.: Generalized dimensions of strange attractors, Phys. Lett. A, 97, 227–230, 1983. a, b, c
Hentschel, H. G. E. and Procaccia, I.: The infinite number of generalized dimensions of fractals and strange attractors, Physica D, 8, 435–444, 1983. a, b, c, d, e, f, g, h, i
Huang, N. E. and Wu, Z.: A review on HilbertHuang transform: Method and its applications to geophysical studies, Rev. Geophys., 46, RG2006, https://doi.org/10.1029/2007RG000228, 2008. a, b
Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng, Q., Yen, N.C., Tung, C. C., and Liu, H. H.: The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis, Proc. R. Soc. Lon. Ser.A, 454, 903, https://doi.org/10.1098/rspa.1998.0193, 1998. a, b, c, d, e, f, g, h, i, j, k
Jin, F.F.: Tropical OceanAtmosphere interaction, the Pacific cold tongue, and the ElNiñoSouthern Oscillation, Science, 274, 76–78, 1996. a
Kaplan, J. L. and Yorke, J. A.: Chaotic behavior of multidimensional difference equations, in: Functional Differential Equations and Approximation of Fixed Points, Lect. Notes Math., edited by: Peitgen, H.O. and Walther, H.O., 730, 204–227, Springer, Berlin, Heidelberg, https://doi.org/10.1007/BFb0064319, 1979. a
Kravtsov, S., Dewar, W. K., Berloff, P., Ghil, M. and McWilliams, J. C.: A highly nonlinear coupled mode of decadal variability in a midlatitude oceanatmosphere model, Dyn. Atmos. Oceans, 43, 123–150, 2007. a
Legras, B. and Ghil, M.: Persistent anomalies, blocking, and variations in atmospheric predictability, J. the Atmos. Sci., 42, 433–471, 1985. a
L'Hévéder, B., Codron, F., and Ghil, M.: Impact of anomalous northward oceanic heat transport on global climate in a slabocean setting, J. Climate, 28, 2650–2664, 2014. a
Liu, Z.: Dynamics of interdecadal climate variability: A historical perspective, J. Climate, 25, 1963–1995, 2012. a, b
Lorenz, E. N.: Formulation of a loworder model of a moist general circulation, J. Atmos. Sci., 41, 1933–1945, 1984. a
Lovejoy, S.: The halforder energy balance equation – Part 1: The homogeneous HEBE and long memories, Earth Syst. Dynam., 12, 469–487, https://doi.org/10.5194/esd124692021, 2021. a
Lovejoy, S. and Schertzer, D.: The Weather and Climate: Emergent Laws and Multifractal Cascades, Cambridge University Press, Cambridge, 496 pp., 2013. a, b, c, d, e
Lovejoy, S., Procyk, R., Hébert, R., and Del Rio Amador, L.: The fractional energy balance equation, Q. J. Roy. Meteorol. Soc., 147, 1964–1988, https://doi.org/10.1002/qj.4005, 2021. a
Mandelbrot, B.: How Long Is the Coast of Britain? Statistical SelfSimilarity and Fractional Dimension, Science, 156, 636–638, https://doi.org/10.1126/science.156.3775.636, 1967. a
Meehl, G. A., Arblaster, J. M., and Loschnigg, J.: Coupled Ocean–Atmosphere Dynamical Processes in the Tropical Indian and Pacific Oceans and the TBO, J. Climate, 16, 2138–2158, 2003. a
Mosedale, T. J., Stephenson, D. B., Collins, M., and Mills, T. C.: Granger Causality of Coupled Climate Processes: Ocean Feedback on the North Atlantic Oscillation, J. Climate, 19, 1182–1194, 2006. a
Neelin, J. D., Latif, M., and Jin, F.F.: Dynamics of Coupled OceanAtmosphere Models: The Tropical Problem, Annu. Rev. Fluid Mech., 26:1, 617–659, 1994. a
Nese, J. M. and Dutton, J. A.: Quantifying predictability variations in a loworder oceanatmosphere model: A dynamical system approach. J. Climate, 6, 185–203, 1993. a
Ott, E.: Chaos in Dynamical Systems – 2nd Edn., Cambridge University Press, Cambridge, 490 pp., https://doi.org/10.2277/0521811961, 2002. a, b, c
Pawelzik, K. and Schuster, H. G.: Generalized dimensions and entropies from a measured time series, Phys. Rev. A, 35, 481–484, 1987. a
Philander, S. G. H.: El Niño and the Southern Oscillation, Academic Press, New York, 1990. a
Poli, P., Hersbach, H., Tan, D. G. H., Dee, D. P., Thépaut, J.N., Simmons, A., Peubey, C., Laloyaux, P., Komori, T., Berrisford, P., Dragani, R., Trémolet, Y., Hólm, E. V., Bonavita, M., Isaksen, L., and Fisher, M.: The data assimilation system and initial performance evaluation of the ECMWF pilot reanalysis of the 20thcentury assimilating surface observations only, ERA report series, 14, European Centre for Mediumrange Weather Forecasts, Reading, UK, 2015. a, b
Primavera, L. and Florio, E.: Parallel Algorithms for Multifractal Analysis of River Networks, in: Numerical Computations: Theory and Algorithms, Lect. Notes Comput. Sc., edited by: Sergeyev Y. and Kvasov D. 11973, Springer, Cham, 2020. a
Rehman, N. and Mandic, D. P.: Multivariate empirical mode decomposition, P. Roy. Soc. AMath. Phy., 466, 1291–1302, 2010. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o
Roebber, P. J.: Climate variability in a loworder coupled atmosphereocean model, Tellus A, 47, 473–494, 1995. a
Steinhaus, H.: Length, shape and area, Colloq. Math., 3, 1–13, 1954. a
Timmermann, A., Jin, F.F., and Abshagen, J.: A nonlinear theory for El Niño bursting, J. Atmos. Sci., 60, 152–165, 2003. a
Van der Avoird, E., Dijkstra, H. A., Nauw, J. J., and Schuurmans, C. J. E.: Nonlinearly induced lowfrequency variability in a midlatitude coupled oceanatmosphere model of intermediate complexity, Clim. Dynam., 19, 303–320, 2002. a
Vannitsem, S.: The role of the ocean mixed layer on the development of the North Atlantic Oscillation: A dynamical system's perspective, Geophys. Res. Lett., 42, 8615, https://doi.org/10.1002/2015GL065974, 2015. a, b, c
Vannitsem, S.: Predictability of largescale atmospheric motions: Lyapunov exponents and error dynamics, Chaos, 27, 032101, https://doi.org/10.1063/1.4979042, 2017. a, b, c, d
Vannitsem, S.: Time series used in the manuscript “Causal dependences between the coupled oceanatmosphere dynamics over the Tropical Pacific, the North Pacific and the North Atlantic” [Data set], Zenodo, https://doi.org/10.5281/zenodo.1135134, 2018. a
Vannitsem, S. and Ekelmans, P.: Causal dependences between the coupled ocean–atmosphere dynamics over the tropical Pacific, the North Pacific and the North Atlantic, Earth Syst. Dynam., 9, 1063–1083, https://doi.org/10.5194/esd910632018, 2018. a, b, c, d
Vannitsem, S. and Ghil, M.: Evidence of coupling in ocean atmosphere dynamics over the North Atlantic, Geophys. Res. Lett., 44, 2016–2026, https://doi.org/10.1002/2016GL072229, 2017. a
Vannitsem, S., Demaeyer, J., de Cruz, L., and Ghil, M.: Lowfrequency variability and heat transport in a loworder nonlinear coupled oceanatmosphere model, Physica D, 309, 71–85, 2015. a, b, c, d, e, f, g, h
Vannitsem, S., Demaeyer, J., and Ghil, M.: Extratropical lowfrequency variability with ENSO forcing: A reducedorder coupled model study. J. Adv. Model. Earth Sy., 13, e2021MS002530, https://doi.org/10.1029/2021MS002530, 2021. a
Van Veen, L.: Overturning and wind driven circulation in a loworder ocean–atmosphere model, Dynam. Atmos. Oceans, 37, 197–221, 2003. a
Wang, C.: Threeocean interactions and climate variability: a review and perspective, Clim. Dynam. 53, 5119–5136, 2019. a
Wunsch, C. and Ferrari, R.: Vertical mixing, energy, and the general circulation of the ocean, Annu. Rev. Fluid Mech., 36, 281–314, 2004. a
Xue, P., Malanotte Rizzoli, P., Wei, J., and Eltahir, E. A. B.: Coupled ocean atmosphere modeling over the Maritime Continent: A review, J. Geophys. Res.Oceans, 125, e2019JC014978, https://doi.org/10.1029/2019JC014978, 2020. a, b
Another direct measure is the socalled codimension of the mean $c=d{D}_{\mathrm{0}}$ where d is the dimension of the phase space (e.g., Lovejoy and Schertzer, 2013). For the sake of simplicity we prefer to use here only the multifractal width since it can be easily derived from D_{q,τ} and not to introduce an additional alternative concept.