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 5 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/multivariate generalized fractal dimensions. With these two approaches, we show that the oceanatmosphere 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 10 different, while for strong coupling for which coupled modes develop, the scaling properties are more similar especially at longer time scales. 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 topological and geometric features for different regions, being related to the different drivers and processes occurring at different timescales in the coupled atmosphere-ocean system. 15 Our approach can therefore be used to diagnose the strength of coupling in real applications.


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 extra-tropical Northern 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 modelling (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 25 et al., 2020, and reference therein), highlighting how the atmospheric low-frequency variability (LFV) is strictly 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 systems point of view by developing suitable reduced order ocean-atmosphere models 30 dealing with the modelling 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.
The current work presents an investigation on how a recently introduced concept of multiscale generalized fractal dimensions can be used to analyze the topological and geometric properties of attractors in coupled ocean-atmosphere systems (Alberti et 35 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 (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 40 on partitioning the phase-space into hypercubes of size to define a suitable invariant measure through the filling probability of the i−th 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) 45 the information dimension and the correlation dimension with Θ(· · · ) being the Heaviside function. All these fractal dimension measures, as well as their higher order extensions D q , 50 have been used to characterize the global dynamical, topological, and geometric properties of a given system (Hentschel and Procaccia, 1983), however, without exploring how these properties evolve at different scales (Alberti et al., 2020a). More 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 time-scales.
Here we used the original version of the model  where the four relevant fields, i.e., the barotropic and baroclinic atmospheric streamfunctions, the ocean streamfunction and the ocean temperature, are given by 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 95 larger values of C.

Methods
Traditional multivariate and/or spatiotemporal data analysis methods are commonly based on fixing an orthogonal decomposition basis, satisfying certain mathematical properties of completeness, convergence, linearity, and stationarity (Chatfield, 2016). However, these conditions are not usually met when real-world geophysical data are analyzed, which calls for more 100 adaptive methods (Huang et al., 1998). 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 (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 105 fractal dimensions (Hentschel and Procaccia, 1983).

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 110 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 needs 115 to be carefully implemented for multivariate processes (Rehman and Mandic, 2010), since it is based on the cubic spline 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 {u(t)}| t∈T and the lower {l(t)}| t∈T envelopes; The previous steps are iteratively repeated until the obtained candidate MIMF {h(t)}| t∈T (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)}| t∈T can be filtered out from the data. Hence, we can write 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)}| t∈T which form the multivariate decomposition basis, plus a multivariate residue {r(t)}| t∈T .
For each MIMF we can define a k −variate mean timescale as representing the typical oscillation scale of the j−th mode for the k -th univariate component c j,k extracted from the mul- Similarly, by ensemble averaging over the k-dimensional space we can introduce the concept of a multivariate mean timescale as with · · · k denoting an ensemble average over the k-dimensional space. Thus, the k −variate timescale τ j,k is evaluated for each mode and for each k −dimensional data, while the multivariate mean timescale τ j is the mean over all k ∈ [1, k].
Moreover, as for univariate EMD (Huang et al., 1998), we can introduce the concepts of instantaneous amplitudes {a j (t)}| t∈T and phases {φ j (t)}| t∈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)}| t∈T = {a j (t)}| 2 t∈T . Thereby, we can characterize the 150 spectral content by introducing an alternative yet equivalent definition of the power spectral density (PSD) as with σ 2 (τ ) 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)}| t∈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 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).

