In fluid flow, turbulent mixing is characterized by chaotic interaction. It is always accompanied by an energy cascade, with eddies gradually moving across scales until the energy is exhausted by viscosity. Eddy interaction is a key factor in the redistribution of energy, momentum and carbon, thereby determining the oceanic general circulation [Ferrari and Wunsch (2009); Busecke and Abernathey (2019)]. However, the mechanism of the energy cascade, such as the critical scaling features, is still poorly understood [Yeung et al. (2015); Buaria et al. (2020)].

According to the celebrated poem derived from Richardson atmosphere analysis [Richardson and Lynch (1922)], turbulence energy is injected from eddies at large scales and transferred to finer scales. The fundamental concept of energy cascade is first quantified by the seminal K41 theory [Kolmogorov (1991)], in which the injected energy is redistributed among eddies of different scales. Kolmogorov (1962) first proposed the first multiplicative cascade model, and it describes a successive uneven split process in energy cascade. However, the energy cascade mechanism was not rigorously demonstrated. Furthermore, a significant deviation from the scaling law was discovered for high orders. Subsequent experiments showed that the energy cascade is intermittent rather than homogeneous, with strong bursts of energy occurring during relatively stable periods [Batchelor and Townsend (1949)]. In essence, the energy cascade is intermittent, with strong bursts of activity occurring at irregular time intervals and in localized patches of space [McMillan et al. (2016)]. As a result, researchers have begun to consider the effects of intermittency with great enthusiasm, proposing the fractal concept [Biferale (2003)] to characterize the heterogeneity of energy transfer. Here, intermittency is explained as the concentration of turbulent energy in ‘active eddies’ while cascading towards finer scales. In the β model [Frisch et al. (1978), Frisch et al. (2006)], self-similar structures of generated sizes fill only a fixed part of the space of a given scale in the transfer. However, the monotonous scaling contradicts the results of the experiments. A continuous infinity of scaling exponents is expected to exhibit in the energy cascade. Furthermore, it appears that the intermittency resulting from the interactions of multi-scale eddies cannot be properly characterized by similarly uniform evolution features. Subsequently, the simple β model is generalized to random fractions [Benzi et al. (1984)], and it enriches the geometric dimension of the cascade evolution and introduces the multifractal spectrum [Sreenivasan (1991)]. According to this, the multifractal model translates the energy cascade into a spontaneous process in which large scale energetic vortices are supposed to simultaneously generate a multitude of eddies on all spatial scales [Carbone et al. (1995), Zhao (2003); Dupont et al. (2020)], accompanied by energy transfer until exhaustion. The superiority of fractal models is clearly due to the combination of self-similarity and intermittency. It is precisely this that has contributed to the understanding of multi-scale behavior. However, the inhomogeneity of scale interaction is partly ignored, which hides the full recognition of the energy cascade.