Multivariate and multiscale generalized fractal dimensions
The behavior of complex systems usually consists of a collection of scales whose dynamical features determine their collective 160 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 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 i−th hypercube B s * ,τ ( ) of size centered at the point {s * } on the considered (d−dimensional) phase-space of {s 1 (t)}| t∈T can be defined as 170 By defining a q−th 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 phase-space properties at different 175 scales by deriving the generalized dimensions (Alberti et al., 2020a). Summarizing: 1. we extract multiscale components from {s(t)}| t∈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) with j = 1, . . . , N j (by construction, MIMFs are ordered from short to long scales, i.e., τ j < τ j if j < j ); 4. we evaluate the generalized dimensions D q,τ from F j (t) for each j (i.e., for each scale τ j ),; 5. we evaluate the singularities and singularity spectrum f 185 From Eq. (13) we can inspect 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. Finally, it is expected (for ensuring convergence) that when j → N j then D q,τ → D q , with D q being the standard generalized fractal dimensions proposed by Hentschel and Procaccia (1983). Figure 1 reports the 3-D projection of the full system attractor in the subspace (T o,2 , Ψ o,2 , ψ a,1 ) for two representative values of the friction coefficient C (0.008 and 0.015 kg m −2 s −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 195 transport within the ocean (Ψ o,2 ), and by the vertically averaged zonal flow within the atmosphere (ψ a,1 ), respectively.

Multivariate Empirical Mode Decomposition
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, 2015Vannitsem, , 2017.

200
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 205 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 with a purely stochastic dynamics. Indeed, in both 210 cases we have used N = 10 5 data points, thus the expected number of MIMFs is N noise j = 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 is 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 dif-215 ferent 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 220 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 Fig. 3. It can be clearly noted that the oceanic 225 variability mainly contributes to the low-frequency dynamics (E χ,τ > 95% for χ = {Ψ o,i , T o,i } and τ 10 4 days), while the atmosphere is mainly characterized by short-term variability for C = 0.008 (E χ,τ > 95% for χ = {ψ a,i , θ a,i } and τ 10 days) 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.

230
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 235 of the full system attractor onto the subspace (T o,2 , Ψ o,2 , ψ a,1 ) and compare it with the projection obtained by considering all timescales (as in Figure 1). Thus, for each variable χ = {ψ a,i , θ a,i , Ψ o,i , T o,i } we can define a reconstruction based on empirical modes, R χ,95% , as with {c χ,j (t)}| t∈T being the j −th multivariate empirical mode extracted via the MEMD of the variables χ. The 3-D pro-240 jection onto the subspace (T o,2 , Ψ o,2 , ψ a,1 ) of R χ,95% is shown in Fig. 4, while Tab. 1 summarizes the mode indices j and corresponding k −variate timescales τ j ,k (see Eq. (5)) used for the reconstruction.

Multiscale generalized fractal dimensions
Under general conditions, the complexity of a dynamical system can be conveniently investigated by means of the nonlinear 250 properties of its phase-space trajectory (e.g., its attractor or repellor 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 255 Procaccia, 1983;Ott, 2002), while multiscale dynamical systems can be characterized by topological properties changing as different scales are considered (Alberti et al., 2020a). For this purpose, we investigate the topological properties of the attractor of the coupled ocean-atmosphere model by evaluating the multiscale generalized fractal dimensions described in Section 3.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 modesviolet symbol). Each panel also shows the 95% confidence intervals as error bars.
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 2,τ → 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 265 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. However, the role of C clearly emerges in determining the values of D 2,τ , being lower for the larger C value. Indeed, D 2,τ ∼ 8 for C = 0.008, while D 2,τ ∼ 1.5 for C = 0.015. This reflects the different topological properties of the attractor of the full system associated with a different dynamical behavior of the model variables (Faranda et al., 2019). However, the most interesting features emerge 270 when the different variables of both atmosphere and ocean are separately investigated by means of the multiscale generalized fractal dimensions.
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 , D q,τ > D q ,τ ∀τ (Alberti et al., 2020a;Hentschel and Procaccia, 1983). Moreover, when evaluating the multiscale generalized fractal dimensions for each MIMF separately (e.g., Figs. 7(a) and 8(a)) a decreasing value for D j q is found as τ increases, with all D j q converging towards the same value of 1 at large timescales. Conversely, when the D q,τ are evaluated for reconstructions based on MIMFs a completely different behavior emerges between the oceanic and the atmospheric variables. In this case, the atmospheric variables are characterized by scale-independent D q,τ , being representative of a high-dimensional system whose 280 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 shortterm one for the ocean variability.
By estimating the Lyapunov spectra (not shown) separately for the ocean and the atmosphere we obtained that for C = 0.008 285 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 a bit larger than 4. Following the Kaplan-Yorke conjecture (Kaplan and Yorke, 1979), the Lyapunov dimension can be used as a proxy of the Hausdorff and, hence, capacity dimension. 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 0,τ ≈ 4. This quantitative disagreement could be related to the fact that 290 the ocean can be viewed as a relatively stable system perturbed by high-frequency "noise" provided the atmosphere. Deeper investigations will be devoted to clarify this point in future research.    Figure 8. Same as in Fig. 7, but for C = 0.015.
As a further step, we evaluate the full spectrum of generalized fractal dimensions for each MIMF by considering all moments q ∈ [−20, 20], thus providing an estimate of the asymptotic values D ±∞,τ . This analysis allows characterizing how the 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., Franzke et al., 2020). The D q,τ spectrum for q ∈ [−20, 20] 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, the full system exhibits signatures of multifractality at all timescales, especially at very short and very long timescales. By defining the multifractal width as ∆ ∞ . = D −∞,τ −D +∞,τ we observe (see 300 Fig. 10(a,b), black circles) that ∆ ∞ = 3 for τ ∈ [τ S , τ L ] days, while ∆ ∞ > 3 for both τ < τ S and τ > τ L , with τ S ∼ 20 days 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 by considering all orders q ∈ [−20, 20] for each subsystem (i.e., atmosphere and ocean) individually. For both values of C, the corresponding results are shown in Fig. 11. 305 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,τ1 ∈ [5, 8] to D q,τ17 ∈ [1, 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 310 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. 11(d)). 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 −∞,τ and D +∞,τ we can clearly see that larger values, of the order of 4, are found for the atmosphere, at almost all timescales (and especially at shorter timescales), for both values of C. Conversely, larger values are 315 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 2, 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).

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 320 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 325 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 (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 1 = √ 2 cos (πy/L y ), and (ii) φ 2 = 2 sin (πx/L x ) sin (2πy/L y ) (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), and a third 335 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 up a 3-D 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 340 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 345 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 gen- ming up from j = 1 to N j (D j q ). Figure 13 shows the corresponding results for the North Atlantic region, the North Pacific region, and the Tropical Pacific region, 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 355 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 the one hand, both the North Atlantic and the North Pacific regions (see Fig. 13(d,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 360 scale ( Fig. 13(a,b)) we find evidence for a decreasing ∆ ∞ as τ increases, being representative of a transition from a shortterm multifractal nature to 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. 13(c)), at annual/multi-annual timescales (i.e., τ ∼ 1 − 8 years), being also characterized by 365 the largest values of the multiscale generalized fractal dimensions (see Fig. 13(f)). 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 oceanatmosphere model and the reanalysis data show some qualitatively similar behavior of the multiscale generalized fractal di-370 mensions, although 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 topological and geometric features 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 375 dependency in more detail. Such an in-depth investigation is outlined as a part of our future work.

Conclusions
We have provided a first time systematic investigation of the multiscale dynamics of a reduced order coupled ocean-atmosphere model  as described by means of the topological and geometric features (Alberti et al., 2020a). Second, by exploiting the novel concept of multiscale/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 390 scaling properties are more similar especially at longer time scales. 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 395 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 400 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 the scale τ . Interestingly, the multiscale 405 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., 410 τ ∼ 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 other extratropical regions.
Our findings for both the model and the reanalysis data suggest that our approach can be used to diagnose the strength of cou-415 pling in the ocean-atmosphere system and to investigate the topological features of the system. We have demonstrated that the