Common methods for investigating energy cascade include structure function, multifractal, scaling exponent, etc. Conventional Navier-Stokes equations are used to analyze energy cascade and intermittency in Fourier space [Dascaliuc and Grujic (2011); Sahoo and Biferale (2018)]. A morphing continuum theory is used to evaluate the contribution of small-scale structures to the turbulence energy transfer [Cheikh et al. (2019). The partial differential equation is employed in the prototype model to investigate the irreversibility of the energy cascade, and the nonstationary interaction is introduced in it with the application of the nonlinear term [Josserand et al. (2017)]. The generalized Holder means are used to study the relationship between the energy cascade and the strong burst events in the dissipative range [Vela-Martin (2022)]. The multifractal method is utilized to assess the dependence of energy cascade with dissipation in the Northwest Atlantic Ocean [Isern-Fontanet and Turiel (2021)]. However, a set of conditions corresponding to the detailed balance in turbulent systems must be satisfied a prior [Lee et al. (2018); McKeown et al. (2020)]. Subtle adjustments of statistical terms are required for the physical equations, which can only be checked by a posteriori [Cerbus and Chakraborty (2017); Xie and Buhler (2018)]. Meanwhile, they focus only on the source or sink terms that result in energy transfer at a single point [Cardesa et al. (2015)], which makes the dynamic perception in practical ocean systems impossible. Furthermore, estimating the energy transfer rate based on numerical simulations is extremely challenging. The number of degrees of freedom required for flow configuration increases sharply with increasing Reynolds number. In addition, the sampled grids of numerical simulations prevent measurements over a relatively large region [Klein et al. (2019)].

As a result, there is an urgent need to investigate and interpret the multi-scale eddy interaction using appropriate methods to achieve a better description of the energy cascade dynamics based on the amount of turbulent data collected in experimental observations [Evans et al. (2022)]. The complex network is considered to be an effective method to study the chaotic system. In general, a complex network is superior to abstract and simplify the stochastic process, by mapping the underlying dynamics into topological elements [Shirazi et al. (2009)]. Some studies have reported the study of turbulent and vortex flows on topological networks [Sorriso-Valvo et al. (2007); Iacobello et al. (2021a); Iacobello et al. (2021b)]. It has been possible to relate topological features to the corresponding physical behavior and to explore the spatial details of turbulent flow [Charakopoulos et al. (2014); Iacobello et al. (2018)]. A single global network [Scarsoglio et al. (2016)] was constructed to investigate the spatial patterns in forced isotropic turbulent fluid. Self-similarity of energy dissipation rate series was identified using the visibility algorithm in three-dimensional fully developed turbulence [Liu et al. (2010)]. A small-world network was proposed to characterize energy transfer in the turbulent cascade [Gurcan (2021)]. However, the evolution of the energy cascade was driven by the additional degrees of freedom provided by this small-world network. Furthermore, turbulence always encounters eddies of various scales [Cardesa et al. (2015)], and it is extremely difficult to unravel the scale effects unambiguously using only the single layer network. Nowadays, to better understand the mechanism of the chaotic system, people are starting to study the ascending effect of multi-layer networks composed of relevant single-layer networks. Gao [Gao et al. (2018); Gao et al. (2020)] developed a multi-layer network based on multiple entropy to characterize the nonlinear behavior and evolution of gas-liquid flow. These results paved the way for the realization of multivariate fusion as well as the identification, categorization, and exploration of dynamic turbulence features based on topological properties in network architectures.

In this work, an energy cascade multi-layer network (ECMN) model is proposed to reveal the underlying intermittency and inhomogeneous eddy interactions in the energy cascade. We focus on the classical forward energy cascade and all of the turbulence energy is ideally assumed to be transferred from large scale to small scale structures. The microstructure data were measured in the South Yellow Sea (SYS) and the South China Sea (SCS). Physical dynamics are mapped onto topological structures through effective node definition, valid edge determination and structure aggregation. Properties of the ECMN and weighted single-layer network (WSLN) allow the exploration of multi-scale interaction in the energy cascade. Researches have supported the relationship between dissipation and energy cascade [Cardesa et al. (2015); Cardesa et al. (2017); Ballouz et al. (2020); Vela-Martin (2022)]. Therefore, a network structure coefficient κ based on topological framework is proposed to estimate the energy transfer strength in turbulent mixing and to validate the effectiveness of ECMN. The results show that the use of a multi-layer network can effectively characterize the evolution of energy cascade and quantify the energy transfer, providing a novel insight for future research.

The chaotic system is transformed into a topological structure using complex network theory. The ECMN is constructed using nodes definition and edges determination. The dynamical variables of the system are encoded in nodes, while the edges represent energy transfer across different scales. The ECMN’s intricate structure is then condensed into the WSLN, which is used to examine the energy cascade’s scale properties. Meanwhile, features of ECMN are utilized to parameterize energy transfer strength, providing a topological approach to scrutinize turbulent mixing. The skeleton is presented in Figure 1.

Because of their crucial role in the ocean circulation, data measured from the SYS [Pang et al. (2017); Teng et al. (2017); Yuhua et al. (2017)] and the SCS [Chunhua et al. (2019); Tan et al. (2021)], which are characterized by energetic currents (Figure 2A), were analyzed. The measurements were taken on 22-24 September 2021 at St. A (35°32'N,121°14'E, the mean water depth around 38 m) and on 12 December 2020 at St. B (19°39'N,115°27'E, with the mean water depth close to 1500 m). Due to the limitation of cable length, only the fluid statistics over 300 m were observed at St. B.

Vertical microstructure and mixing signatures were examined using hydrographic data. The vertical microstructure was investigated using the MicroRider instrument (Figure 2B), manufactured by the Rockland Scientific International (RSI). The MicroRider nose cone contains a pair of orthogonal velocity shear probes, a high-resolution temperature probe, an accelerometer and a pressure sensor, with sampling rates all at 512 Hz. During the observations, the MicroRider instrument was released from the stern deck, which was loosely tethered with a tether (Figure 2C).

The dissipation rate of turbulent kinetic energy ε [Wolk et al. (2002), Gregg (1999), Roget et al. (2006)], which quantifies the intensity of the turbulent mixing, can be obtained under the assumption of isotropic turbulence by integrating the power spectra over the wavenumber k.

where ν is the kinematic molecular viscosity (nearly 1 × 10 − 6   m 2 s − 1 ( ∂ u / ∂ z ) 2 ¯ is the variance of vertical shear, and Φ(k) is the power spectrum of shear in Equation 5. The power spectrum is integrated within the inertial domain ( k m i n <   k   <   k m a x ) . Here, k_min and k_max are the lower and upper wavenumber limits for integration, respectively. In the operation, the k_min is assumed to be 1 cpm (cycles per meter). The low wavenumber range ( k   <   k m i n ) cannot be detected by the MicroRider. The apparent high wavenumber peaks ( k   >   k m a x ) are caused by instrument vibration which interferes with the shear signals.

Figure 3 shows the vertical variation. There are peaks in the pitch variation corresponding to different types of disturbance. Above a depth of 5m, wave action caused instrument disturbance and slight pitching of the MicroRider (Figure 3A). The instrument entered a stable dive mode with an average rate of descent of 0.6 m/s (Figure 3B). A sharp drop in temperature indicates a strong pycnocline at 25m at St A and 105m at St B (Figure 3C). The pycnocline separates the water column into well-mixed surface, thermocline and bottom boundary layers in the horizontal direction. The stratification inhibits the turbulent mixing and vertical redistribution of marine material [Sharples and Simpson (2012)], as indicated by the sharp decrease in dissipation rate in Figure 3F. Below the pycnocline, a seasonal cold water mass occupies the water column and dominates the hydrological flow evolution [Fangli et al. (2011); Kyung-Hee et al. (2012); Fan (2016), Brown et al. (1999)], which is growing in spring, maturing in summer, decreasing in autumn and winter. Meanwhile, strong stratification causes relative tilting and disturbances in the vicinity of junctions. Therefore, the Goodman coherent noise reduction algorithm is applied to reduce the contamination in the shear measurements and to calibrate for the instrument vibration (Goodman et al. (2006)). The perturbation variation of shear is first extracted by comparing the original shear and its mean variation. The contamination was then removed by subtracting all of the coherent fluctuations between the shear perturbation variation (Figure 3E) and the accelerometer variables (Figure 3D). Figure 3G shows the shear power spectrum at stable depth, the blue and orange curves represent the shear power spectrum, the black line indicates the Nasmyth spectrum and the solid triangle represents the cutoff wavenumber. Within the cutoff wavenumber limit, the shear curves are found to almost overlap with the Nasmyth curves (Figure 3G), and this means that the shear spectrum is well matched to the standard Nasmyth spectrum, indicating effective data for the following turbulence analysis. Figure 3F illustrates the vertical turbulent mixing. Under the combination of wind and terrain frictions, the turbulent mixing intensifies in the surface and boundary layers of the SYS water column with a dissipation rate magnitude close to 10⁻⁷. Meanwhile, turbulent mixing is inhibited in the well-mixed thermocline layer with an average dissipation magnitude of 10⁻⁹. The vertical variation of the dissipation rate in the SYS is consistent with the description in previous researches [Xu et al. (2020); Song et al. (2021)]. Owing to the wind stress and buoyancy flux, the depth of mixing layer depth in the SCS is much deeper in winter. Due to the limited cable length, only the microstructure data in the surface and stratification layers were measured. The thermocline also impedes turbulent mixing in the water column in the SCS. Meanwhile, the wind stress stimulates turbulent mixing in the surface and its corresponding mixing magnitude is nearly 10⁻⁷, which is two orders higher than that in the thermocline (10⁻⁹) [Iossif (2012); Wang et al. (2012); Xiao et al. (2013)].

The entire shear signal is divided into uniform segments with 1 m to characterize the turbulent evolution in each vertical profile. Considering the diving stage of the instrument, depth segments of 6 m, 16 m and 26 m derived from the profile collected at St. A are selected to show their corresponding multi-layer frameworks. For the same reason, network structures are inferred from depth segments of 30 m, 80 m, 130 m, 180 m, 230 m and 280 m of the profile measured at St. B. These depth segments represent, in order, the initial phase, the sustained phase of stable diving, the turbulent phase crossing the thermocline and the near-bottom phase, respectively.

The energy cascade process is mapped by the network structure and topological features of the ECMN. Nodes and edges in the ECMN are depicted in a three-dimensional space with the following axes: x, which represents the time series; y, which indicates the intrinsic scale properties, network layers are arranged along the y axis, and the time scale decreases along the y−direction; and z, which indicates the local intermittency, indicating the energy accumulation of each node. Furthermore, a pair of simultaneous nodes (i, α) in layer α and (j, β) in layer β are adjacent to each other by an inter-layer edge e_αβ if they satisfy the qualification of edge determination. The inter-layer edge e_αβ represents the occurrence of energy transport between scales S_α and S_β .

Figures 4A–C show the corresponding ECMN built from shear segments at 6 m, 16 m, and 26 m, respectively. Shear signal in each segment is decomposed into several IMFs, and physical feature statistics are extracted at each sampling point of the IMF. Each node is a data collection containing time variable, scale variable, local intermittency variable and phase variable. Scattered nodes are assigned to a particular layer according to the scale feature, i.e. nodes derived from the same IMF are located at the identical layer. Here, layers are represented by transparent gray planes and nodes located on the identical layer are colored the same. Nodes with the shortest scale are located in the front right plane and are colored blue (Figure 4). Other nodes are located in the back layers and to the left as scale increases. Nodes with the largest scale feature are colored red in Figures 4A, C. Meanwhile, nodes with the largest scale are colored gray in Figure 4B. Furthermore, the location of nodes in a particular layer varies according to the time (x−direction) and energy properties (z−direction). In addition, the black dashed lines represent the connection between pairs of simultaneous nodes in different layers and they indicate the occurrence of energy transfer across scales.

Figures 5A–F show the corresponding ECMN constructed based on shear segments at 30 m, 80 m, 130 m, 180 m, 230 m and 280 m, respectively. Taking the Figure 5A as an example, the underlying energy cascade process in the 30 m depth segment is converted into a topological framework. The shear signal s(t) in the 30 m depth segment is decomposed into six IMFs. Wavelet analysis, LIM and HT algorithms are performed on IMF1-IMF6 to evaluate the underlying characteristics. For each IMF, we will obtain three variations of time, local intermittency and phase, and a fixed value of scale. Physical statistics obtained at the same sampling point and the scale value together form topological nodes. A total of 412 nodes are obtained in this segment, and they are spread out in different locations according to their feature variables. Meanwhile, the feature variables in each pair of nodes i and j from different layers are evaluated to estimate whether they satisfy the requirements and determine the existence of edge.

We assume that different layers have different physical meanings, which are associated with different powers and frequencies. The marginal Hilbert spectra (Figure 6) for each depth segment indicate the contribution and distribution of the turbulent energy across frequencies. According to the marginal Hilbert spectra, larger scale layers generally have higher peaks (such as the red curves in Figures 6A, C, G, I, and the grey curves in Figures 6B, D-F, H), representing the mean flow with overwhelmingly large energy and sustaining the primary ocean circulation. In addition, smaller scale layers are usually found with relatively low energy. Herein, we ideally assume that all of the turbulence energy is transferred from the large scale to the small scale. Thus, the larger scale layer is the energy-producing layer and the smaller scale layer corresponds to the energy-receiving scale.

The details of the ECMN are shown more clearly in Figures 7–10 (others in the Supplementary Material), where the overlapping edges between layers and their connected layers in Figures 4, 5 are separated into these sparse subgraphs, respectively. The sparse distribution of inter-layer edges indicates the inhomogeneity of energy cascade evolution. The energy transfer interaction between Layer 2 and its energy-receiving layers is relatively intense at the 6 m segment (Figure 7B). On the other hand, few inter-layer edges are found at the 26 m segment (Figure 9B). Meanwhile, the energy in Layer 1 (Figure 10A) is very stable and less energy is detected to be transferred from Layer 1. However, Layer 1 in Figure 5B is isolated from other layers, which maintain relatively quiescent energy, and no energy is transferred to small scales. In addition, the inter-layer edge distribution indicates the inhomogeneity in scale interaction. The inter-layer edges adjacent to Layer 3 and the its energy-receiving layers are shown in Figure 8C. It describes an energy transfer process, in which a small amount of energy in Layer 3 cascades to Layer 4, Layer 5 and Layer 6. In contrast, the scales derived from Layer 3 differ from those derived from Layer 2 (Figure 8B) or Layer 4 (Figure 8D). Energy in Layer 2 (Figure 8B) in found to be transferred to Layer 3-6 and turbulence energy in Layer 4 (Figure 8D) is found to be transferred to Layer 5-6. The inter-edges intensity increases as the energy of the scales decreases. Meanwhile, dramatically dense existences of inter-layer edges are found at the ends of subgraphs, indicating strong energy transfer between these miniature scales. The distribution of inter-layer edges and specific energy-receiving scales in each ECMN indicates that the characteristics of evolution and scale are dramatically intermittent and inhomogeneous.

Subsequently, the aggregation method is applied to investigate the scale property of the energy cascade under the multi-layer framework. Each depth segment is divided into two equal sections. The first half and the second half of the ECMN are aggregated into corresponding WSLN, respectively (Figures 11, 12). Here, the node sizes vary with the corresponding energy and the line width represents the weight of edge in the aggregated WSLN. Fine edges can be seen extending from the larger nodes (Figure 12A1). In particular, Layer1 is separated from the others in Figure 12A2. This implies that these energetic structures have settled into a relatively stationary state, with little energy observed from them. Meanwhile, nodes of intermediate energy are identified as relatively active and some of the energy contained within them is transferred to finer structures via inter-layer edges (Figure 12E1). The largest weight, indicating enormous inter-layer edges, appears in the pair of smallest nodes, e.g., Layer5 and Layer6 in Figure 11B1. And it is much more visible in Figure 11B2. Moreover, almost all of the nodes are connected to these small or miniature structures. In each step, the staple concentration of energy in the ‘active eddies’ tends to cascade directly to small or miniature structures [Josserand et al. (2017)] compared to other generated structures (e₁₆ in Figure 11B1, e₁₆ in Figure 12C2, e₁₅ in Figure 12E2). The same results can be obtained from the degree normalization histogram (Figure 13). Small layer-nodes have dramatic peaks of overwhelming magnitude. Therefore, these nodes are characterized as radical fighters that contribute significantly to the energy transfer. In other words, larger scale energy structures are stable participants, controlling the relatively monolithic stability of the turbulence system. Nevertheless, these small eddies disturb this stable system and promote the evolution of turbulent mixing.

Cascade is traditionally defined as a process in which energy is continuously transferred from one scale to the next in a decreasing order of magnitude Cardesa et al. (2015). However, the WSLN framework and its offshoots suggest that the mechanism of energy cascade differs significantly from that of homogeneous cascade models. And it suggests a synchronous multi-scale energy cascade pattern in which energy in one scale is synchronously transferred to all or part or none of the structures with smaller scales. Energy in node Layer2 (Figures 11A1, A2) is transferred to all layers of smaller scale. Energy in node Layer3 (Figures 11B1, B2) is transferred to all of the nodes with smaller scale. However, node Layer2 is only connected to node Layer5 in the first half of the 30 m depth segment (Figure 12A1) and node Layer2 is connected to node Layer3 and node Layer5 in the first half of the 80 m depth segment (Figure 12B1), indicating a large scale gap in energy transfer. Node Layer2 is connected to all smaller scales (Layer4, Layer5, Layer6) in Figure 12C2. However, the node Layer2 becomes isolated, and impedes energy transfer in the previous one (Figure 12C1). Similarly, Layer2 is adjacent to Layer4 and Layer5 in Figure 12E1, but it is connected to Layer5 and Layer6 in Figure 12E2, indicating a scale interruption. Furthermore, despite the similarity of the connections of the node pairs, the weight of the corresponding edges is completely different, implying the different intensity of energy transfer (Figures 11, 12), indicating the strong intermittency and inhomogeneity in the evolution of the energy cascade.

In addition, the network framework, based on marginal Hilbert spectra, covers both the inertial and dispersive regions of the turbulent power spectrum, allowing full identification of structures at all scales generated by the energy cascade. Meanwhile, syncretic signals are separated into isolated elements with relatively narrow frequency ranges and it characterizes full details of energy transfer between these structures of varying scales, providing a novel approach to turbulence analysis in the future.

The network structure coefficient κ is estimated to parameterize the energy transfer strength based on the specific layer scale and the number of inter-layer edges. The performance of is estimated with the microstructure profile. The variations of κ (shown in Figure 14) are in good agreement with the vertical profile of the dissipation rates. The shading represents the standard deviation of the dissipation rates. The correlation coefficients between κ and ϵ are 0.82 and 0.91 respectively, indicating a strong positive relationship. Thus, the ECMN framework is shown to be a valid structure for uncovering the underlying mechanism of the energy cascade, and κ is verified as an effective parameter for quantifying energy transfer. Furthermore, due to the strong stratification the κ estimates exhibit an abrupt increase around the pycnocline and are qualified to detect turbulent mixing oscillations. Thus, it shows generally good agreement in both strong mixing and well mixed regions, and the ECMN can be verified as an effective model for characterizing the underlying evolution of turbulent mixing.

Larger eddy decomposition and smaller vortex generation promote ocean circulation and turbulent mixing. The study of energy cascade is significant in oceanic studies. In the present work, the complex network analysis is used for the study of the chaotic interaction and scale property in energy cascade. An energy cascade multi-layer network (ECMN) was constructed, where nodes represent physical multivariate of turbulent fluid and inter-layer edges are indicators of energy transfer between scales. In the network, the underlying interactions of the energy cascade are transformed into topological elements, and the dynamic evolution of turbulent mixing can be revealed based on the topological property. In addition, the network framework covers both the inertial and dispersive regions of the turbulent power spectrum, allowing full identification of structures at all scales generated by the energy cascade. For a more detailed consideration of highly connected nodes, effective components in a multi-layer structure are aggregated into a single-layer network for both network analyses.

According to the topological framework of ECMN, the energy cascade process contradicts the homogeneous Richardson model, in which large eddies split into medium-sized eddies and gradually break down into small eddies, step by step. Energy is transferred between multi-scale structures. And the synchronous energy cascade pattern is demonstrated based on the topological characteristics of the network. It is also confirmed that the energy cascade process is chaotic, non-uniform and intermittent. We find that participants in energy transfer usually have different scale characteristics. Even when these participants have the same scale, the amount of energy exchanged is extremely different from each other. Meanwhile, the distribution of inter-layer edges to be irregular, indicating an intermittent evolution of the energy cascade. Furthermore, the large-scale vortices are shown to be stable structures that maintain the stability of the fluid flow. And these small eddies are the primary result of disturbance, facilitating turbulent mixing process.

We have also developed the network structure coefficient κ to assess the strength energy transfer. And it shows strong positive correlation between κ and intensity of turbulent mixing represented by the dissipation rate ϵ. The characteristics of the proposed network model is directly related to turbulent mixing. The results indicate that the proposed network model can reveal the chaotic property of energy cascade and effectively evaluate intermittent energy interaction underlying turbulent mixing.

The network framework offers novel insights for turbulence investigation and render the multi-layer network-based method efficient for analyzing chaotic evolution and inhomogeneous attributes in nonlinear and complex systems. Several obvious extensions, such as the direction of the energy cascade, are left for future studies. We believe that by considering the scale partitioning and focusing on the multi-scale interaction, we can conduct a complete research and propose a coherent link between turbulence and topology in the future.

The data presented in this study are available upon request from the corresponding author.

BM: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. HY: Conceptualization, Formal analysis, Funding acquisition, Project administration, Resources, Supervision, Validation, Visualization, Writing – review & editing. DS: Resources, Supervision, Writing – review & editing. JL: Formal analysis, Supervision, Validation, Writing – review & editing. WS: Resources, Writing – review & editing. XL: Funding acquisition, Writing – review & editing.