Keywords: turbulence, intermittence, multifractalité, le modèle SGC, loi de décroissance de l Probability Density Function xxxvi

Enfin, nous revisitons à l'aide du SGC la décroissance d'énergie d'une turbulence libre en tenant compte de l'intermittence. En raison de cette dernière, la mise en évidence du terme de rétrodiffusion d'énergie est plus complexe à démontrer, ainsi que ses impacts sur la loi de décroissance de l'énergie. Mais la phénoménologie reste la même, bien qu'avec des effets d'intermittence, par exemple, l'énergie est stockée à grande échelle par bouffées, et non plus de manière continue.

Résumé

La turbulence est l'un des problèmes fondamentaux non résolus de la physique classique, en dépit de sa manifestation dans de nombreux domaines, y compris l'ingénierie, par exemple dans l'énergie éolienne. Ceci est lié à notre manque de connaissance des propriétés des équations déterministes 3D de Navier-Stokes (NS) aussi fondamentales que l'existence et l'unicité de ses solutions. Cela n'empêche pas les chercheurs et les ingénieurs de l'utiliser. À l'aide de méthodes statistiques, le mécanisme de la turbulence a été partiellement révélé, comme le processus de transfert d'énergie. Par exemple, des modèles de fermeture de la turbulence, tels que le modèle quasi-normal à amortissement turbulent markovianisé (EDQNM), ont été introduits pour tenir compte partiellement de la hiérarchie infinie des équations des moments statistiques générée par le terme non linéaire des équations de Navier-Stokes. Ils ont mis en évidence la possibilité d'une rétrodiffusion de l'énergie du pic du spectre énergétique vers les plus gros tourbillons et donc une modification de la loi de décroissance de l'énergie de la turbulence. Cependant, ces avancées ne prennent pas en compte la caractéristique fondamentale de la turbulence : l'intermittence, qui signifie que la turbulence est extrêmement hétérogène et conduit à une grande divergence entre les preuves empiriques et ces modèles. Nous avons donc choisi le modèle déterministe Scaling Gyroscope Cascade (SGC), pour étudier la multifractalité de l'intermittence.

iii Le SGC est basé sur une discrétisation parcimonieuse de la forme de Bernoulli des équations de NS dans l'espace de Fourier qui préserve bien l'interaction de la triade d'un tourbillon parent et de ses tourbillons enfants, générant étape par étape une forte intermittence. Tout d'abord, les codes Python pour les trois méthodes de simulation numérique explicite -la méthode d'Euler, la méthode de Runge-Kutta du quatrième ordre et la méthode d'Adams-Bashforth asservie -sont présentés et testés afin de déterminer l'approche de simulation numérique la plus efficace pour le modèle SGC. Il ressort que la méthode d'Euler est la méthode de simulation numérique la plus efficace en comparant le temps d'exécution et la mémoire maximale. En outre, la structure spatiale du modèle SGC suggère que la complexité de calcul augmente de manière exponentielle avec le nombre d'étapes de la cascade.

Ensuite, l'intermittence du modèle SGC à de grands pas de cascade est étudiée en injectant différents forçages. L'existence de fluctuations spatio-temporelles est confirmée à l'aide d'une analyse statistique du flux d'énergie dans le domaine inertiel.

La distribution de probabilité de ces fluctuations présente des queues beaucoup plus épaisses que celles d'une distribution gaussienne. Pour obtenir un aperçu plus détaillé, l'analyse est poursuivie dans le cadre de l'Universal Multifractal (UM), basé sur des cascades multiplicatives stochastiques qui sont à la fois stables et attractives. Ces cascades sont déterminées par seulement quelques paramètres UM qui sont physiquement significatifs pour tout modèle de cascade, y compris le SGC.

Parmi les différents résultats obtenus, nous démontrons que l'indice clé de multifractalité est significativement inférieur à 2, remettant ainsi en cause le modèle log-normal encore souvent utilisé pour la turbulence hydrodynamique. from the energy spectrum peak to the largest eddies and therefore modification of the energy decay law of turbulence. However, these advancements do not take into account the fundamental characteristic of turbulence: intermittency, which means that turbulence is extremely heterogeneous and leads to a large discrepancy between the empirical evidence and these models. We therefore chose the deterministic Scaling Gyroscope Cascade (SGC) model, to investigate the multifractality of intermittency.

SGC is based on a parsimonious discretisation of the Bernoulli's form of the NS vi equations in Fourier space that well preserves the triad interaction of a parent eddy and its child eddies, generating step by step a strong intermittency.

Firstly, the Python codes for the three explicit numerical simulation methods -the Euler method, the fourth order Runge-Kutta method and the slaved Adams-Bashforth method -are presented and tested to determine the most efficient numerical simulation approach for the SGC model. It comes out that the Euler method is the most effective numerical simulation method by comparing the running time and maximum memory. Besides, the spatial structure of the SGC model suggests that the computing complexity increases exponentially with the number of cascade steps.

Then, the intermittency of SGC model at large cascade steps is investigated by injecting various forcings. The existence of spatial-temporal fluctuations is confirmed with the help of a statistical analysis of the energy flux in the inertial range. The probability distribution of these fluctuations has tails that are much heavier than those of a Gaussian distribution. To get more detailed insights, the analysis is pursued in the Universal Multifractal (UM) framework, based on stochastic multiplicative cascades that are both stable and attractive. These cascades are determined by only a few UM parameters that are physically meaningful for any cascade models, including the SGC.

Among the various obtained results, we demonstrate that the key multifractality index is significantly less than 2 thus questioning the log-normal model still often used for hydrodynamic turbulence.

Last but not least, we revisit with the help of SGC the energy decay of a free turbulence taking into account the intermittency. Due to the latter, the evidence of the energy backscatter term is more complex to demonstrate, as well as its impacts vii on the energy decay law. But the phenomenology remains the same, although with intermittency effects, e.g., energy is stored at large scales by puffs, no longer in a continuous manner. Turbulence occur everywhere from life to scientific research, such as the mixing of cream in a coffee cup, flight, and chemical experiments. Rapid development of technologies motivates the study of turbulence to predict and control it. For example, when designing the structural load and safety in marine engineering, one of factors causing the movement of sediments and resuspension is turbulence in the bottom layer of the ocean. Besides, turbulence will enhance sound scattering, making studying sound transportation in the ocean more complex. Turbulence in the atmosphere is more important than ever to many fields since recent research [START_REF] Williams | Intensification of winter transatlantic aviation turbulence in response to climate change[END_REF] predicts a high frequency of atmospheric turbulence occurring and it will strengthen with climate change to which turbulent heat flux may play an important role, though difficult to estimate [START_REF] Song | Turbulent heat flux reconstruction in the north pacific from 1921 to 2014[END_REF]. Turbulence-related injuries in the flight, according to the International Air Transport Association [START_REF] London | International air transport association[END_REF], reach 25%. On the other side, Gestionnaire du Réseau de Transport d'Electricité (RTE) claims that renewable resourcewind energy in France has remarkable growth [START_REF] France | Gestionnaire du réseau de transport d'electricité[END_REF], as wind energy in 2020 generated 7.9% of electricity and the maximum average rate of coverage of consumption in 2021 reached 31.35%. This renewable resource is the kinetic energy available in the large-scale movements of air in the atmosphere [START_REF] Burton | Wind energy handbook[END_REF], which is transported from the ground to the height and captured by wind turbines Fig 1 .1 1 . It requires quantifying turbulence caused by converting the kinetic energy of wind into thermal energy [START_REF] Manwell | Wind energy explained: theory, design and application[END_REF].

However, the mechanism of turbulence is neither completely understood nor precisely predictable [START_REF] Meneveau | Big wind power: seven questions for turbulence research[END_REF]. In the absence of computer technique, turbulence research at the beginning has been largely focused on the analytical approximations through deterministic and stochastic approaches, since the Navier-Stokes (NS) equation (See 1.2.1), which is the deterministic equation of turbulence, is still very challenging to guarantee and obtain the global existence of unique solutions. After the statistical method is proposed, significant achievements were achieved in the 20th century. The statistical description (See 1.2.2), which relies on the moments of velocity fluctuations from the NS equation, enables to study the dynamic of turbulence and gives generic properties, such as the energy spectrum and decaying energy. However, one of the challenges is an infinite hierarchy of enclosed moment equations caused by the nonlinear term, as the partial time derivative of moments n is determined by moments n + 1. In order to obtain a finite and closed set of equations, which is called closure scheme, this difficulty could be handed by truncating by a model.

Based on the assumption of homogeneous isotropic turbulence [START_REF] Taylor | Statistical theory of turbulence-ii[END_REF] and stationary statistics, the Kármán-Howarth [START_REF] De Karman | On the statistical theory of isotropic turbulence[END_REF] equations are introduced, which inspires statistic theory of turbulence and lays the foundation of "closure" approximations.

Kolmogorov proposed the well-known assumption of turbulence (See 1.4, 1.5) relying on the cascade structure (See 1.3). This assumption has been used to testify the following closure models.

The Quasi-Normal approximation (QNA) [START_REF] Millionshtchikov | On the role of the third moments in isotropic turbulence[END_REF] by assuming that the probability distributions of fluctuating velocity are normal distributions, closes the second and third-order moments equations. The fourth-order moments of the cumulant, which indicates the difference between the velocity distribution and normal distribution, are zero. Scientists have taken notice of this promising analytical approach because of its ability to achieve high Reynolds numbers and relative simplicity.

Kraichnan [START_REF] Kraichnan | Irreversible statistical mechanics of incompressible hydromagnetic turbulence[END_REF] noticed realizability problem of QNA that it violates the basic probabilistic inequalities and the energy spectrum obtained by the QNA model is unphysical. So he proposed a field-theoretic closure called Direct Interaction Approximation (DIA) to solve this realizability problem, which provides the exact statistical solution to a stochastic model and saves many properties shown in the NS equation.

The pity of DIA is that it is not consistent with Kolmogorov theory. Plus, Kraichnan [START_REF] Kraichnan | Inertial-range transfer in two-and three-dimensional turbulence[END_REF] investigated the energy transfer process and revealed the importance of nonlocal triad interaction (See 1.6) in two-dimensional turbulence that transfers most of energy from large scale to small scale [START_REF] Ohkitani | Triad interactions in a forced turbulence[END_REF][START_REF] Kida | Triad interactions in the dissipation range[END_REF]. The nonlocal triad interaction provides an approach to study the inner detail of the energy transferring in the scales of turbulence, but it questions Kolmogorov's theory due to the assumption of statistical independence of cascade motion.

Eddy Damped Quasi-Normal Markovian (EDQNM) Model [START_REF] Orszag | Analytical theories of turbulence[END_REF], a DIA-based closure model, is the currently most accepted closure model. By introducing a phenomenological eddy damping rate, EDQNM avoids unphysical features, especially the realizability guaranteed by Markovianization assumption, so that it is consistent with Kolmogorov theory. It has been proven as an effective model to study the dynamics and transport properties of turbulence [START_REF] Gréa | The effects of variable viscosity on the decay of homogeneous isotropic turbulence[END_REF]. In the framework of EDQNM, some deep insights have been proposed through the nonlocal interactions. One of them is the backscatter term [START_REF] Lesieur | Amortissement autosimilaire d'une turbulence à grand nombre de reynolds[END_REF] becomes a non-neglected factor in the energy decay of turbulence, which reveals the mechanism of inverse energy [START_REF] Métais | Statistical predictability of decaying turbulence[END_REF] and provides a more accurate energy decay process.

Unfortunately, these advances haven't taken one of the most fundamental phenomena of natural turbulence: intermittency [START_REF] Batchelor | The nature of turbulent motion at large wave-numbers[END_REF][START_REF] Frisch | Intermittency in nonlinear dynamics and singularities at complex times[END_REF] (See 1.7) into account, resulting in strongly non-Gaussian fluctuations over a wide range of space-time scales. Even though the traditional assuring scale invariance method by neglecting intermittency can describe turbulence energy transport, which involves only the second-order moments, description for the higher-order statistical moments urgently needs to be cor-

Turbulence

Turbulence are disorder flows over time and space domain, and the motion is random and believed to be highly unpredictable. There is no precise mathematical definition, but a physical definition was first proposed in Reynolds experiment [START_REF] Reynolds | Xxix. an experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels[END_REF].

In this experiment, the flow was pumped into a long, straight circular pipe to get rid of vibration from the left side, which is shown in Fig 1. where L is the scale of motion which is D in Reynolds experiment. The critical value R c suggested by Reynolds was the order of 2000. For R > R c , the mixture of colored stream and surrounding water on the right side is irregular and random, which is turbulent flow. On the contrary, the unmixed stream is laminar flow.

Re = U L v ∼ T he inertial f orce T he viscosity f orce (1.2.1)

The deterministic equation: NS equation

It's generally accepted that the equation governing turbulence is Navier-Stokes equation [START_REF] Navier | Mémoire sur les lois du mouvement des fluides[END_REF][START_REF] Bistafa | On the development of the navier-stokes equation by navier[END_REF]:

∂u ∂t + (u • ∇)u = v∇ 2 u + f - 1 ρ ∇p (1.2.2)
where u is velocity vector which depends on time t and location vector r; p is pressure; ρ is pressure; f is body force; ∇ represents gradient. For incompressible fluid, the divergence of velocity is

∇ • u = 0. (1.2.3)
It is of great importance that it characterizes the constant transfer of physical quantities in the fluid. The unsteady-state term, which is the first item on the left side of Eq 1.2.2, is the rate of changed momentum of fluid per unit volume. The advection term, which is the second item on the left side of Eq 1.2.2, represents the momentum that flows in and flows out of control surface. The right side of Eq 1.2.2 is pressure term, external forcing, and viscous dissipation, which contributes to momentum of dissipation of Newtonian viscous fluids. Unfortunately, the existence of rigorous mathematical unique solution for incompressible Navier-Stokes equation hasn't been obtained so far due to its complexity [START_REF] Robinson | The three-dimensional Navier-Stokes equations: Classical theory[END_REF] resulting from the non-linear term and infinite-dimensional problem.

NS equation follows three basic physical laws: the conservation of mass, Newton's second law, and energy conservation. Energy transfer rate [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF] is defined as the energy per mass per time that passes through a spherical shell in Fourier space:

ε = - 1 2 ∂ ∂t (u • u). (1.2.4)
It is conserved by the nonlinear term of NS equation assuming v = 0, f = 0 and density is constant.

Moreover, NS equation is symmetrical so that many properties are invariant under these conditions. These symmetries are:

1. Space-translations: t, r, u → t, r + l, u, l ∈ R 3 2. Time-translations: t, r, u → t + τ, r, u, τ ∈ R 3. Galilean transformations: t, r, u → t, r + U t, u + U , U ∈ R 3
4. Parity: t, r, u → t, -r, -u 5. Rotations: t, r, u → t, Rr, Ru, R ∈ SO(3, R)

6. Scaling: t, r, u → e 1-a t, er, e a u, e ∈ R + , a ∈ R

The first four symmetries can be explained as Galilean invariance and scaling symmetry is obvious under the limit v → 0, f → 0. A hidden scaling symmetry [START_REF] Mailybaev | Hidden scale invariance in navier-stokes intermittency[END_REF] of NS equation (v → 0) is introduced from rescaling around suitably defined Lagrangian scaling centres, which is different from Galilean invariance.

Statistical approach

When Reynolds number of NS equation is high, the degrees of freedom are huge and turbulence is very sensitive to small disturbances and errors. It's impossible to trace and predict turbulence because of its random instantaneous motion. After

Taylor [START_REF] Taylor | Mechanism of the production of small eddies from large ones[END_REF] suggested that the motion of turbulence can be studied as the stochastic process, statistical method [START_REF] Krommes | Fundamental statistical descriptions of plasma turbulence in magnetic fields[END_REF] is introduced as one of the important research approaches to quantify turbulence, e.g., the averages, probability distribution functions, spectra, correlations. One of statistical methods is correlation function:

R ij (l) = ⟨u i (r) • u j (r + l)⟩. (1.2.5)
This correlation function describes the relationship of velocity at different locations.

Kármán and Howarth in 1938 [START_REF] De Karman | On the statistical theory of isotropic turbulence[END_REF] proposed a self-similar law for this spatial correlation function of velocity during studying decay of three-dimensional isotropic turbulence at a high Reynolds number.

Many physical quantities of turbulence are related to the correlation function. For example, energy spectrum [5] which is under assumption of homogeneous, describes the distribution of average kinetic energy wave number space k:

R ii (l) = ∞ -∞ e -ikl E(t)dk ⇒ E(k) ≡ 1 2 |k| ϕ ii (k)dσ(k) (1.2.6)
where ϕ ii (k) =

+∞

-∞ e ik•r R ii (r)dr is three-dimensional. To avoid the difficulty to determine probability density function with few parameters, the moment of high order q [3] providing approximation is:

⟨[u(r + l) -u(r)] q ⟩ = ⟨(δ r u(l)) q ⟩ = S q (l), (1.2.7)
where

δ r u(l) = u(r + l) -u(r) (1.2.8)
is velocity increment. Assuming that turbulence are maintained at a steady state by external forces, the right sides of Eq 1.2.7 under assumption of homogeneous and isotropic is independent of time, direction, and selected spatial point r, so it is simplified to a velocity component relating only to the magnitude of distance l, which is named as structure function:

S q (l) ∝ l ζ(q) , (1.2.9)
where ζ(q) is the exponent of the qth-order moment structure function.

Fruitful conclusions were built on assumption of statistical theory, such as the well-known K41 theory. What's more, strongly anisotropic systems are studied by Generalized Scale Invariance (GSI) [START_REF] Schertzer | Generalised scale invariance in turbulent phenomena[END_REF][START_REF] Lovejoy | Generalized scale invariance in the atmosphere and fractal models of rain[END_REF]. It's actually self-affine generated by the pullback transforms and differential dilations/contractions operators. In two dimensions, the operators could be represented by off-diagonal matrices.

Energy cascade

Richardson [START_REF] Richardson | Weather prediction by numerical processes[END_REF] is considered that laid the foundation for turbulence structure as he first introduced the energy cascade to describe scaling law of turbulence. Part proposed a statistical theory of turbulence, which measures statistical quantities of turbulence. Moreover, achievements in multiplicative cascade models suggest the importance of the cascade is far more significant.

Kolmogorov assumption

Here, the well-known Kolmogorov assumption, which quantitatively describes energy distribution between scales and is considered as the fundamental advancement of turbulence research, is explained. To start, several hypotheses are proposed to get rid of large anisotropic eddies, which are affected by the boundary condition.

Lemma 1

The hypothesis of local isotropy: When Reynolds number is high enough, the motion of turbulence at small scales l < l I is statistically isotropic, where l I is the scale dividing isotropic eddies from those anisotropic large eddy.

ϵ ∝ u 2 τ = u 2 l/u = u 3 l , (1.4.1)
where τ is eddy turnover time. Eq 1.4.1 yields that the mean square of velocity increment between two points is two-thirds powers of the distance.

Lemma 3 Structure function: The second-order structure function depends on length scale l and mean dissipation rate ϵ:

S 2 (l) ∝ l 2H * ∼ (ϵl) 2/3 , (1.4.2)
where H * is a measure of correlation.

Therefore, famous -5/3 energy spectrum [START_REF] Obukhov | On the distribution of energy in the spectrum of turbulent flow[END_REF] is obtained from the definition of energy spectrum Eq 1.2.6.

Lemma 4 Kolmogorov-Obukhov spectrum: Scaling law of spectrum in inertial range is:

E(k, t) ∝ ϵ 2/3 k -β * = ϵ 2/3 k -5/3 , (1.4.3)
where β * is power spectrum exponent.

Besides, relationship between the second-order structure function and energy spectrum is

β * = 1 + ζ(2) ≈ 1 + 2H * , (1.4.4)
where ζ( 2) is the second-order-moment structure function exponent.

Thesis -X. LI; Ecole des Ponts ParisTech -Hydrology Meteorology et Complexity

Variety of experiments [START_REF] Saddoughi | Local isotropy in turbulent boundary layers at high reynolds number[END_REF][START_REF] Grant | A spectrum of turbulence at very high reynolds number[END_REF] and theories [START_REF] Calif | Modeling of atmospheric wind speed sequence using a lognormal continuous stochastic equation[END_REF] since the 1950s have confirmed that the micro-structure of turbulence is what predicted by Kolmogorov and Obukhovthe second-order structure function of turbulence obeys 2 3 law, and the one-dimensional spectrum of turbulence obeys -5 3 power law.

Energy decay law

The dynamics of turbulence are generally referred to energy transferring from large eddies to small eddies, which maintains motion of turbulence. Without external forces, the kinetic energy (TKE) of fully developed turbulence will be consumed by the viscosity term. Therefore, free energy decay law is one of the key characteristics to reveal the mechanism of turbulence, which could be widely applied to rapid expansion of wind energy [START_REF] Hvelplund | Renewable energy and the need for local energy markets[END_REF]. Until now, plenty of research on the free energy decay laws has been carried out. Even though there is a general agreement on the kinetic energy of freely decaying turbulence that decays with time, the exponent of energy decay law is still debatable.

According to the definition of energy spectrum, kinetic energy is

K ≡ 1 2 R ii (0) = +∞ 0 E(k, t)dk. (1.5.1)
Kármán and Howarth [START_REF] De Karman | On the statistical theory of isotropic turbulence[END_REF] firstly suggest energy decay is self-similar for scales outside the dissipation range. Batchelor and Townsend [START_REF] Batchelor | Decay of isotropic turbulence in the initial period[END_REF] investigated the decay of homogeneous isotropic three-dimensional turbulence at a high Reynolds number and predicated energy decay law in the initial period when viscous dissipation and inertia forces are of comparable importance, is

K(t) ∝ t -1 , (1.5.2)
whereas energy decay law in the final period after which inertia forces are negligible is

K(t) ∝ t -5 2 .
(1.5.3)

Comte-Bellot and Corrsin [START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of grid-generated turbulence[END_REF] obtained that energy decay obeys a power law

K(t) ∝ t αe (1.5.4)
where exponent α e is less than -1. George in 1992 [START_REF] George | The decay of homogeneous isotropic turbulence[END_REF] confirmed this conclusion that exponent α e should be between -5 2 and -1 and explicit value is decided by the behavior of energy spectrum near k → 0. By assuming Loitsiansky integral

I(t) = -r 2 R ii (r, t)dr (1.5.5)
is an invariant, Kolmogorov [START_REF] Kolmogorov | On degeneration (decay) of isotropic turbulence in an incompressible viscous liquid[END_REF] derived

K(t) ∝ t -10 7 . (1.5.6)
Besides, the freely decaying turbulence with this assumption is commonly known as Batchelor turbulence.

However, Proudman and Reid in 1954 [START_REF] Proudman | On the decay of a normally distributed and homogenous turbulent velocity field[END_REF] and Batchelor and Proudman in 1956 [START_REF] Batchelor | The large-scale structure of homogenous turbulence[END_REF][START_REF] Métais | Statistical predictability of decaying turbulence[END_REF] proved that Loitsiansky scale isn't an invariant, which brings questions to -10 energy decay law. Saffman [START_REF] Saffman | The large-scale structure of homogeneous turbulence[END_REF] therefore proposed that the second integral moment of velocity correlation is an invariant:

L s (t) = R ii (r, t)dr ̸ = 0, (1.5.7)
which is named as Saffman invariant, and the exponent α e associated with this hypothesis which is referred as Saffman turbulence is -6 5 . Both measurement data and experimental data of freely decaying grid-generated turbulence [START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of grid-generated turbulence[END_REF][START_REF] Krogstad | Is grid turbulence saffman turbulence?[END_REF][START_REF] Morize | Energy decay of rotating turbulence with confinement effects[END_REF] have verified it.

Actually, energy decay process is dominated by contributions from small wave numbers as a consequence of exponential factor in the integrand, so its power-law decay during the final period results from a behavior of E(k, 0) (k → 0). Energy decay law α e = - 10 7 corresponds to the energy spectrum E(k, 0) ∼ k 4 at large scale, while energy decay law α e = - 6 5 is related to spectrum E(k, 0) ∼ k 2 . Meanwhile, integral scale l measuring the correlation distance between two terms of distance or time has a power law:

L(t) = ∞ 0 k -1 E(k, t)dk ∞ 0 E(k, t)dk ∝ t α l . (1.5.8)
Exponent α l in Batchelor turbulence is 2 7 obtaining from constant u 2 l 5 = const, whereas it is 2 5 in Saffman turbulence resulting from constant u 2 l 3 = const. More detailed research on energy decay process will be discussed in Chapter 5.

The nonlocal triad interaction

Although a clear big picture of turbulence in which energy spectrum and energy decay obey self-similar law has been obtained, the microcosm mechanism dominated energy transferring in the scales of turbulence is also of great importance. Actually, triads interaction preserves detailed conservation of kinetic energy that provides significant information about the way nonlinear interactions occur between these modes. Fourier analysis of NS equation, an infinite set of ordinary differential equations, could be projected into a plane perpendicular to the wave vector due to incompressibility eliminating longitudinal motion. Thus, the quadratic nonlinearity decomposes into three Fourier modes interactions, which is named as triad interaction by Kraichnan [START_REF] Kraichnan | Irreversible statistical mechanics of incompressible hydromagnetic turbulence[END_REF].

Velocity u(x, t) of homogeneous turbulence in Fourier space is:

u(k, t) = 1 (2π) 3 u(x, t)e
κ • u = 0, (1.6.3)
where κ is wave number. When any vector V decomposes into:

V = V ∥ + V ⊥ , (1.6.4) where component V ∥ is parallel to κ and component V ⊥ is normal to κ. Then component of V ∥ becomes V ∥ = e(e • V ) = κ(κ • V )/κ 2 ⇐⇒ V ∥ j = κ j κ k κ 2 V k , (1.6.5)
where the unit vector e equals e = κ/κ. Unit vector Due to Eq 1.6.4, component

V ⊥ has V ⊥ = V -κ(κ • V )/κ 2 ⇐⇒ V ⊥ j = P ik V k , (1.6.6)
where P ij is projection tensor allowing that component Eq 1.6.6 is the projection of V onto the plane normal to κ:

P ij (κ) ≡ δ jk - κ j κ k κ 2 .
(1.6.7)

Projection tensor

Regarding momentum conservation on an fluid element, it yields [START_REF] Zhou | Turbulence theories and statistical closure approaches[END_REF]

∂u j ∂t + ∂u j u k ∂x k = v ∂ 2 u j ∂x k ∂x k - 1 ρ ∂p ∂x j , (1.6.8)
the equation for time evolution of velocity vector is obtained by its Fourier transform:

d u j dt + vκ 2 u j = -iκ j p -V j (1.6.9)
where V j is derived from Fourier transform of nonliear term of Eq 1.6.8. Due to Eq 1.6.3, V j has:

κ 2 p = iκ j V j ⇒ V ∥ j = κ j κ k κ 2 V k = -iκ j p, (1.6.10)
which indicates that -V ∥ j is balanced by pressure term. Hence, Eq 1.6.9 turns into:

d u j dt + vκ 2 u j = -(δ ij - κ j κ k κ 2 ) V k = -P ij V k = -V ⊥ j (1.6.11)
Now it's time to give out V j :

V j (κ, t) = iκ k κ ′ u j (κ ′ ) u k (κ -κ ′ ),
(1.6.12) Eq 1.6.11 becomes

( d dt + vκ 2 ) u j (κ, t) = -iκ l P jk (κ) κ ′ u k (κ ′ , t) u k (κ -κ ′ , t). (1.6.13)
It shows time evolution of velocity u(κ, t) are related to velocities at wave number p = κ ′ and q = κκ ′ (κ = k). These three wave vectors (k, p, q) in Fig 1 .5 forms a triangle k + p + q = 0.

(1.6.14)

Hence, the elementary interactions between this triad (k, p, q) resulting from the quadratic nonlinearity conserves energy.

q ≪ p ∼ k, (1.6.15a) p ≪ q ∼ k, (1.6.15b) k ≪ p ∼ q.
(1.6.15c)

The nonlocal triad interaction makes energy locally exchanged between two larger wave numbers. Plus, the nonlocal triad interaction in Eq.1.6.15c could be represented by k-independent viscosity v turb (k) [START_REF] Kraichnan | Eddy viscosity in two and three dimensions[END_REF], which has widely applied to the subgrid models for homogeneous turbulence. Kraichnan [START_REF] Kraichnan | Inertial-range transfer in two-and three-dimensional turbulence[END_REF] investigated energy transfer process and confirmed the importance of nonlocal triad interaction, which has been demonstrated by more detailed studies [START_REF] Domaradzki | Analysis of energy transfer in direct numerical simulations of isotropic turbulence[END_REF][START_REF] Domaradzki | Local energy transfer and nonlocal interactions in homogeneous, isotropic turbulence[END_REF] through the DNS of Navier-Stokes turbulence and motives further research of energy transfer process [START_REF] Domaradzki | Nonlocal triad interactions and the dissipation range of isotropic turbulence[END_REF][START_REF] Zhou | Interacting scales and energy transfer in isotropic turbulence[END_REF][START_REF] Fouquet | Evolution of high reynolds number two-dimensional turbulence[END_REF][START_REF] Ohkitani | Triad interactions in a forced turbulence[END_REF].

For instance, the parameter λ k to separate the local and nonlocal interactions [START_REF] Lesieur | Turbulence in fluids: stochastic and numerical modeling[END_REF] max(k, p, q) min(k, p, q) ≥ λ k , (1.6.16)

and the parameter that has been used by Lesieur is 2. Moreover, the backscatter term and intermittency based on nonlocal interactions will be discussed later.

The intermittency

Further investigation showed that structure function proposed by K41 theory, as well as scaling laws of energy dissipation rates at higher order moments, don't agree with experimental results. These significant deviations are caused by intermittency, which breaks self-similarity properties [START_REF] Frisch | Fully developed turbulence and intermittency[END_REF] and brings the scale invariance under doubt. For example, one of the shortcomings of K41 theory was that energy dissipation rate ϵ is uniform in space and constant in time. Moreover, it is not only the fundamental characteristic of natural turbulence [START_REF] Schertzer | The dimension and intermittency of atmospheric dynamics[END_REF][START_REF] Douady | Direct observation of the intermittency of intense vorticity filaments in turbulence[END_REF] but also a general phenomenon in geophysics [START_REF] Schertzer | Nonlinear variability in geophysics: Multifractal simulations and analysis[END_REF], raising doubts about previous research.

The general concept of intermittency is that high activate part of field concentrates Intermittency was first introduced by Batchelor and Townsend [START_REF] Batchelor | The nature of turbulent motion at large wave-numbers[END_REF] to describe spatial inhomogeneous energy caused by the stretching of vortex filaments that the fraction of available space for energy reduces as scale of eddies decreases. This spottiness of eddies whose scale is small is the internal intermittency of turbulence. Landau [START_REF] Landau | Mechanics of fluids[END_REF] noticed that intermittency should be taken into consideration in energy transferring process and claimed dissipation rate ε is a function of space and time. Later Kolmogorov [START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number[END_REF] refined K41 theory and ε in K62 hypothesis is assumed to have a log-normal distribution

σ 2 logε ∼ A + alog(L/l) (1.7.1)
where a is some universal constant, σ 2 logε is the variance and L is the largest scale.

The second-order moment scaling function now is

S 2 (l) ∝ ε 2/3 l 2/3 (1.7.2)
which is also known as refined similarity hypothesis.

More research about intermittency will be explained in Chapter 2.

Summary

Turbulence is unquestionably complicated, and given how important it is, its mechanism and properties have to be exposed. Some basic knowledge and significant advancements of turbulence that have been described above provide a brief comprehension of the dynamic of turbulence.

The fractal geometry

As well known, fractal theory characterizing the simplest scale invariance-similarity degree, laid the foundation to overcome the challenge faced by traditional geometry in describing complex irregular events in the nature, and becomes a powerful tool for studying physical phenomena. It motivates the β-model that take intermittency into consideration.

Despite the fact that many studies significantly boosted the development of fractal theory [START_REF] Hausdorff | Der wertvorrat einer bilinearform[END_REF], fractal came to light after Mandelbrot [START_REF] Mandelbrot | How long is the coast of britain? statistical self-similarity and fractional dimension[END_REF][START_REF] Richardson | Atmospheric diffusion shown on a distance-neighbour graph[END_REF] explained Richardson's empirical finding that the length of coast will be greatly exceeded under some dimensional measurements, as the lengths are different under different scale measurements.

When the scale approaches 0, the length will be positive infinity. Therefore, he gave

Hausdorff dimension [START_REF] Hausdorff | Dimension und äußeres maß[END_REF] based on self-similar method and named as fractal dimension D F . Fractal dimension D F is the measure of geometric irregularity degree. Fractal dimension of the Britain coastline is 1.25. It is less than Euclidean dimension of the plane 2 and larger than Euclidean dimension of straight line 1. Moreover, fractal dimension describes self-similarity that is repeated at each scale. The fractal geometry allows to quantify the complexity of a chaotic system by fractal dimension of the converged attractor of trajectories in their phase space. Moreover, it has been widely applied to model the natural phenomena, for example, topographical surfaces, wind activity [START_REF] Chang | Fractal dimension of wind speed time series[END_REF], and rainfall fields.

A fractal set F [START_REF] Falconer | Fractal geometry: mathematical foundations and applications[END_REF] satisfies the followings:

(i) F has a fine structure, i.e. details on arbitrarily small scales;

(ii) F is too irregular to be described in traditional geometry, both locally and globally;

(iii) F has some form of self-similarity, perhaps approximate or statistical;

(iv) F in most cases is defined in a very simple way, perhaps recursively. There are several methods to estimate fractal dimension D F , such as the box counting method. The most commonly theoretical definition is given as below [START_REF] Mandelbrot | How long is the coast of britain? statistical self-similarity and fractional dimension[END_REF][START_REF] Mandelbrot | The fractal geometry of nature[END_REF]:

D F = lim r→0 - logN log r(N ) ⇐⇒ N = r(N ) -D F (2.1.1)
where N is the number of parts that object exactly decomposed into, and r(N ) is the similarity ratio that divides object or inverse resolution. In this way, fractal dimension D F is no longer an integer number. And a large fractal dimension D F indicates the structure is highly irregular.

Another decisive parameter in fractal geometry is co-dimension c, which is the intermittency exponent or Hurst exponent:

c = D -D F (2.1.2)
where D is Euclidean dimension. Contrary to fractal dimension D F , the co-dimension c describes the degree of sparseness.

Here, one-dimensional fractal set Cantor set Fig. which suggests that the co-dimension is c ≈ 0.37.

D F = lim n→∞ - log(N (l n )) log( 1 λ n ) = log2 log3 ≈ 0.
Typically, the difficulty of directly obtaining fractal dimension requires estimating the co-dimension to describe irregular properties of fields.

Multiplicative Process

Following L.F.Richardson, the notion of multiplicative cascade was developed in statistical study of turbulence by Kolmogorov as a phenomenological framework to consider intermittency. Yaglom [START_REF] Yaglom | The influence of fluctuations in energy dissipation on the shape of turbulence characteristics in the inertial interval[END_REF] exhibits an explicit cascade to take intermittency into account. It's a discrete procedure in which the breakdown of turbulent eddies occurs, which is named as "self-similar breakdown of turbulence eddies". He proposed that the distribution of ε is log-normal distribution. This result also testifies to one of the shortcomings of K41 theory. Although log-normal approximation is considered as a good solution for the problem of intermittency in hydrodynamic turbulence, its application has been questioned by other researches. The most common challenge [START_REF] Mandelbrot | Possible refinement of the lognormal hypothesis concerning the distribution of energy dissipation in intermittent turbulence[END_REF][START_REF] Mandelbrot | Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier[END_REF] is that tails of probability distribution generated by nonlinear cascade process [START_REF] Nelkin | Scaling theory of hydrodynamic turbulence[END_REF][START_REF] Kraichnan | Convection of a passive scalar by a quasi-uniform random straining field[END_REF] deviate from the log-normal distribution, especially for highorder moments. And it was until 1980, the fact that these multiplicative cascade models [START_REF] Yaglom | The influence of fluctuations in energy dissipation on the shape of turbulence characteristics in the inertial interval[END_REF][START_REF] Novikov | The intermittency of turbulence and the spectrum of energy dissipation fluctuations[END_REF][START_REF] Gurvich | Breakdown of eddies and probability distributions for small-scale turbulence[END_REF] actually is general multifractal process, is realized.

Energy transfers from the large scale l 0 to small scale l by successive steps. Separate ratio in the whole process is assumed as a fixed value λ. After separating n steps, the length of small scales is

l n = l 0 λ n (2.2.1)
Energy transfer rate ε ln at scale l n should be calculated by the multiplicative method as shown in Fig. 2.3:

ε n = ε ln = µ n • ε n-1 . (2.2.2)
where µ n is usually assumed to be a scale-independent increment µ.

ε n = µ n µ n-1 • • • µ 1 ε 0 = ( n m=1 µ m )ε 0 = µ n ε 0 (2.2.3)
Since turbulence energy is conserved in inertial range, the total transfer rate is conserved

⟨ε n ⟩ = ⟨ε 0 ⟩. (2.2.4)
Besides, the scaling moment of energy flux density ε n is:

⟨(ε n ) q ⟩ = n m=1 ⟨(µ m ) q ⟩⟨ε 0 ⟩ q = ⟨µ q 1 ⟩ n ⟨ε 0 ⟩ q (2.2.5)

β-model

Based on Novikov-Stewart model [START_REF] Novikov | The intermittency of turbulence and the spectrum of energy dissipation fluctuations[END_REF], β-model [START_REF] Frisch | A simple dynamical model of intermittent fully developed turbulence[END_REF] which is one of the fundamental cascade models considering intermittency is introduced. It has scaling law associated with fractal dimension and therefore provides more physical insight than the traditional approach based on probabilistic models of the dissipation. The idea is that part of the total space in scales is excited.

In the β-model, as shown in

Thesis -X. LI; Ecole des Ponts ParisTech -Hydrology Meteorology et Complexity

The fraction of alive eddies at scale l n is

β n = β n , (2.3.1)
where β is the fraction of alive eddies decreasing by each step. This is the reason why it's named as β-model.

According to Eq.2.1.1, the alive fraction β n can be expressed in fractal form:

β n = N N (l n ) = r(N ) -D F r(N ) -D = λ n•D F λ n•D = λ n•(-c) . (2.3.2)
The alive fraction β at each step from Eq 2.3.1 is

β = λ -c . (2.3.3)
Considering this model in the probabilistic form, there are two states of the random increments µ. The probability of alive increment is:

P r(µ ̸ = 0) = λ -c . (2.3.4)
The probability of dead increment is:

P r(µ = 0) = 1 -λ -c . (2.3.5)
Here, the nonzero increment

µ = λ c , (2.3.6)
which is assumed from the alive fraction β at each step, is larger than 1 to guarantee the conservation of ε in whole process:

⟨µ⟩ = λ c • P r(µ = λ c ) + 0 • P r(µ = 0) = 1 ⇐⇒ ⟨ε n ⟩ = ⟨ε 0 ⟩. (2.3.7)
After iterating n steps, there are two possible states at solution r n = λ n . The probability of alive state is:

P r(ε n = µε n-1 = (r n ) c = λ nc ) = λ -nc . (2.3.8)
The probability of dead state is: Due to the length of scales, there are two kinds of scale changes: up-scaling and down-scaling. The length of scales changes with a fixed number that is generally related to the scale ratio λ.

P r(ε n = µε n-1 = 0) = 1 -λ -nc . ( 2

The probabilistic conservation of upscaled β-model

The up-scaling process is to study the field whose resolution is Λ at the smaller resolution Λ λ n until reaching 1. As an illustration, a two-dimensional β-model is upscaled. The field in Fig 2 .6a has a resolution of Λ = 2 5 . The resolution of upscaled field in Fig 2 .6b is 2 3 . Each event in the up-scaled field consists of all events at maximum resolution in the same location so that the mean of these events ε Λ at maximum resolution Λ is the event in the up-scaled field. In order to simplify writing, A n represents the survival eddy at nth step whose

energy flux density is ε n = λ nc A n : ε n = λ nc , (2.3.10)
while A c n is complementary event of A n corresponding to the dead eddy at nth step whose energy flux density is ε n = 0

A c n : ε n = 0. (2.3.11)
The number of alive eddies at nth step is denoted as r and s is the number of dead eddies. Eddies at nth step are independent. Upscaled steps is m.

Therefore, A r n (A c n ) s
is defined as possible upscaled event at (nm)th step, which is generated by r alive eddies and s dead eddies at nth step (r +s = 2 m ). The energy flux density of

A r n (A c n ) s is A r n (A c n ) s : ε n-m = r • ε n (A n ) + s • ε n (A c n ) 2 m = r • 2 nc 2 m , ∀r + s = 2 m . (2.3.12)
After upscaled one step (m = 1), there are three types of upscaled states:

A 2 n : ε n-1 = 2 nc ; (2.3.13a) A n A c n : ε n-1 = 2 nc 2 ; (2.3.13b) (A c n ) 2 : ε n-1 = 0. (2.3.13c)
According to the up-scaling process and independence, their probabilities are:

P (A 2 n ) = P (A n ) 2 ; (2.3.14a) P ((A c n ) 2 ) = P (A c n ) 2 ;
(2.3.14b)

P (A n A c n ) = 2P (A n )P (A c n ). (2.3.14c)
Here, upscaled event A n A c n could be produced by two possible arrangements of eddies at nth step: (A n , A c n ), (A c n , A n ), so the coefficient of its probability in Eq 2.3.14c should be 2.

Moreover, the sum of upscaled events' probabilities equals 1, which ensures probabilistic conservation in upscaling process, as shown in following:

P (A 2 n ) + P (A n A c n ) + P ((A c n ) 2 ) = (P (A n ) + P (A c n )) 2 = 1. (2.3.15)
For other cases (m > 1), the number of possible upscaled event A r n (A c n ) s increases and probability can be directly derived by the binomial coefficient theory, which indicates possible arrangements of eddies at nth step should be taken into consideration.

Therefore, probability of upscaled event A r n (A c n ) s is given by:

P (A r n (A c n ) s ) = ( 2 m r )P (A n ) r P (A c n ) s . (2.3.16)
The possibility field of up-scaled events is also conserved:

r+s=2 m P (A r n (A c n ) s ) = r+s=2 m ( 2 m r )P (A n ) r P (A c n ) s = (P (A n ) + P (A c n )) 2 m = 1. (2.3.17) Fig 2.
7 is the probability of events in two cascade step upscaled field. Energy flux density is normalized by

ε n = λ nc so that event 1 in Fig 2.7 is A n and event 0 is A c n .
Figure 2.7: The probability of events in two cascade step upscaled field. The original field is one-dimensional β-model whose resolution is 2 8 . The codimension is 0.2. Now the statistical moments of up-scaled event ε n-m could be calculated as follow:

E(ε q n-m ) = r+s=2 m (ε n-m (A r N (A c N ) s )) q P (A r N (A c N ) s ) = r+s=2 m ( r 2 m • 2 nc ) q ( 2 m r )P (A N ) r P (A c N ) s . (2.3.18)
For the case q = 1:

E(ε n-m ) = 2 m r=1 r 2 m • 2 nc ( 2 m r )P (A N ) r P (A c N ) s = 2 m r=1 ( 2 m -1 r-1 )P (A N ) r-1 P (A c N ) 2 m -r = (P (A N ) + P (A c N )) 2 m -1 = 1. (2.3.19)
The high statistical moments in study to describe distribution have q ≤ 4. The upscaled step m is assumed large than 2 and it yields q < 2 m . From Eq 2.3.19, the statistical moments function can be substituted into:

E(ε q n-m ) = 2 m r=1 ( r 2 m • 2 nc ) q-1 ( 2 m -1 r-1 )P (A N ) r-1 P (A c N ) s = 2 m r=1 ( r 2 m ) q-1 ( 2 m -1 r-1 )P (A N ) r-q P (A c N ) s = 2 m r=1 ( r 2 m ) q-1 • (2 m -1)! (2 m -r)! • (r -1)! P (A N ) r-q P (A c N ) s .
(2.3.20)

According to the inequality

1 2 m < 1 2 m -1 < 1 2 m -2 , (2.3.21)
Eq 2.3.20 becomes

E(ε q n-m ) < 2 m r=1 r q-1 • (2 m -q)! (2 m -r)! • (r -1)! P (A N ) r-q P (A c N ) s < (2 m -q)! (2 m -1)! P (A N ) 1-q P (A c N ) 2 m -1 + 2 m r=2 r q-1 • (2 m -q)! (2 m -r)! • (r -1)! P (A N ) r-q • P (A c N ) s . (2.3.22)
The right hand side of Eq 2.3.22 is simply denoted as E max (ε q n-m ). Considering

r = 1 + r -1, E max (ε q n-m ) turns into E max (ε q n-m ) = (2 m -q)! (2 m -1)! P (A N ) 1-q P (A c N ) 2 m -1 + 2 m r=2 r q-2 • (1 + r -1) • (2 m -q)! (2 m -r)! • (r -1)! • P (A N ) r-q P (A c N ) s = (2 m -q)! (2 m -1)! P (A N ) 1-q P (A c N ) 2 m -1 + 2 m r=2 r q-2 • (2 m -q)! (2 m -r)! • (r -1)! P (A N ) r-q • P (A c N ) s + 2 m r=2 r q-2 • (2 m -q)! (2 m -r)! • (r -2)! P (A N ) r-q P (A c N ) s . . . = (2 m -q)! (2 m -1)! P (A N ) 1-q P (A c N ) 2 m -1 + 2 m r=2 r q-2 • (2 m -q)! (2 m -r)! • (r -1)! P (A N ) r-p • P (A c N ) s + • • • + 2 m r=q (2 m -q)! (2 m -r)! • (r -q)! P (A N ) r-p P (A c N ) s .
(2.3.23)

Obviously the last term in Eq 2.3.23 equals 1 and

E(ε q n-m ) < E max (ε q n-m ), E max (ε q n-m ) > 1. (2.3.24)
When statistical moment is q = 1, the upscaling process of β-model is conserved. When the statistical moments q is larger than 1, the value of upscaled fields

E max (ε q n-m
) suggests the weight of the alive eddies in the up-scaled field gets higher.

The probabilistic conservation of β-model in the downscaling

This part will discuss the conservation of down-scaled field whose probability is determined by the hierarchical path in downscaling.

The ancestor eddy must be alive ε 0 = 1 due to the fact that alive eddies cannot be generated by dead mother eddy. The probabilities of increments are simplified as:

P r(µ = λ c ) = λ -c = P 1 ; (2.3.25a) P r(µ = 0) = 1 -λ -c = 1 -P 1 . (2.3.25b)
The number of alive eddies at nth step is denoted by r n in downscaling, which corresponds to energy flux density ε n = r n • 2 nc , preventing to mislead with the alive number at other steps.

When n = 1, there are three kinds of possible states at 1st cascade step:

All small eddies are alive (r 1 = 2) :

ε 1 = 2 • 2 c ; (2.3.26a)
One of the small eddies is alive (r 1 = 1) :

ε 1 = 1 • 2 c ; (2.3.26b)
N one of the small eddies is alive (r 1 = 0) :

ε 1 = 0 • 2 c . (2.3.26c)
The probability of state Eq.2.3.26b has to take arrangements of events into consideration too. So their probabilities are:

P (ε 1 = 2 • 2 c ) = ( 2 2 ) • P 2 1 ;
(2.3.27a)

P (ε 1 = 1 • 2 c ) = ( 2 1 ) • P 1 1 • (1 -P 1 ) 1 ; (2.3.27b) P (ε 1 = 0 • 2 c ) = ( 2 0 ) • (1 -P 1 ) 2 . (2.3.27c)
Plus, the sum of these events' probabilities is 1 ensuring probability conservation.

When n > 1, energy flux density at nth step that depends on states at former steps can be generated by different events, suggesting that probabilities should be obtained by conditional probability theory

P (A|B) = P (AB) P (B) . (2.3.28)
Event at nth step can be generated by numerous different states at (n -1)th step when only the state at (n -1)th step is taken into account. Therefore, events at the nth step can be generated by many possible paths by considering hierarchical former steps.

Alive number r n could be either the odd number or the even number. These situations will be separately discussed. r n = 2 n -2j (0 ≤ j < 2 n-1 ) alive eddies at nth step, which is even, should be generated from event r n-1 = 2 n-1 or event r n-1 = 2 n-1 -j if we only consider the (n-1)th step. Meanwhile, event r n-1 = 2 n-1 -j could be generated through different paths from 1st step to (n -2)th step which exhibits hierarchical property in the downscaling. Because dead eddies only generate dead eddies instead of alive eddies that has an impact on the event at nth step. As a result, 2j dead eddies at nth step, which are generated by different former layers, could consist of dead eddies produced by the last step or the former steps. In this way, there are plenty of arrangements of dead eddies which will be represented by coefficient a(m). h (z) dead eddies which are produced at (nz)th step except those eddies generated by former steps could be referred as binomial coefficient ( dead eddies at (nz)th step which is caused by former step. So h (z) • 2 z = 2j is the number of dead eddies at nth step, whereas the sum of h (z) is the index of the probabilities of dead eddy.

The following gives a clear formula of events probabilities:

P (ε n = (2 n -2j) • 2 nc ) =P (r n = 2 n -2j) =P (r n-1 = 2 n-1 ) • ( 2 n 2j ) • P 2 n -2j 1 • (1 -P 0 ) 2j + P (r n-1 = 2 n-1 -1)• ( 2 n -2 2j-2 ) • P 2 n -2j 1 • (1 -P 1 ) 2j-2 + • • • + P (r n-1 = 2 n-1 -j) • P 2 n -2j 1 = j m ′ =0 P (r n-1 =2 n-1 -m ′ ) • ( 2 n -2m ′ 2j-2m ′ ) • P 2 n -2j 1 • (1 -P 1 ) 2j-2m ′ = a(m) • P w 1 • (1 -P 1 ) v (2.3.29)
where a(m)

= n z=0 ( 2 n-z -g (z) h (z)
), in which

h (z) = v, h (z) • 2 z = 2j, g (z) = m ′ =z+1 h (m ′ ) • 2 m ′ -z , (2 n-z -g (z) -h (z) ) = w.
Another case is p n = 2 n -2j -1 alive eddies at nth step. The same procedure is easily adapted to obtain its probability.

P (ε n = (2 n -2j -1) • 2 nc ) =P (p n = 2 n -2j -1) =P (p n-1 = 2 n-1 ) • ( 2 n 2j+1 ) • P 2 n -2j-1 1 • (1 -P 1 ) 2j+1 + P (p n-1 = 2 n-1 -1) • ( 2 n -2 2j-2+1 ) • P 2 n -2j-1 1 • (1 -P 1 ) 2j-2+1 + • • • + P (p n-1 = 2 n-1 -j)( 2 n -2j 1 ) • P 2 n -2j-1 1 • (1 -P 1 ) = j t=0 P (p=2 n-1 -m ′ ) • ( 2 n -2m ′ 2j+1-2m ′ ) • P 2 n -2j-1 1 • (1 -P 1 ) 2j+1-2m ′ = a(m) • P w 1 • (1 -P 1 ) v (2.3.30)
), in which

h (z) = v, h (z) • 2 z = 2j + 1, g (z) = m ′ =z+1 h (m ′ ) • 2 m ′ -z , 2 n-z -g (z) -h (z) = w.
Moreover, the sum of probabilities of energy density at nth step equals 1 suggesting hierarchical theory exits in downscaling. Above two probabilities equations of downscaled β-model provide probability without computing probabilities of former events.

The detailed probability of up-scaled and down-scaled β-model are presented. It

shows that the measure of probability is conserved when scale changes.

α-model

As known, two states are used to roughly describe eddies in β-model. However, Schertzer and Lovejoy in 1984 discovered that various thresholds is allowed for a hierarchy of many fractals to characterize irregular field, which is the idea of multifractal.

Instead of that eddy is either "alive" or "dead", α-model is introduced to reveal multifractal field by considering intensity of activity that eddy is going to act more active or less active, as shown in where

γ + = c α ; γ -= - c α ′ ; 1 α + 1 α ′ = 1. (2.4.2) When γ -→ -∞, γ + = c which is β-model.
And there is only one unique singularity in β-model.

Due to Eq.2.3.4, 2.3.5, probabilities of increments are:

P r(µ = λ γ + ) = λ -c ; (2.4.3a) P r(µ = λ γ -) = 1 -λ -c . (2.4.3b)
As the averaged increments should be conversed ⟨µ⟩ = 1, α-model requires

⟨µ⟩ = λ γ + • P r(µ = λ γ + ) + λ γ -• P r(µ = λ γ -) = 1, (2.4.4)
where parameters are (γ + , γ -, c) and two of them could be chosen freely. The singularities in separating process indicate α-model is multifratal. Obviously, β-model is a special case that α = 1 and α ′ = 0 , which suggests β-model is a fractal model.

After iterating n = n + + n -step, probability of states at nth step is:

P r(ε n = (λ n ) γn ) = ( n n + )λ -cn + (1 -λ -c ) n -. (2.4.5)
where

γ n = n + γ + + n -γ - n + + n - ; (2.4.6)
n + is number of more active steps and n -is number of less active steps.

For those intensity states, probability is summed as the following:

P r(ε n ≥ (λ n ) γ i ) = j p ij (λ n ) -c ij (2.4.7)
where p ij is the prefactors in Eq 2.4.5, c ij is the sub-codimensions.

Multifractal field

By re-normalizing α-model, multifractal field is obtained. These n steps with two increment states cascade are replaced by a single step cascade with n + 1 states.

The new scale ratio λ ′ replaces λ n and λ ′ will eventually approach ∞. Here, c i is defined by

c i = min{c ij } = c(γ i ). (2.5.1)
And Eq 2.4.7 turns into For uniform and continuous distribution, Eq 2.5.2 could be written as

P r(ε λ ′ ≥ (λ ′ ) γ i ) = p i (λ ′ ) -c i . ( 2
P r(ε λ ≥ λ γ ) ≈ λ -c(γ) , (2.5.3)
where p i is absorbed by ≈. It describes probabilities of events whose intensity exceeds the scale-dependant threshold λ γ and yields a scale-invariant singularity γ, as shown in Fig 2 .10. Since Eq 2.5.2 indicates there is an unique co-dimension corresponding to scaling threshold, c(γ) in Eq 2.5.3 is an increasing function [START_REF] Schertzer | Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes[END_REF] of the singularity γ in Fig. Now let's move to scaling statistical moments of energy dissipation rate. The statistical moments of energy dissipation rate ⟨(ε λ ) q ⟩ has been given in Eq 2.2.5.

Schertzer and Lovejoy [START_REF] Lovejoy | The weather and climate: emergent laws and multifractal cascades[END_REF] obtained that the qth moment of energy dissipation rate ⟨(ε λ ) q ⟩ can be expressed with the resolution λ

⟨(ε λ ) q ⟩ ≈ λ K(q) (2.5.4)
where K(q) is moment scaling function which is convex. Moreover, Parisi and Frisch in 1985 [START_REF] Parisi | On the singularity spectrum of fully developed turbulence (1985) turbulence and predictability in geophysical fluid dynamics[END_REF] discovered the connection of scaling moment function K(q) and the codimension function c(γ) through Legendre transform (Appendix A.1):

c(γ) = max q (qγ -K(q)) = q γ γ -K(q γ ); (2.5.5a) K(q) = max γ (qγ -c(γ)) = qγ q -c(γ q ) (2.5.5b)
where q γ is defined by K ′ (q γ ) = γ and γ q is defined by c ′ (γ q ) = q. The statistical moment and singularity are unique corresponded. Because there are only convexity constraints, characteristic scaling function K(q) and codimension function c(γ) may depend on infinite number of parameters.

Universal Multifractal

The idea that only a few parameters are crucial in determining "universality" classes caught scientists' attention in the middle of 1980s. Based on the central limit theorem of additive process, Schertzer and Lovejoy in 1987 [START_REF] Schertzer | Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes[END_REF] obtained that multiplicative cascade processes converge to a "universality" process under quite general conditions, which is named as Universal Multifractal (UM). For conservative processes, these infinite parameters determined scaling moment function K(q) and codimension function c(γ), can be reduced to two.

When scaling moment is q = 1 indicating it is the mean of process, corresponding singularity is γ 1 . Due to Eq 2.5.5, it has

⟨ε λ ⟩ ≈ λ γ 1 -c(γ 1 ) . (2.6.1)
As the field is conserved ⟨ε λ ⟩ = 1, it has

γ 1 = c(γ 1 ) = C 1 (2.6.2)
Therefore, C 1 is the codimension of singularity which related to the mean field. Moreover, Legendre transforms yield

c(γ 1 ) = c(C 1 ) = C 1 , c ′ (γ 1 ) = c ′ (C 1 ) = 1 (2.6.3a) K(1) = γ 1 -c(γ 1 ) = 0, K ′ (1) = γ 1 = C 1 .
(2.6.3b)

Then multifractality can be characterized by the deviations of K(q) from the mean process q = 1. Since moment scaling function K(q) is convex, local radius of curvature R K at q = 1, which is local rate of change of slope, enables to define multifractality index α. The local radius of curvature R K at q = 1 is

R K (q = 1) = (1 + (K ′ (1)) 2 ) 3/2 K ′′ (1) = (1 + C 2 1 ) 3/2 K ′′ (1) (2.6.4)
And multifractality index α is defined as

α = K ′′ (1) K ′ (1) . (2.6.5)
Therefore C 1 is mean codimension measuring mean concentration of the activity.

Typically, C 1 of irregular field couldn't be smaller than 0 (C 1 ≥ 0). When C 1 equals 0, the process is homogeneous, and intermittency doesn't appear. Multifractality index α (0 ≤ α ≤ 2) measures how quickly intermittency evolves when singularity slightly deviates from the singularity obtained from average field. α = 0 is a monofractal field, such as β-model, while α = 2 is the Log-normal model [START_REF] Yaglom | The influence of fluctuations in energy dissipation on the shape of turbulence characteristics in the inertial interval[END_REF]. It's important to emphasize that α in UM is not the same as the one in α-model.

So the moment scaling function K(q) and codimension function c(γ) described by UM parameters [START_REF] Lovejoy | The weather and climate: emergent laws and multifractal cascades[END_REF] are (See in Appendix A.2)

K(q) =      C 1 α -1 (q α -q), α ̸ = 1 C 1 qLn(q), α = 1 (2.6.6a) c(γ) =        C 1 ( γ C 1 α ′ + 1 α ) α ′ , α ̸ = 1 C 1 exp( γ C 1 -1), α = 1 (2.6.6b) where 1 α + 1 α ′ = 1. (2.6.7)
Also included is multifractal impact on structural function S n (l) Eq 1.2.9 [START_REF] Schmitt | Multifractal analysis of the greenland ice-core project climate data[END_REF].

When a field is monofractal, scaling exponent ζ(q) is linear

ζ(q) = qH, (2.6.8)
however structure function of universal multifractal field exhibits deviations from linearity, which can be characterised by UM parameters

ζ(q) = qH - C 1 α -1 (q α -q), (2.6.9)
where H = ζ(1) is scaling Hurst exponent corresponding to the scaling of average absolute fluctuations. UM parameters could be obtained from plot of log 10 (qζ ′ (0)ζ(q)) and log 10 (q).

Double Trace Moment

UM parameters can be directly estimated by Double Trace Moment (DTM) method [START_REF] Lavallée | On the determination of universal multifractal parameters in turbulence[END_REF][START_REF] Veneziano | A modified double trace moment method of multifractal analysis[END_REF]. Considering renormalized ηth power of field whose resolution is λ ε

(η) λ = ε η λ ⟨ε η λ ⟩ , (2.6.10)
and statistical qth moment of ε

(η)
λ obtained from Eq 2.5.4 is

⟨ε (η)q λ ⟩ = ⟨ε ηq λ ⟩ ⟨ε η λ ⟩ q ≈ λ K(ηq) λ qK(η) = λ K(q,η) (2.6.11)
where K(q, η) = K(ηq) -qK(η).

(2.6.12)

In UM framework, moment function has

K(q, η) = η α K(q) (2.6.13)
and K(q, η) turns into

K(q, η) = η α C 1 α -1 (q α -q). (2.6.14)
It suggests the slope of log-log plot of K(q, η) and η is α and

C 1 = exp(b) α-1 q α -q
, where b is the intercept shown in Fig. 2.12. Figure 2.12: The function K(q, η) obtained by DTM method. UM parameters is estimated by linear fit of empirical points.

Sampling dimension

It is crucial to note that statistical properties of multifractal process in above are defined for infinite statistical ensembles. However, empirically statistical analysis normally is supported by finite samples instead of infinities that don't capture the whole rare events. Thus, sample dimension D s is introduced to describe finite samples by generating a maximum singularity γ s and critical moment q s , which is the maximum value that codimension and scaling moment function are reliable statistical estimated.

For a collection of N s samples, the resolution of samples is λ and there are λ D values for each D-dimensional sample. Based on the definition of codimension Eq 2.5.3, there is a maximum value of γ, which is denoted γ s , and could be observed on at least one sample. Hence, it yields:

N s λ D λ c(γs) ≈ 1. (2.6.15)
And the sampling definition is defined as

N s = λ Ds ⇒ D s = log(N s ) log(λ) .
(2.6.16)

From Eq 2.6.15, it has

c(γ s ) = D + D s .
(2.6.17)

According to Fig 2 .11, when the sampling dimension is considered, the maximum singularity gets large indicating that more rare events are studied.

The critical moment q s = c ′ (γ s ) which depends on number of samples, is always used to study scaling moment function. It is the highest order moment after which estimated moment scaling function K(q) will become linear. The critical moment q s = c ′ (γ s ) obtained from the UM framework is:

q s = ( D + D s C 1 ) 1 α .
(2.6.18)

Summary of Chapter 2

In this chapter, the most basic phenomenological cascade β-model considering the intermittency is explained. Since many studies have attempted to understand the intermittency with its' probability distribution and the physical β-model gives the probabilistic representation, the probability conservation under the changing scales is studied. It demonstrates that probability conservation is a scaling invariant, which lays a solid foundation for the further developed multifractal field. However, the alive eddies at high moments will take more weights, causing the probability field to be deformation of NS equation, which is particularly difficult to obtain the general solution. These methods are the simplest explicit Euler method, the classic fourth-order Runge-Kutta method, and two methods considering physical properties. They are the semi-implicit Euler method, which is often used to describe the motion of rigid bodies, and the slaved Adams-Bashforth method, which takes into account the fast damping of waves with high wave numbers. Algorithm analysis is performed, and the results obtained from these three approaches are compared in order to identify effective numerical simulation approach.

Scaling Gyroscope Cascade (SGC) model

Scaling Gyroscopes Cascade (SGC) model [START_REF] Chigirinskaya | Cascade of scaling gyroscopes: Lie structure, universal multifractals and self-organized criticality in turbulence[END_REF][START_REF] Chigirinskaya | Scaling gyroscopes cascade: universal multifractal features of 2-d and 3-d turbulence[END_REF][START_REF] Schertzer | Multifractal cascade dynamics and turbulent intermittency[END_REF] Bernoulli form of NS equation is

( ∂ ∂t -v∆)u(x, t) = u(x, t) ∧ ω(x, t) -grad(α * ); (3.1.2a) ω(x, t) = curl(u(x, t)) (3.1.2b)
where ω(x, t) is vorticity and α * is kinematic pressure. Projector P (∇) corresponding to P (k) in Fourier space

P i,j (∇) = δ i,j -∇ i ∇ j ∆ -1 (δ i,j is Kronecker ′ s δ) (3.1.3)
imposes the imcompressibility condition Eq.3.1.1 restriction on advection term, classic pseudospectral technique typically employed in DNS [START_REF] Pope | Turbulent flows[END_REF][START_REF] Canuto | Spectral methods in fluid dynamics[END_REF], and Eq.3.1.2 in Fourier space turns into:

( ∂ ∂t + vk 2 ) u(k, t) = P (k) • p+q=k u(p, t) ∧ ω(q, t)d d p, (3.1.4)
where ω is vorticity:

ω i m = ik i m ∧ u i m . (3.1.5)
Here, the velocity-vorticity vertex of triad interaction (k, p, q) also only exits under orthogonality condition

k • u(k, t) = 0. (3.1.6)
Assuming arbitrary nonlocalness parameter λ k in Eq.1.6.16, the orthogonality of vorticity-velocity vertex interaction is

|k| ≪ |p| ≈ |q|, p ⊥ k =⇒ ( u(p) ∧ ω(q)) ⊥ k (3.1.7a) |p| ≪ |k| ≈ |q|, u(p) ∥ k =⇒ ( u(p) ∧ ω(q)) ⊥ k, ( u(q) ∧ ω(p)) ∥ k. (3.1.7b)
Therefore, Eq.3.1.4 turns into:

( ∂ ∂t + vk 2 ) u(k, t) = |p|>λ k |k| ( u(p) ∧ ω(q))d d p + ( |p|≤λ -1 k |k| u(p)d d p) ∧ ω(q). (3.1.8)
Since the triad interaction of SGC model is chosen as its tree structure (k i m , k 2i m+1 , k 2i+1 m+1 ), which has orthogonality constraints, Eq. 3.1.8 becomes

( ∂ ∂t + vk 2 ) u i m = u 2i m+1 ∧ ω 2i m+1 + u 2i+1 m+1 ∧ ω 2i+1 m+1 + u a(i) m-1 ∧ ω i m . (3.1.9)
The symmetric property of gyroscope equation yields

u(k, t) = ik ∧ ω(k, t) k 2 . ( 3
k 2i+1 m+1 = -k 2i m+1 , (3.1.11)
which is the last orthogonality triad.

Following the matrix representation of Eq.3.1.9 considering orthogonality, discrete SGC model is as follow:

( ∂ ∂t + vk 2 m )u i m = k m+1 [|u 2i m+1 | 2 -|u 2i+1 m+1 | 2 ] + (-1) i+1 k m u i m u a(i) m-1 (3.1.12)
where k m is wave number at layer m (0 ≤ m ≤ n); k m+1 is wave number at layer

m + 1; u i m is velocity of eddy in location i at layer m (0 ≤ i ≤ 2 m -1); a(i)
(t) ∂t = k 1 [|u 0 1 (t)| 2 -|u 1 1 (t)| 2 ] -vk 2 0 u 0 0 (t). (3.1.13)
Velocities at last layer u i n without generating daughter eddies are

∂u i n (t) ∂t = (-1) i+1 k n u i n (t)u a(i) n-1 (t) -vk 2 n u i n (t). (3.1.14)
According to the invariant of gyroscope equation, turbulence energy corresponds to the square of the angular momentum M 2 conserves.

Waleffe [START_REF] Waleffe | The nature of triad interactions in homogeneous turbulence[END_REF][START_REF] Waleffe | The helical decomposition and the instability assumption[END_REF] reported that the instability assumption for rapidly rotating turbulence, which is typically applicable for any sort of triad interaction, suggests direction of energy transferring -from the poles of the rotation axis to the equator.

Furthermore, this instability assumption is made even more apparent with helical mode decomposition. As an alternative to the Craya-Herring decomposition, helical mode decomposition [START_REF] Cambon | Spectral approach to non-isotropic turbulence subjected to rotation[END_REF][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF] allows for departures from isotropy caused by external rotational forces and the spectral energy to move from parallel to the normal wave vector, which seems to be compatible with trend toward two-dimensionality. Furthermore, Benzi et al. [START_REF] Benzi | Helical shell models for threedimensional turbulence[END_REF] investigated in detail various classes of helical shell models that arise from helical decomposition of interactions, and four different types were discussed. Two of them match the two-dimensional and three-dimensional shell models.

Consequently, this helical decomposition might be applied to SGC model to provide more specific details.

The numerical methods

Four different numerical algorithms are discussed in this section to determine the scheme contributing to time developed velocity field of SGC model. Euler method, the simplest numerical approach, is presented first, followed by widely used, high precision fourth-order Runge-Kutta method. As a next step, slaved Adams-Bashforth method considering high dissipation at large wave numbers and the semi-implicit Euler scheme that maintains the rigid body property are studied. A brief overview of these methods will be provided, and the corresponding Python code and Pseudo-code will be provided in Appendix B.3.

The Euler method

Euler method [START_REF] Biswas | A discussion on euler method: A review[END_REF], first reported by Leonhard Euler in 1768, utilizes the elementary difference methods to obtain approximate solutions to differential equations or initial value problems. Theory underlying Euler's method of using the concept of local linearity of small line segments to approximate solutions is well illustrated by Taylor's theorem, which was published in 1714. Algorithm for Euler method is investigated since it is simple and its discretization of continuous variables can be easily calculated numerically.

Typical form of the first-order differential equation is

dy dt = f (t, y(t)) (3.2.1)
with initial condition y(t 0 ) = y 0 . The classic first-order forward Euler method discretizes Eq.3.2.1 writing as

y(t + ∆t) -y(t) ∆t ∼ = f (t, y(t)) y(t + ∆t) ∼ = y(t) + ∆t • f (t, y(t)), (3.2.2)
where ∆t is the time step.

Numerically simulated velocity equations using Euler's method are categorized into three functions due to the spatial structure of SGC model. First of all, velocity at the top layer u 0 0 is computed by

u 0 0 (t + ∆t) = u 0 0 (t) + ∆t • k 1 [|u 0 1 (t)| 2 -|u 1 1 (t)| 2 ] -∆t • vk 2 0 u 0 0 (t). (3.2.3)
External forcing is added to the code as forced turbulence is going to be studied in the next chapter, which is not necessary to repeat.

Velocities u i n located at the last layer are

u i n (t + ∆t) = u i n (t) • (1 -∆t • vk 2 n ) + ∆t • (-1) i+1 k n u i n (t)u a(j) n-1 (t). (3.2.4)
As for velocities u i m whose cascade step m is between 1 and n -1, they are simulated by

u i m (t + ∆t) =u i m (t) • [1 -∆t • vk 2 m + ∆t • (-1) i+1 k m u a(i) m-1 (t)] + ∆t • k m+1 [|u 2i m+1 (t)| 2 - |u 2i+1 m+1 (t)| 2 ]. (3.2.5)
Furthermore, the numerical convergence of Euler's method has been demonstrated by Augustine Louis Cauchy in 1824, thus it's clear that time steps must be sufficiently small to ensure convergence. The absolute stability region of Euler method is

|1 + f (t, y(t)) y • ∆t| < 1, (3.2.6)

Fourth-order Runge-Kutta method

Runge-Kutta method is a conventional, one-step approach that has been utilized extensively due to its high accuracy. It was initially proposed by Carl Runge in 1895 [START_REF] Runge | Über die numerische auflösung von differentialgleichungen[END_REF] and later extended to solve systems of differential equations by Martin Wilhelm Kutta in 1901 [START_REF] Kutta | Beitrag zur naherungsweisen integration totaler differentialgleichungen[END_REF]. To avoid computing higher order derivatives in Taylor expansion technique while still achieving higher order convergence, Runge-Kutta approach inserts a series of points into range [t, t + ∆t] to anticipate the slopes k 1 , k 2 • • • , and then takes their weighted average. Obviously, Euler method is Runge-Kutta method of the first order. The discretized Eq.3.2.1 for fourth-order Runge-Kutta (RK4) is

y(t + ∆t) = y(t) + ∆t 6 (k * 1 + 2k * 2 + 2k * 3 + k * 4 ) (3.2.7)
where

k * 1 = f (t, y(t)) k * 2 = f (t + ∆t 2 , y(t) + ∆t 2 • k * 1 ) k * 3 = f (t + ∆t 2 , y(t) + ∆t 2 • k * 2 ) k * 4 = f (t + ∆t, y(t) + ∆t • k * 3 ). (3.2.8)
Except that RK4 method has to calculate the slope four times, numerical velocities equation is the same between the Euler and RK4 methods, therefore, Python code for RK4 method is directly provided. The equations describing velocity field of all layers for each time step are written in pseudo-code and Python code, which are quite similar to that of Euler and in Appendix B.3.2.

Besides, the absolute stability region of RK4 method is

|1 + f (t, y(t)) y • dt + 1 2! • ( f (t, y(t)) y • dt) 2 + 1 3! • ( f (t, y(t)) y • dt) 3 + 1 4! • ( f (t, y(t)) y • dt) 4 | < 1.
(3.2.9)

Slaved Adams-Bashforth Method

The linear multi-step method is another approach that is typically used to resolve ordinary differential equations. Unlike the one-step method, which simply discards former data in favor of the most recently acquired data, this method uses a linear combination of former points and derivative values to predict the next point. One of the classic multi-step approaches is Adams-Bashforth method, commonly known as explicit Adam method. In order to account for the quick damping at high wave numbers, "slaved frog" scheme [START_REF] Frisch | Viscoelastic behaviour of cellular solutions to the kuramoto-sivashinsky model[END_REF] has been implemented in two-step Adams-Bashforth method, which is slaved Adams-Bashforth (ABF) approach and has been applied to the investigation of shell model [START_REF] Pisarenko | Further results on multifractality in shell models[END_REF].

A simple linear model that displays fast damping at high wave numbers is

dU dt = -aU + g(t). (3.2.10)
After integrating, it turns into

U (t + δt) = e -aδt U (t) + t+δt e -a(t+δt-s) g(s)ds U (t + δt) = e -aδt U (t) + 1 -e -aδt a g(t) (3.2.11)
Hence, the slaved Adams-Bashforth Method is:

U (t + ∆t) = e -a•∆t u(t) + 1 -e -a•∆t a ( 3 2 g(t) - 1 2 g(t -∆t)). (3.2.12)
The initial time t 0 is estimated by g(t 0 -∆t) = g(t 0 ).

Velocity equation u 0 0 obtained from slaved ABF method is

u 0 0 (t + ∆t) = e -a 0 •∆t u 0 0 (t) + 1 -e -a 0 •∆t a 0 ( 3 2 g 0 (t) - 1 2 g 0 (t -∆t)), (3.2.13)
where a 0 = vk 2 0 and g 0

(t) = k 1 [|u 0 1 (t)| 2 -|u 1 1 (t)| 2 ].
Velocities at the end layer u i n are computed by

u i n (t + ∆t) = e -an•∆t u i n (t) (3.2.14)
where

a n = vk 2 n -(-1) i+1 k n u a(j) n-1 (t).
Velocities

u i m (1 ≤ m ≤ n -1) are computed from u i m (t + ∆t) = e -am•dt u i m (t) + 1 -e -am•∆t a m ( 3 2 g m (t) - 1 2 g m (t -∆t)) (3.2.15)
where

a m = vk 2 m -(-1) i+1 k m u a(j) m-1 (t) and g m (t) = k m+1 [|u 2i m+1 (t)| 2 -|u 2i+1 m+1 (t)| 2 ].
Despite the fact that two schemes are applied to increase the precision of turbulence numerical simulations, slaved ABF technique still has to adhere to the stability of ABF approach. The characteristic equation of ABF method satisfies root condition, suggesting that it is highly stable and avoids the round-off error from growing exponentially.

Semi-implicit Euler method

Semi-implicit Euler method known as Störmer-Verlet and symplectic-Euler [START_REF] Hairer | Structure-preserving algorithms for ordinary differential equations[END_REF] is commonly applied to the numerical simulations of rigid bodies, as it's more stable than classic explicit Euler schemes and preserves the nature of volume conservation of phase space volume element in physical system. Hence, semi-implicit Euler scheme is studied in the simulation of SGC model, which is obtained from the gyroscope equation-Euler's equation for a rigid body. Semi-implicit Euler scheme for SGC model, in which linear dissipation is implicit and nonlinear transfer term is explicit, is as follows:

(

1 ∆t + vk 2 m )u i m (t + ∆t) - u i m (t) ∆t = g i m (t) u i m (t + ∆t) = g i m (t)∆t + u i m (t) 1 + vk 2 m ∆t (3.2.16)
where the nonlinear transfer term is wrote as

g i m (t) = k m+1 [|u 2i m+1 | 2 -|u 2i+1 m+1 | 2 ] + (-1) i+1 k m u i m u a(i) m-1 .
Velocity u 0 0 (t + ∆t) at the top layer is computed by:

u 0 0 (t + ∆t) = k 1 [|u 0 1 (t)| 2 -|u 1 1 (t)| 2 ] • ∆t + u 0 0 (t) 1 + vk 2 0 • ∆t (3.2.17)
Velocities u i n located at last layer are:

u i m (t + ∆t) = (-1) i+1 k n u i n (t)u a(i) n-1 (t) • ∆t + u i m (t) 1 + vk 2 m • ∆t (3.2.18)
For velocities u i m whose cascade step m is between 1 and n -1, they are simulated by

u i m (t + ∆t) = (k m+1 [|u 2i m+1 (t)| 2 -|u 2i+1 m+1 (t)| 2 ] + (-1) i+1 k n u i n (t)u a(i) n-1 (t)) • ∆t + u i m (t) 1 + vk 2 m ∆t (3.2.19)

Result analysis

In order to determine an optimal numerical method for SGC model, it is essential to compare velocity fields generated by different approaches, such as energy flux utilized to study intermittency and it's multifractality analyzed by UM framework.

Moreover, complexity of the principal function as well as the specific running time and space of these algorithms will be calculated in next section to evaluate numerical methods.

To begin, the specific equation for energy flux and required simulation parameters

Π(k m ) = - m-1 m ′ =0 T (k m ′ , t) = n m ′ =m T (k m ′ , t) ∝ n m ′ =m u(k m ′ , t) • [ ∂u(k m ′ , t) ∂t + vk 2 m ′ u(k m ′ , t)]. (3.3.3)
The energy flux through wave number k m of SGC model is

Π(k m ) ∝ n j=m 2 m -1 i=0 u i j (t) • [k j+1 [|u 2i j+1 (t)| 2 -|u 2i+1 j+1 (t)| 2 ] + (-1) i+1 k j u i j (t)u a(i) j-1 (t)] (3.3.4)
which is related to the nonlinear transfer term of SGC model.

The initial kinetic energy spectrum is assumed as Kolmogorov spectrum

E(k, 0) = C K ε 2 3 k -5 3 , (3.3.5)
where C k is Kolmogorov constant. C k and ε are taken to be 1 to start simulation.

Initial velocities u i m (0) are obtained from the assumed energy spectrum E(k m , 0) by

|u i m (0)| 2 = E(k i m , 0)k i m ≡ E(k m , 0)k m . (3.3.6)
An artificial viscosity v is required to dissipate the energy and it can be derived from

Reynolds number at Kolmogorov scale Re(l η * ) ≈ 1 For simplicity's sake, the external forcing f is assumed to be 0, which eliminates the complexity of turbulent characteristics brought on by injecting forcing.

Re(l η * ) = u η * l η * v = u η * vk η * ≈ 1. ( 3
Lines 3 2 n+1 -1 O(2 n ) Function 1 in Algorithm 1 2 n+2 -1 O(2 n ) Function 1 in Algorithm 2 3 • 2 n+1 -2 O(2 n )
However, storage space S(n) of slaved ABF method is larger than that of one-step methods, demonstrating that slaved ABF method requires more usage memory for each time step. And the difference in S(n) rapidly grows as cascade step n enlarges as illustrated in Fig. 3.5.

T (n) g(n) = C ̸ = 0, (3.4.2) then T (n) is denoted as T (n) = O(g(n)) and O(g(n)) is asymptotic time complexity.
Below is Table 3 Numerous achievements in turbulence [START_REF] Franzke | The structure of climate variability across scales[END_REF] are challenged by intermittency, which commonly refers to spatially and temporally unevenly distributed velocity fluctuations. Because of "the heavy tails" generated by extreme events, it is well acknowledged that intermittency is non-Gaussian distribution. Advance studies [START_REF] Aitchison | the lognormal distribution[END_REF], such as K62, claim that it follows log-normal distribution, which has been proven in hydrodynamic turbulence. However, log-normal model is criticized by universal behavior [START_REF] Schertzer | Nonlinear variability in geophysics[END_REF][START_REF] Gupta | A statistical analysis of mesoscale rainfall as a random cascade[END_REF], UM framework based on the central limit theorem is then considered a powerful approach to characterize this intermittency. Empirical UM parameters estimated from various turbulence data, especially from laboratory experiments [START_REF] Schmitt | Empirical determination of universal multifractal exponents in turbulent velocity fields[END_REF] and atmospheric in-situ/remotely sensed data [START_REF] Schmitt | Estimation of universal multifractal indices for atmospheric turbulent velocity fields[END_REF], are multifractality index α ≈ 1.5 and mean co-dimension C 1 ≈ 0.25. Moreover, results of SGC model simulated by RK4 method [START_REF] Chigirinskaya | Cascade of scaling gyroscopes: Lie structure, universal multifractals and self-organized criticality in turbulence[END_REF][START_REF] Schertzer | Multifractal cascade dynamics and turbulent intermittency[END_REF][START_REF] Chigirinskaya | Scaling gyroscopes cascade: universal multifractal features of 2-d and 3-d turbulence[END_REF] confirm that forced SGC model whose cascade step is 12 exhibits strong fluctuations in energy transfer process and UM parameters are consistent with those obtained from empirical data. Therefore, the existence of intermittency and its multifractality is firstly investigated by forced SGC model whose cascade step is 12

and then expanded to large cascade steps 14, 15 in this chapter.

The forced SGC model

Forced SGC model instead of freely decaying SGC model is simulated to achieve quasi-equilibrium state, which prolongs turbulence motion and hence exhibits more information, such as strong fluctuations, for studying the intermittency. A quasiequilibrium state [START_REF] Tatsumi | Theory of homogeneous turbulence[END_REF] is a statistically stable condition in which viscous dissipation is not precisely balanced by external energy and system kinetic energy does not remain constant during the whole process, as it does in the equilibrium state. Despite the fact that the system is not self-preserving, it is nevertheless possible to obtain a similar law partly.

There are two fundamental approaches of forcing [START_REF] Sullivan | Deterministic forcing of homogeneous, isotropic turbulence[END_REF]: deterministic forcing and stochastic forcing [START_REF] Alvelius | Random forcing of three-dimensional homogeneous turbulence[END_REF][START_REF] Eswaran | An examination of forcing in direct numerical simulations of turbulence[END_REF]. Evidence [START_REF] Chigirinskaya | Cascade of scaling gyroscopes: Lie structure, universal multifractals and self-organized criticality in turbulence[END_REF] suggests, however, that SGC model isn't largely dependent on the type of external forcing. Specifically, simple deterministic forcing f is utilized in the study. Forced SGC model is

∂u i m ∂t = k m+1 [|u 2i m+1 | 2 -|u 2i+1 m+1 | 2 ] + (-1) i+1 k m u i m u a(j) m-1 -vk 2 m u i m + f. (4.1.1)
If external force is significantly larger than viscous term (f ≫ vk 2 m u i m ), Eq 4.1.1 turns

into ∂u i m ∂t = k m+1 [|u 2i m+1 | 2 -|u 2i+1 m+1 | 2 ] + (-1) i+1 k m u i m u a(j) m-1 + f (4.1.2)
which gives unstable state or even worse explodes the model. Conversely, when forcing is less than viscous term, Eq 4.1.1 is equivalent to freely decaying SGC model. To sum up, an appropriate forcing f is essential to achieve quasi-equilibrium state, and deterministic forcing is performed by merely adding forcing at the top layer. Velocities as well as artificial viscosity are given according to Eq 3.3.6 and Eq 3.3.7 in Chapter 3.

Numerical parameters for the forced SGC model

Typical time unit of turbulence is initial large-eddy turnover time [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF] that allows significant distortion of eddy to transfer energy:

τ e = L(0) (E(0)) 1/2 , (4.2.1)
where integral scale [5] measures the correlation distance between two points in terms of distance or time characterizing energy containing scale:

L(t) = kmax 0 k -1 E(k, t)dk kmax 0 E(k, t)dk (4.2.2)
and system energy is computed by

E(t) = m=n m=0 2 m -1 i=0 |u i m (t)| 2 . (4.2.3)
Due to specific spatial structure of SGC model, time evolution of energy spectrum

E(k m , t) is calculated by E(k m , t) = 1 2 m 2 m -1 i=0 |u i m (t)| 2 ≈ E m (t), (4.2.4)
where

E m (t) = 2 m -1 i=0 |u i m (t)| 2 (4.2.5)
is energy at layer m whose wave number is k m . The following section will give out temporal evolution of energy fluxes injected with various forcings. Energy flux confirms the existence of intermittency and its probability distribution is provided to testify that it's non-Gaussian distribution and has a "heavy" tail. Then UM analysis is performed to describe the multifractality.

Cascade step 12

Intermittency of energy flux

Non-Gaussian distribution and "heavy" tail

Probability distribution is given to confirm that fluctuations in turbulence are non-Gaussian distribution and have a "heavy" tail [START_REF] Franzke | The structure of climate variability across scales[END_REF].

Probability distribution of fluctuations ∆|Π(k 6 )| for Case 4.1 exceeding threshold s is presented according to probability definition of multiplicative cascade Eq.2.5.3. Temporal evolution of UM parameters and critical moment q s obtained from small sub-sequence size clearly reveals more detailed information, such as the trend, since results of size 512 in Fig. 4.5b and size 256 in Fig. 4.5c provide more elaborated descriptions than the analysis of size 1024 in Fig. 4.5a. The trends of UM parameters and q s for size 256 are largely similar to those for size 512 with the exception of α in the final period, which results in changes of q s at the same time. UM parameters obtained by ensemble analysis are almost the same as shown by frequency histograms.

As seen in

Multifractality analyzed by UM framework

α varies as a result of changes in sample sizes, so does q s derived by ensemble analysis, which would not be considered. Frequency histogram of α for size 512 is more comparable that for size 256 rather than size 1024. Almost identical C 1 are found over all three sizes, and its frequency histogram shows nearly predicted double growth.

In summary, size 512 is chosen for UM analysis of energy flux in the following case study.

It should be emphasized that UM parameters1 in above is derived from a prede- The bias results from the fact that convergence of certain computational quantities is determined by field distribution [START_REF] Desaulniers-Soucy | Empirical test of the multifractal continuum limit in rain[END_REF]. Sometimes only a small fraction of the field primarily contributes to certain statistical moments. The convergence of trace moments may only occur at rather large scales and varies with sub-series, indicating that "flat" region at the upper right for each sub-series is different. For this reason, fit correlation coefficient for the plot is limited to a value larger than 0.995 in order to locate optimal range closest to 1 for each sub-series to ensure that slope is well fitted. then floats steadily after N = 30. The tendency of q s thus is decreasing, followed by end period of steadily fluctuating steadily.

α is between 1 and the Levy index in Log-normal model 2, while C 1 ranges between 0 and 0.5. Ensemble analysis considering number of samples yields larger values q se than averaged q s . As for UM parameters, it's typical for discrepancy between ensemble analysis and averaged values. Moreover, time evolution of UM parameters reflects time period close to atmospheric turbulence whose empirical parameters are α ≈ 1.5, C 1 ≈ 0.25. The number of subsequences centered around α ≈ 1.5 and C 1 ≈ 0.25 is visible on frequency plots. q s can also be utilized to study atmospheric turbulence whose q s is around 2.5. imum is not far from 1. C 1 grows over time, but it's less than 1. Both q s suggesting by q s and q se because of large sample size are far away from that of the atmospheric turbulence. In the beginning, mean codimesnion C 1 is around 0.1 as injected forcing

Analysis of

has not yet taken effect and the system is reaching quasi-equilibrium. UM analysis and frequency plots reveal the overlap period that α near 1.5 and C 1 close to 0.25 and subsequences number for ensemble analysis. 2 The first set of results is time evolution of UM parameters and q s obtained from fixed range method, while the second set of results is obtained from optimal range method. For each set, the first column displays time evolution of multifractality index α and its frequency, and mean codimension C 1 is the second column. Critical moment q s is in the last column. The outcomes for remaining cases are organized in the same structure and won't be repeated. Even though two methods yield nearly identical trends for C 1 that float steadily in intermediate phases, mean values and frequency histograms of UM parameters as well as q s testify the difference caused by estimation range. Ensemble analysis yields identical UM parameters and q s for two different range estimations. 11b. C 1 is small at the beginning and has an increment at the end. q s is larger than 2.5, and is close to q se as the number of sample is reduced.

Time evolution of UM parameters indicates that the intermediate sub-series are close

to that of atmospheric turbulence, while frequency plots give the number of specific sub-series. turbulence, making it difficult to conduct ensemble analysis. This analysis concludes that forcing 0.9 is excessive not only stopping SGC model but also generating subseries far away from atmospheric turbulence. Scaling moment equation K(q) is indistinguishable from theoretical scaling moment equation until critical moment q s = 2.67, after which empirical curve becomes linear.

Cascade step 14

Intermittency of forced SGC model with a large cascade step n = 14 is investigated, and the range of external forcings is from 0. The probability distribution compared with corresponding Gaussian distribution displays "heavy" tails for Case 4.2. Hence, UM framework is applied to study this extreme events causing "heavy" tails.

Multifractality described by UM framework

UM parameters for Case 4. obtained from ensemble analysis are close to mean values, while q s obtained from ensemble analysis are much larger than mean values, which is consistent with the definition of q s taking into account sample dimension.

Non-Gaussian distribution and "heavy" tail

Probability distribution of fluctuations for Case 4.3 is given in this part. The definition of q s implies that it is normal for Case 3.6 that q s yielded from two methods are not the same. The UM parameters obtained by the optimal range are more reasonable and accurate, as shown by the fact that the mean value of the UM parameters in Fig. 1.5b is higher than that in Fig. 1.5a. In addition,

9

(b) Adjusted range estimation Two range estimation methods yield nearly identical trend curves for UM parameters and q s . Frequency plots suggest that the means haven't changed significantly, but a few modifications have been made to C 1 . Changes in α lead to adjustments of the frequency plot and the mean of q s q s obtained from ensemble analysis is larger than the mean due to sample number. UM parameters in the early phase have one value α greater than 2 and several C 1 lower than 0.1 as the model is trying to reach quasi-equilibrium before injected forcing is not yet acting. When forcing is large, the time scale for analysis cannot reach expected time, preventing it from providing more information. Considering that C 1 has only two subseries around 0.25, there are no more than three subseries approaching atmospheric turbulence for ensemble analysis. The number of subseries with α ≈ 1.5 and C 1 ≈ 0.25, as shown by frequency histograms, is summarized in Table 4 Chapter 5

Backscatter term and energy decay law 5.1 Global scaling and the decay law of free turbulence Since at least the empirical observations made by Leonardo da Vinci, turbulence has been known to be dissipative, i.e. without external forcing it will decay. At the mathematical level, this is due to the presence of dissipative term ν∆u in the Navier-Stokes equation. However, it has also been observed that this decay is slow, in fact much slower than an exponential law that would result from neglecting all the other terms of these equations. More precisely it was empirically found and theoretically argued that the decay of free turbulence follows a power law for the whole period a non-negligible viscosity. On the contrary, the often studied case s = 4, see below, requires a significantly modified dissipativity α D (s) = 1.75 to be in the situation of strict global self-similarity.

Stationarity of big eddies and energy backscattering

The stationarity of big eddies was largely discussed by Kolmogorov [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF] and Landau and Lifchitz [START_REF] Landau | Lehrbuch der theoretischen physik[END_REF] in relation with the so-called Loitsianski's invariant or Loitsianski's integral Eq.1.5.5:

I =< u 2 (t) > ∞ 0 r 4 f (r)dr (5.2.1)
where f (r) is the longitudinal correlation. Some physical arguments [START_REF] Davidson | Was loitsyansky correct? a review of the arguments[END_REF] have been put forward, notably by [69] to ensure that this integral is indeed time-invariant, however, no mathematical proof has yet been provided and there is reason to believe that it will be not done. Indeed, a basic feature of Loitsianski's integral is that it is proportional to the prefactor C 4 (t) of k 4 term of energy spectrum expansion at k → 0, while the analytical closures, starting with Quasi-Normal approximation [START_REF] Millionshchikov | On the theory of homogeneous isotropic turbulence[END_REF] point out [START_REF] Proudman | On the decay of a normally distributed and homogenous turbulent velocity field[END_REF] that:

T (k, t) = A(t)k 4 + O(k s+1 ), k → 0; s > -1; A > 0 (5.2.2)
The fact that A(t) is positive implies that the largest eddies are growing and storing energy and thus have the potential to slow down turbulence decay. This is very effective if initially s ≥ 4, particularly if s ≫ 4, modifying the scaling exponent expressions (Eq.5.1.4). On the contrary, it will have a negligible impact for s < 4 and the scaling exponent expressions (Eq.5.1.4) will apply unchanged. The nonrealisability of Quasi-Normal approximation at long times does not bring into question this behavior for short times. This is more than being confirmed by Lesieur and

Schertzer [START_REF] Lesieur | Amortissement autosimilaire d'une turbulence à grand nombre de reynolds[END_REF] with the help of Eddy Damped Quasi-Nomal Markovianized model (EDNM) [START_REF] Orszag | Statistical theory of turbulence in fluid dynamics 1973 les houches summer school of theoretical physics[END_REF]. They demonstrated that for all times:

A(t) = 7 15 ∞ k/a θ k,p,p E(p, t) 2 p 2 dp, (5.2.3)
where θ k,p,q is the correlation relaxation time of the triad k, p, q. The main contribution to A(t) comes from energy-containing eddies of wave numbers p ∼ ℓ -1 and therefore does not depend on either of wavenumber k or the non-locality parameter a (k ≤ ap, q ∼ p). This term was initially called 'beating term' because it corresponds to interactions between two wavenumbers very close to each other p ∼ q. But, due to the fact that it scatters energy from energy-containing eddies back to the largest eddies, it became known as 'backscatter term'.

For initial s ≥ 4, Lesieur and Schertzer [START_REF] Lesieur | Amortissement autosimilaire d'une turbulence à grand nombre de reynolds[END_REF] introduced the scaling exponent γ > 0 of the Loitsianski's integral and the corresponding prefactor C 4 (t):

C 4 (t) = C 4 (t 0 )(t/t 0 ) γ , (5.2.4)
which modifies the scaling exponents of v and ℓ (Eq.5.1.4) into:

α v = - 5 -γ 7 ; α ℓ = 2 + γ 7 , (5.2.5)
with the following numerical estimates:

α e = -2α v = 1.
38; α ℓ = 0.31; γ ≈ 0.17.

(5.2.6)

Lesieur and Schertzer [START_REF] Lesieur | Amortissement autosimilaire d'une turbulence à grand nombre de reynolds[END_REF] also argued that the fact that backscatter is in k 4 is also obtained by the renormalization group [START_REF] Forster | Large-distance and long-time properties of a randomly stirred fluid[END_REF]. It is also worthwhile to note that most of the aforementioned results on the velocity field were extended by Herring et al., [START_REF] Herring | A comparative assessment of spectral closures as applied to passive scalar diffusion[END_REF] to the passive scalar field, not only with the help of EDQNM, but also with the Test Field Model (TFM), which can be understood as somewhat more elaborated. Frisch et al [START_REF] Frisch | Comments on the quasi-normal markovian approximation for fully-developed turbulence[END_REF] showed also these results could be also obtained with the less markovianised version of the Quasi-Normal model (QNM) [START_REF] Tatsumi | The multiple-scale cumulant expansion for isotropic turbulence[END_REF], which is fully realizable contrary to the original QN, but still fails to renormalize the Green functions, contrary to the Test-Field-Model(TFM).

Other studies of backscatter term

More detailed research achievements are given in the following since global scaling law and the role of backscatter term in the stationarity of big eddies have been explained in the previous two sections.

The permanence of large eddies has been extensively investigated using the classic two-point closure theories, and thus indirectly confirms the existence of energy backscatter term. Except for closure analysis methods, the existence of energy backscatter is studied by the commonly numerical approach for turbulence. Using large eddy simulation (LES) [START_REF] Lesieur | New trends in large-eddy simulations of turbulence[END_REF], Ossia and Lesieur [START_REF] Ossia | Energy backscatter in large-eddy simulations of threedimensional incompressible isotropic turbulence[END_REF] examined energy backscattering in free decay of three-dimensional incompressible isotropic turbulence and confirmed the non-negligible energy backscattering term Eq.5.2.2 for initial power spectrum s = 4. The exponent of energy decay law is -1.4 which is between -1.38

proposed by EDQNM and -1.43 predicted by Kolmogorov. However, the quantitative measure of energy backscatter intensity γ relies on the resolution and integration time. Chasnov [START_REF] Chasnov | Computation of the loitsianski integral in decaying isotropic turbulence[END_REF] observed the continuous energy backscattering term and obtained the estimated exponent γ ≈ 0.25 from spectral velocity field, when the finite size effects are negligible, through computing 1024 independent cases and large resolution 643 to improve statistical convergence of the infrared modes.

Another interesting argument is the threshold for destabilizing stationarity of large eddies, which is commonly acknowledged as an integer 4. With EDQNM framework, Eyink and Thomson [START_REF] Eyink | Free decay of turbulence and breakdown of self-similarity[END_REF] suggested the threshold is s c ≈ 3.45, above which backscatter term k 4 can be clearly observed and a new energy decay law is derived. Lesieur and Ossia [START_REF] Lesieur | 3d isotropic turbulence at very high reynolds numbers: Edqnm study[END_REF] conducted EDQNM analysis of three-dimensional isotropic turbulence (without helicity) to study this backscatter term issue with cases whose noninteger s is from 3.2 to 3.9. The computation doesn't yield significant observed backscatter term indicating that big eddies are stationary when s < 3, which is consistent with energy decay exponent α e and scaling law Eq.5.1.4. The slope of energy spectrum with non-integer exponent 3 < s < 4 over time evolution is somewhat higher than s and backscatter does gradually observed. The exponent of energy decay law ranges from -4 3 to -1.38. However, they didn't find evidence of backscatter term k 4 above threshold s = 3.45 suggested by Eyink and Thomson. Furthermore, Ossia and Lesieur [START_REF] Ossia | Energy backscatter in large-eddy simulations of threedimensional incompressible isotropic turbulence[END_REF] utilized LES approach and claimed that backscatter term k 4 did not appear at s = 3.5 despite an increase in the energy spectrum. Briard et al. [START_REF] Briard | Spectral modelling for passive scalar dynamics in homogeneous anisotropic turbulence[END_REF] developed an extension of EDQNM to anisotropic flows in relation to the interplay between thirdorder statistics and anisotropy of turbulent flow [START_REF] Cambon | Third-order statistics and the dynamics of strongly anisotropic turbulent flows[END_REF]. See also the analysis of freely decaying rotating turbulence by Delache et al. 2014 [START_REF] Delache | Scale by scale anisotropy in freely decaying rotating turbulence[END_REF].

New challenge

The crucial role of backscatter term in the stationarity of big eddies [START_REF] Schertzer | Comportements auto-similaires en turbulence homogene et isotrope[END_REF], which leads to the reversed energy cascade phenomenon and modifies free energy decay process, has been revealed within the framework of EDQNM. Even though the backscatter term has been observed and quantified its impact in many kinds of research, a strong limitation of EDQNM theory is that it is based on a quasi-normal approximation preventing it from being able to generate intermittency. This fundamental property [START_REF] Sreenivasan | On the fine-scale intermittency of turbulence[END_REF] poses questions about the dynamics of turbulence, in particular, the stationarity of large eddies and scaling law. SGC model preserving the nonlocal triad interaction, which can generate intermittency proven in the last chapter, is applied to investigate the existence of backscatter terms taking intermittency into account. More precisely, the selected NS interactions for SGC model are semi-local interactions rather than strongly nonlocal interactions generating backscatter term. Therefore, this study will provide some insight into whether semi-local interactions can indirectly establish the nonlocal interactions to observe backscatter term in simulation.

SGC model and its backscatter term

Theoretical insights

Due to the triads interaction of SGC model, the evidence for backscatter term is theoretically hypothesized to be derived directly from the nonlinear transfer term.

Since the nonlocal transfer term revealing backscatter term is generated by the nonlocal triad interaction Eq.1.6.15c, the nonlinear transfer term of eddy u i m (Eq.3.1.9) generated by same triad interaction which actually is semi-local is:

T N L U i m (k, t) = u 2i m+1 (k, t) ∧ w 2i m+1 (k, t) + u 2i+1 m+1 (k, t) ∧ w 2i+1 m+1 (k, t). (5.5.1)
Hence, the nonlocal transfer term in wave number space can be obtained from the velocity equation Eq.5.5.1 and its conjugate, as

∂E(k, t) ∂t = u(k, t) • ∂ u(k, t) ∂t + u(k, t) • ∂ u(k, t) ∂t . (5.5.2)
The velocity field in Fourier space has: Then, the conjugation of

∀ k : u(-k) = u(k), (5.5
T N L U i m (k, t) is: T N L U i m (k, t) = u 2i m+1 (k, t) ∧ w 2i m+1 (-k, t) + u 2i+1 m+1 (k, t) ∧ w 2i+1 m+1 (-k, t).
(5.5.6)

The nonlocal transfer term T N L (k, t) producing backscatter term for spatial eddy

u i m is T N L (k i m , t) = u i m (k, t) • u 2i m+1 (k, t) ∧ w 2i m+1 (k, t) + u i m (k, t) • u 2i+1 m+1 (k, t) ∧ w 2i+1 m+1 (k, t)+ u i m (k, t) • u 2i m+1 (k, t) ∧ w 2i m+1 (-k, t) + u i m (k, t) • u 2i+1 m+1 (k, t) ∧ w 2i+1 m+1 (-k, t).
(5.5.7)

After the matrix representation, the nonlocal transfer term T N L (k i m , t) is

T N L (k i m , t) =u i m (k, t) • k m+1 • (|u 2i m+1 (k, t)| 2 -|u 2i+1 m+1 (k, t)| 2 )+ u i m (k, t) • k m+1 • (|u 2i m+1 (k, t)| 2 -|u 2i+1 m+1 (k, t)| 2 ) =2u i m (k, t) • k m+1 • |u 2i m+1 (k, t)| 2 -2u i m (k, t) • k m+1 • |u 2i+1 m+1 (k, t)| 2 (5.5.8) For power spectrum k s (k → 0), T N L (k i m , t) is T N L (k i m , t) ≈2k s/2 m • k m+1 E(k 2i m+1 , t) -2k s/2 m • k m+1 E(k 2i+1 m+1 , t) ≈2k s/2 m • k s+1 m+1 -2k s/2 m • k m+1 E(k 2i+1 m+1 , t).
(5.5.9)

Unfortunately, Eq.5.5.9 fails to provide backscatter term yet. It is important to note that energy lost during the elimination of initial ah-hoc state is not considered, which might be a determining factor. Tree structure of SGC model whose semi-local triad interaction is (k i m , k 2i m+1 , k 2i+1 m+1 ) makes the separating of backscatter term more complex. But Eq.5.5.9 suggests there is a positive term that will slow down the energy decay. The simulation thus will be carried out in the next section to study the energy decay law considering intermittency.

Numerical results

Initial condition for free decaying SGC model

In this part, energy spectrum and initial simulation parameters are presented to investigate the impact of intermittency on backscatter term as well as energy decay law. The energy spectrum has to include both power spectrum k s resulting in backscatter term and k -5 3 inertial spectrum exhibiting extreme intermittency. Energy spectrum shown in Fig 5 .1 is But the energy spectrum for SGC model in simulation is

E(k, 0) =      C s k s , 0 < k m ≤ k L , k -5 3 , k L ≤ k m ≤ k n , (5
E(k, 0) =      C s k s , 1 ≤ k m ≤ k peak , k -5 3 , k peak ≤ k m ≤ k n .
(5.5.11)

It requires the numerical simulated energy spectrum Eq.5.5.11 to be normalized by

E(k peak , 0): E r (k m , t) = E(k m , t) E(k peak , 0) , ( 5
.5.12) which is the time evolution of energy spectrum Eq.5.5.10. Initial velocities for spatial eddies are obtained from energy spectrum Eq 5.5.11, as described in Chapter 4. So do the artificial viscosity, time step, and time units. The classic Euler method is chosen in the simulation of this Chapter.

Corresponding compensated spectrum, as well as the local Reynolds number, can be calculated. The compensated spectrum for inertial spectrum is defined as:

E c (k m /k peak , t) = k 5 3 E r (k m , t). ( 5
Re L (k m /k peak , t) = U L v = (k m /k peak E r (k m , t)) 1/2 vk m /k peak . (5.5.14)
The high Reynolds number implies that the inertial force is many times greater than the viscosity force, so the viscosity force can be neglected and the flow is fully developed turbulent.

Meanwhile, two crucial properties in the dynamic of turbulence remain the same after normalization, which are the exponent of integral scale and energy decay law.

According to the definition of integral scale Eq.4.2.2, integral scale is

L(t) = n m=0 k -1 m K m n m=0 K m = k -1 peak • n m=0 (k m /k peak ) -1 K m /K peak n m=0 K m /K peak , (5.5.15)
which indicates the exponent of integral scale obtained from energy spectrum Eq.5.5.11

and energy spectrum Eq.5.5.10 is equivalent. The deduction of energy decay law is not going to repeat. In this way, the exponent of integral scale and energy decay can be directly estimated from energy spectrum Eq.5.5.11 instead of using the normalized value.

Free energy decay law of spectrum k -5 3

Before studying the role of intermittency, free energy decay law of SGC model without power spectrum is checked to compare with that of energy spectrum Eq.5.5.10.

n peak is assumed as 0 to ensure energy spectrum Eq.5.5.10 without power spectrum, and initial parameters are listed in Table 5.1. Figure 5.17: Exponent of infrared spectrum of Case 5.6.

Reynolds number and compensated spectrum for Case 5.6 after time 500 are infrared ranges and inertial ranges. According to energy decay exponent α e of Case 5.4 and Case 5.5, the energy of spectrum Eq.5.5.10 decays slower than that of case without infrared spectrum, which is also supported by Case 5.6. The infrared range for Case 5.3 is too small, resulting in strange value for α e and integral scale exponent.

Chapter 6 Conclusion and perspectives

In this chapter, the most relevant findings and challenges encountered in this research are reported below, along with a list of future study directions.

With the advancement of turbulence study, the significance of intermittency of turbulence has been progressively realized. However, there are some disagreements regarding its properties and contribution to the mechanism of turbulence, for instance, it was initially assumed to be uniform, and later investigations suggested that it is a log-normal distribution. Complexity analysis of SGC model suggests that it's challenging to improve the execution time as well as memory massively from algorithms. An effective eddy viscosity [START_REF] Kraichnan | Eddy viscosity in two and three dimensions[END_REF][START_REF] Chollet | Parameterization of small scales of threedimensional isotropic turbulence utilizing spectral closures[END_REF] which is a typical two-point closures technique, applies to predict the large scales without information beyond cutoff wave number. However, this method doesn't consider intermittency. It would be meaningful to improve this method, allowing to shorten the running time and memory of large cascade steps to provide more cases to study the role of intermittency in scaling laws.

Regarding intermittency, forced SGC model in the study has been injected the deterministic forcing at the top layer. Recently the linear forcing proposed by Lundgren [START_REF] Lundgren | Linearly forces isotropic turbulence[END_REF] has been considered as a useful alternative [START_REF] Rosales | Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties[END_REF] to classic forcing, as it's simple and can be easily applied to numerical codes in physics, but provides the same results. This linear forcing raises questions about the ability to maintain intermittency, so adding linear forcing to SGC model to test the properties of intermittency will be intriguing. It is important to emphasize that multifractality analysis is performed on the entire time series, including the initial ah-hoc states. This instability state actually has an impact on UM parameters. A hypothesis that UM parameters will be well in agreement with empirical values if this instability state is removed, which has been supported by some cases. Further detailed study will be carried out soon. Another study worth further consideration is reducing the sensitivity of multifractality index estimates caused by algorithms since semi-implicit Euler method is as efficient as classical Euler method. The simulation had finished, as demonstrated in Chapter In conclusion, this research provides some insight on the effective numerical simulation method for SGC model. The existence of spatial-temporal intermittency is confirmed and the multifractality of intermittency questions the log-normal model.

Moreover, some findings-the role of intermittency on backscatter term and scaling laws is investigated. We must acknowledge, nevertheless, that much more research is needed to fully understand the impact of intermittency on energy decay, even if it is challenging, which will be crucial for empirical utilization of energy in the future.

where [ , ] is defined by [E i , E j ] = C k ij E k ; E i is the basis; C ij is the structure constant.

An invariant is the square of the angular momentum M 2

B.3 Algorithms

This part gives the pseudo-code for the numerical methods of the SGC model and the specific Python code for the main function contributing to the velocity field.

B.3.1 Euler method

The pseudo-code is provided first to better comprehend the following python code. for m = 0, . . . , n do 5:

if m = 0 then ▷ Selecting velocity equation at layer 0 6:

Compute Q[0] = U 0 0 (t + ∆t) according to Eq. 3.2. The Python codes for velocities equations using semi-implicit Euler's scheme are given. First of all, the associated Python code for the velocity at top layer u 0 0 Eq.3.2.17 The pseudo-code for the RK4 method is the same as the Euler method except for the equations for computing the velocity field. So the specific Python code instead of the pseudo-code is shown.

The Python code for the velocity at top layer is: Compute Q[l m -1 + i] = U i m (t + ∆t) according to Eq.3.2.15 for i = 0, . . . , K[m] -1 do ▷ i is the location at layer n 13:

# ## v0
Compute Q[l m -1 + i] = U i n (t + ∆t) according to Eq.3.2.14

Figure 1 . 1 :

 11 Figure 1.1: Wind turbines in Aveyron department of the Midi-Pyrenees region of France

  rected. Universal Multifractals (UM) based on multiplicative cascade models (See Chapter 2) have been introduced as a powerful tool to describe intermittency in a scale-invariant (scaling) framework, which emphasized extreme nonlinear variability. The multifractal analysis of empirical data, particularly from lab experiments and atmospheric in-situ/remotely sensed data, has rather consistently generated estimated UM parameters: multifractality index α ≈ 1.5 and mean codimension C 1 ≈ 0.25. NS equation and stochastic cascades are sometimes disagreed, nevertheless. The Scaling Gyroscope Cascade (SGC) model [20] (See 3.1), which is based on a parsimonious discretization of Fourier transform of the Bernoulli's form of NS equation and preserves majority of the fundamental properties of NS equation, is therefore studied since it provides a more close statistical description of natural turbulence and shows intermittency.This thesis initially investigates the efficient numerical approach for SGC model.The intermittency is next verified by energy flux of numerical simulated SGC model, and its temporal multifractality is analyzed by UM framework (See Chapter 4). Free energy decay process considering intermittency is studied by SGC model (See Chapter 5) to fill part of gap between scaling laws obtained by EDQNM and energy transferring process considering intermittency. Here, SGC model is simulated by direct numerical simulation (DNS), which is an ideal method to explore the dynamic and transitional process of turbulence in great detail under the rapid development of computing technology. It gets rid of the limitations of closure problems and analytical approaches and provides a great amount of information on turbulence. Related basic knowledge is provided in the following.

2 2 .Figure 1 . 2 :

 212 Figure 1.2: The Reynolds experiment

of

  Richardson's poem describing turbulence is "Big whorls have little whorls, Which feed on their velocity; And little whorls have lesser whorls, And so on to viscosity." It suggests that turbulence are full of eddies of different sizes. Large-scale eddies in Fig 1.3, which are like tanks and contain energy continuously, obtain energy from the outside and transport energy to small-scale eddies by the inertia of fluid, while small-scale eddies whirling faster will consume energy and turn it into heat by viscous fraction. The inertia range, which depends on Reynolds number, is between largescale energy storage and small-scale energy consumption. The dissipation exists at all scales but could be neglected in the inertial range.

Figure 1 . 3 :

 13 Figure 1.3: The cascade

Lemma 2

 2 The similarity hypothesis: At sufficiently high Reynolds number, the motions of turbulence at small scale (l η < l < l I ) are in a statistically universal form as shown in Fig 1.4, which is uniquely determined by dissipation rate ε and independent of v. l η is Kolmogorov scale (Re η = 1) where dissipation starts.

Figure 1 . 4 :

 14 Figure 1.4: The inertial range. k is wave number, which has k ∝ 1 l .

Figure 1 . 5 :Figure 1 . 6 :

 1516 Figure 1.5: The triad interaction

  into a small fraction, which is non-Gaussian distribution. As shown in Fig 1.7a, vertical velocity [46] 3 shows high values in the time 10h-15h, which is more obviously observed in Fig 1.7b.

Figure 1 . 7 :

 17 Figure 1.7: The wind velocity exhibiting high intermittent. The time unit is an hour and resolution is 1Hz.

Figure 2 . 1 :

 21 Figure 2.1: Two fractal sets

  2.2 provides an illustration of fractal geometry and clarifies fractal dimension D F .

Figure 2 . 2 :

 22 Figure 2.2: Cantor set

Fig 2 . 4 ,

 24 small eddies are either alive or dead, which is consistent with the disappearance of small eddies[START_REF] Batchelor | The nature of turbulent motion at large wave-numbers[END_REF][START_REF] Zhou | Turbulence theories and statistical closure approaches[END_REF] known as internal intermittency of turbulence.

Figure 2 . 4 :

 24 Figure 2.4: β-cascade. Dash line represents dead eddies and the solid one is alive eddies.

.3. 9 )

 9 In the β-model, the activity of field ε depends on scale invariant fractal codimension. The larger co-dimension or high intermittent exponent c indicates more sparseness ( less occupied fraction) of the set (See the black area of Fig 2.5).

( a ) 4 Figure 2 . 5 :

 a425 Figure 2.5: Examples of two dimensional β-model. The resolution is 2 5 . The white area means alive area, while the black area is dead area.

Figure 2 . 6 :

 26 Figure 2.6: The upscaling process of two dimensional β-model whose codimension is 0.2. The maximum resolution is 2 5 .

  Fig 2.8.

Figure 2 . 8 :

 28 Figure 2.8: α-model. Compared with β-model, some eddies act more active and the rest acts less active.

Figure 2 . 9 :

 29 Figure 2.9: The separating process of α-model when λ equals 2.

.5. 2 )Figure 2 .

 22 Figure 2.10: A irregular field showing the codimensions c(γ) is associated to different singularities γ.

2 . 11 .Figure 2 . 11 :

 211211 Figure 2.11: The codimension function c(γ).

  is a space-time gyroscope step-step model depicting evolution of the superposition coupling of daughter eddies in Fig.3.1, which is based on a parsimonious discretization of Bernoulli's form of NS equation. The model doesn't undergo an infinite number of iterations, and the maximum finite cascade step is n.

Figure 3 . 1 :

 31 Figure 3.1: Gyroscope cascade whose separate ratio is λ = 2.

  .1.10) and the orthogonality in Appendix B.2 leads to

  are presented. According to Kolmogorov's local stationary assumption, the inertial range has ∂E(k, t) ∂t = 0, (3.3.1) which suggests Π(k m ) ≈ ϵ(t). (3.3.2) Energy flux, which refers to energy transfer rate from all wave number k < k m to the other wave number k ≥ k m or the cumulative energy transfer rate above the wave number k m , is calculated from

2 .

 2 It reveals that energy fluxes at the initial time period are indistinguishable, whereas they become inconsistent as time goes on. The occurrence of this disparity is expected due to the chaotic nature and introduction of deviations in the velocity field by different numerical algorithms.

Figure 3 . 2 :Fig. 3 . 2 .

 3232 Figure 3.2: Energy flux Π(k 6 ) for Case 3.1

Figure 3 . 3 :

 33 Figure 3.3: UM analysis of energy fluxes Π(k 6 ) for Case 3.1. The UM parameters is obtained from ensemble analysis of samples whose resolution is 512.

Figure 3 . 4 :

 34 Figure 3.4: Comparing commonly used functions in the big O notation when n ≤ 15.

Figure 3 . 5 :

 35 Figure 3.5: The storage space S(n).

Due to Eq 3 . 3 . 2 ,

 332 details of numerical simulations of SGC model are presented, as energy flux equaling kinetic energy dissipation rate reveals significant local fluctuation of turbulence. Energy spectrum is Kolmogorov spectrum assuming C k and ε to be 1.

Energy flux Π(k 6 )Fig 4 .

 64 Fig 4.1i and Fig 4.1j, prevent energy fluxes of Case 4.1.8-Case 4.1.10 from reaching expected time scale 23810. Moreover, the larger forcing f is, the earlier energy flux Π(k 6 ) stops.

Figure 4 . 1 :

 41 Figure 4.1: Energy flux Π(k 6 ) of Case 4.1. Time unit is initial eddy turn over time.Vertical coordinate of energy fluxes is not uniform.

Figure 4 . 2 :Fig. 4 .

 424 Figure 4.2: Energy flux Π(k 6 ) of Case 4.1 getting rid of extremely high values.

Figure 4 . 3 :

 43 Figure 4.3: Energy flux Π(k 5 ) of Case 4.1.1-4.1.3 (t : 1000 -23810).

Fig 4 . 4 ,

 44 probability distribution of Case 4.1 decays much more slowly than corresponding Gaussian distribution with the same mean and variance. The tails of PDF for Case 4.1 are heavier than those of the equivalent Gaussian distribution, which implies that extreme events occur far more frequently than those in Gaussian distribution. Besides, Fig 4.4 demonstrates that the mean and standard deviation resulting in different Gaussian distributions vary depending on the case. Nonetheless, variation in Gaussian distributions is not the focus of this thesis. The extreme events causing "heavy" tail are going to be studied next.

Figure 4 . 4 :

 44 Figure 4.4: Log-log plot of probability P r(∆|Π(k 6 )|) > s) exceeding a fixed reflectivity threshold s of Case 4.1.

Figure 4 . 5 :

 45 Figure 4.5: Temporal evolution of UM parameters and critical moment q s for Case 4.1.6 under size 1024 (a), 512 (b) and 256 (c).

fined range η ≈ 1 .

 1 Nevertheless, this statistical predefined range estimation does not precisely capture multifractality of intermittency, given energy flux varies dramatically over time. For instance, three sub-series of Case 4.1.6 are shown in Fig.4.6 whose UM parameters are determined from the fixed range η ≈ 1. Fig.4.6a and Fig.4.6c display well-fitted UM parameters, whereas UM parameter estimations in Fig. 4.6bdemonstrates that estimation derived from the fixed range is grossly unfit.

Figure 4 . 6 :

 46 Figure 4.6: Bias of UM parameters for Case 4.1.6.

Fig 4 . 7

 47 illustrates the adjusted analysis for subseries 23. Compared to UM parameters estimated from the fixed range in Fig 4.6b, this approach reduce the bias effectively.

Figure 4 . 7 :

 47 Figure 4.7: Adjustment of Sub-series 23 for Case 4.1.6.

  (a) Multifractality index α (b) Mean codimesnion C 1 (c) Critical moment q s

Figure 4 . 8 :

 48 Figure 4.8: Time evolution of UM parameters and critical moment q s for Case 4.1.6 obtained by limiting correlation coefficient. Sub-series size is 512. The horizontal coordinate is sub-series number.

Case 4. 1 UM

 1 Fig.4.9 2 . Due to the range used to estimate UM parameters, difference between Fig.4.9a and Fig.4.9b is observed from trends of α and q s . α decreases with time in Fig.4.9a, but fluctuates steadily in Fig.4.9b. In contrast to trend in Fig.4.9b, critical moment q s in Fig.4.9a varies steadily after unstable initial state. Evidences from frequency histograms and mean values testify the necessity of adjusted range approach. Ensemble analysis gives the same UM parameters and q s for these two different range estimation. The maximum α for Case 4.1.7 in Fig 4.9b does not exceed 2 and the min-

Figure 4 . 9 :

 49 Figure 4.9: Time evolution of UM parameters and critical moment q s for Case 4.1.7.Sample size is 512.

  (a) Fixed range estimation (b) Adjusted range estimation

Figure 4 . 10 :

 410 Figure 4.10: Time evolution of UM parameters and critical moment of Case 4.1.8.

Figure 4 .

 4 Figure 4.11: Sub-series N = 6 in Case 4.1.8.

Fig 4 .

 4 Fig 4.12 is UM analysis for Case 4.1.9 with a large forcing resulting in an enormous energy flux and ending SGC model at time 12000. UM analysis is therefore performed from 0 to 10000. The trends of C 1 and q s are quite similar when comparing Fig 4.12a and Fig 4.12b. Difference in α, which includes trend, the mean, and frequency plot, contributes to alteration of q s . It indicates that the bias for Case 4.1.9 generated by fixed range estimation is non-negligible. In Fig 4.12b, multifractality index α varies from 1.3 to 1.8 and mean codimension C 1 ranges from 0.05 to 0.25. Due to the limitation of C 1 , there is only one sub-series for Case 4.1.9 close to atmospheric

Figure 4 . 12 :

 412 Figure 4.12: Time evolution of UM parameters and critical moment of Case 4.1.9.

Figure 4 . 13 :

 413 Figure 4.13: Time evolution of UM parameters and critical moment of Case 4.1.10.

Figure 4 . 14 :

 414 Figure 4.14: Time evolution of UM parameters and critical moment of Case 4.1.5.

Figure 4 . 15 :

 415 Figure 4.15: Subseries N = 42 in Case 4.1.5.

4. 3 . 4

 34 Time scale closing to atmospheric turbulence Ensemble analysis of sub-series N = 12 -14 of Case 4.1.3 is shown in Fig.4.21 with estimated UM parameters α = 1.54, C 1 = 0.27.

  (a) UM parameters (b) Scaling moment function K(q)

Figure 4 . 21 :

 421 Figure 4.21: Ensemble analysis of Case 4.1.3.

Figure 4 . 22 :

 422 Figure 4.22: Ensemble analysis of Case 4.1.4.

  (a) UM parameters (b) Scaling moment function K(q)

Figure 4 . 23 :

 423 Figure 4.23: Ensemble analysis of Case 4.1.6.

Figure 4 . 24 :

 424 Figure 4.24: Ensemble analysis of Case 4.1.7.

Fig 4 .

 4 Fig 4.25 gives out UM parameters of Case 4.1.8 estimated by ensemble analysis α ≈ 1.52, C 1 ≈ 0.23, which is performed on subseries N = 12 -17.

Figure 4 . 25 :

 425 Figure 4.25: Ensemble analysis of Case 4.1.8.

4. 3 . 5

 35 Summary of cascade step 12Different forcings are injected into SGC model with cascade step n = 12 to study the intermittency of natural turbulence. Fluctuations in energy flux occur and gradually vanish when forcing is between 0.1 and 0.3, whereas extraordinarily large energy flux generated by forcing f larger than 0.8 stopped SGC model, which could be as a result of spatial structure of SGC model amplifying energy exchange. All energy fluxes whose forcing is between 0.4 and 0.8 display highly intermittent. Probability distribution function of these fluctuations has a "heavy" tail, resulting in non-Gaussian distribution and allowing UM framework as an effective tool to study multifractality.According to UM analysis of Case 4.1.6, sample size 512 is chosen to display properties of fluctuations over the entire time series. UM parameters and critical moment q s confirm the difference brought on by fixed range estimation because of different fluctuations for each subseries. Hence, adjusted range estimation method yields better fitted UM parameters. There aren't many subseries whose UM parameters are α ≈ 1.5 and C 1 ≈ 0.25 when forcing is either small or large, but number of subseries closing to empirical UM parameters for atmospheric turbulence increases for the rest cases. Temporal evolution of UM parameters reveals the time of energy flux which is consistent with the atmospheric turbulence, and frequency plots suggests UM parameters for most cases centered on α ≈ 1.5 and C 1 ≈ 0.25. Averaged UM parameters are α ≈ 1.55 ± 0.05 and C 1 ≈ 0.25 ± 0.05 except forcing 0.1. At last, ensemble analyze of energy flux whose parameters consistent to empirical UM parameters is carried out.

Figure 4 . 26 :

 426 Figure 4.26: Energy flux Π(k 8 ) of Case 4.2.

  (a) Case 4.2.1 (t : 500-) (b) Case 4.2.2 (t : 500-) (c) Case 4.2.3 (t : 500-) (d) Case 4.2.4 (t : 500-) (e) Case 4.2.5 (t : 500-) (f) Case 4.2.6 (t : 500-) (g) Case 4.2.7 (t : 500 -19500) (h) Case 4.2.8 (t : 500 -13500) (i) Case 4.2.9 (t : 1000 -20000) (j) Case 4.2.10 (t : 500 -17500) (k) Case 4.2.11 (t : 0 -6800)

Figure 4 . 27 :

 427 Figure 4.27: Energy flux Π(k 8 ) of Case 4.2 removing unstable time.

4. 4 . 2

 42 Non-Gaussian distribution and "heavy" tail After observing fluctuations in energy flux of Case 4.2, probability distribution of ∆|Π(k 8 )| exceeding a fixed reflectivity threshold s is shown in Fig 4.28. (a) Case 4.2.1 (b

Figure 4 . 28 :

 428 Figure 4.28: Log-log plot of probability P r(∆|Π(k 8 )| > s) exceeding a fixed reflectivity threshold s of Case 4.2

  2 are provided below to investigate multifractality of fluctuations over time. Energy flux is analysed in the same way as in cascade step 12 and analysis is presented from the largest forcing to the smallest forcing. UM analysis for Case 4.2.11 is displayed in Fig 4.29.

Figure 4 . 30 :

 430 Figure 4.30: Time evolution of UM parameters and critical moment of Case 4.2.10.

Figure 4 . 31 :

 431 Figure 4.31: Time evolution of UM parameters and critical moment of Case 4.2.9.

Fig 4 .

 4 Fig 4.31a shows a decreasing multifractality index α, whereas Fig 4.31b displays α varying around the average. Comparing these two range estimations, both the average and frequency plot of α change. Similar tendencies are seen in Fig 4.31 for mean codimension C 1 .As a result, q s exhibits significant variation. This significantly

Figure 4 . 33 :

 433 Figure 4.33: Time evolution of UM parameters and critical moment of Case 4.2.7.

Fig. 4 .

 4 Fig.4.33 is the UM analysis of Case 4.2.7. α in Fig 4.33b compared to Fig 4.33a increases, as indicated by mean value. C 1 is greater than 0.8 at the end in Fig 4.33b, whereas it is close to 0.5 in Fig 4.33a. Thus, q s in Fig 4.33b is without a doubt different from those in Fig 4.33a, as shown by mean values and frequency plots, even

Figure 4 . 34 :

 434 Figure 4.34: Time evolution of UM parameters and critical moment of Case 4.2.6.

Figure 4 . 35 : 6 in

 4356 Figure 4.35: Time evolution of UM parameters and critical moment of Case 4.2.5.

Figure 4 . 36 :

 436 Figure 4.36: Time evolution of UM parameters and critical moment of Case 4.2.4.

  (a) Fix range method (b) Adjusted range method

Figure 4 . 37 :

 437 Figure 4.37: Time evolution of UM parameters and critical moment of Case 4.2.3.

  (a) Fix range method (b) Adjusted range method

Figure 4 . 38 :

 438 Figure 4.38: Time evolution of UM parameters and critical moment of Case 4.2.2.

Fig 4 .

 4 Fig 4.39 displays UM analysis of Case 4.2.1. The existence of range estimationgenerated difference is proved by UM analysis comparison. C 1 suggests that energy flux of Case 4.2.1 does not have subseries to approximate empirical UM parameters, despite the frequency plot in Fig 4.39b indicating subseries satisfying α ≈ 1.5.

Figure 4 . 39 :

 439 Figure 4.39: Time evolution of UM parameters and critical moment of Case 4.2.1.

4. 5

 5 Cascade step 15 Numerical simulations of forced SGC model which has a large cascade step n = 15 are displayed with forcing ranging from 0.5 to 1.1. Simulation time grows exponentially with the maximum cascade step n in SGC model because of its spatial structure (see Chapter 3), so simulation is limited by computing devices. The initial condition is presented in Table.4.5.

4. 5 . 1

 51 Intermittency of energy flux Energy flux of Case 4.3 is shown in Fig 4.43 to study the existence of fluctuations. In Case 4.3.1, Case 4.3.2, and Case 4.3.4, the simulation is able to complete computation within expected time frame. For Case 4.3.6 and Case 4.3.7 with large forcing, the simulation yields enormously large energy flux and ends running. Regardless, extreme values in either the initial unstable state or the end state, make it difficult to clearly observe fluctuations in Fig 4.43.

Figure 4 . 43 :

 443 Figure 4.43: Energy flux Π(k 9 ) whose cascade step n = 15.

Figure 4 . 44 :

 444 Figure 4.44: Energy flux Π(k 9 ) of Case 4.3.

  Fig 4.45 verify "heavy" tail for Case 4.3 by comparing probability distribution of P r(∆|Π(k 9 )| > s) with Gaussian distribution that has the same mean and variance.

Figure 4 . 45 :

 445 Figure 4.45: Log-log plot of probability P r(∆|Π(k 9 )| > s) exceeding a fixed reflectivity threshold s of Case 4.3

Figure 1 . 4 :Fig. 1 .

 141 Figure 1.4: Time evolution of UM parameters and critical moment of Case 3.7.

Figure 4 . 46 :

 446 Figure 4.46: Time evolution of UM parameters and critical moment of Case 4.3.7.

  Analysis of Case 4.3.6 is shown in Fig 4.47, whose time scale reaches its maximum around 20000 and is adequate to study intermittency.

Figure 4 . 47 :

 447 Figure 4.47: Time evolution of UM parameters and critical moment of Case 4.3.6.

Figure 4 . 48 :

 448 Figure 4.48: Time evolution of UM parameters and critical moment of Case 4.3.5.

Fig 4 .

 4 Fig 4.48 illustrates UM analysis of Case 4.3.5, which demonstrates the observed changes in α and q s , including trends, frequency plots, and mean values. For instance, the trend of α after N = 20 in Fig 4.48b is much more stable than that in Fig 4.48a. This confirms the difference caused by range estimation. Besides, time evolution

Figure 4 . 50 :

 450 Figure 4.50: Time evolution of UM parameters and critical moment of Case 4.3.3.

Figure 4 . 51 :Fig. 4 .

 4514 Figure 4.51: Time evolution of UM parameters and critical moment of Case 4.3.2.

Figure 4 . 52 :

 452 Figure 4.52: Time evolution of UM parameters and critical moment of Case 4.3.1.

Fig 4 . 53 .

 453 Fig 4.55b.

Figure 4 . 53 :

 453 Figure 4.53: Ensemble analysis of Case 4.3.4.

Fig 4 .

 4 Fig 4.54 is ensemble analysis of subseries N = 30 -38 in Case 4.3.5. UM parameters are α = 1.47 and C 1 = 0.27 with critical moment function is q s = 2.96.

Figure 4 . 54 :

 454 Figure 4.54: Ensemble analysis of Case 4.3.5.

Figure 4 . 55 :

 455 Figure 4.55: Ensemble analysis of Case 4.3.6.

. 3 )

 3 and Eq3.1.10 considering the complex conjugation symmetry of SGC model yields∀ k : w(-k) = w(k);(5.5.4)∀ k : u ∧ w = u ∧ w.(5.5.5)

  .5.10) where C s is k -(s+ 5 3 ) L required by connected wave number k L . Comte-Bellot, G and Corrsin, S [23] also proposed energy spectrum Eq.5.5.10 and studied the connection between -10/7 decay law and integral scale 2/7 law.

Figure 5 . 1 :

 51 Figure 5.1: Energy spetcrum Eq 5.5.10.

Figure 5 . 2 :

 52 Figure 5.2: Log-log plot of energy E(t) of Case 5.1 and Case 5.2. Time unit of time t is the initial large eddy turnover time. K is system energy E(t).

Fig. 5 .

 5 Fig.5.2 presents the log-log plot of system energy for these two cases without power spectrum. Apart from the initial unstable phase, exponent of energy decay law at time t = 30 -300 in Fig.5.2a is α e = 1.41, and α e = 1.47 at time t between 50 -300 in Fig.5.2b. The following are the other characteristics of same time period.

Figure 5 . 3 :Fig. 5 .

 535 Figure 5.3: Local Reynolds number Re of Case 5.1 and Case 5.2.

Figure 5 . 4 : 2 7. 1 and

 5421 Figure 5.4: Energy spectrum of Case 5.1 and Case 5.2, as well as their compensated spectrum.

Figure 5 . 8 :

 58 Figure 5.8: Log-log plot of system energy and integral scale of Case 5.3

( a )Figure 5 . 9 :

 a59 Figure 5.9: Ensemble averaged spectrum E r (k m , t) of Case 5.4.

Figure 5 . 10 :

 510 Figure 5.10: Ak 4 energy spectrum of Case 5.4

Figure 5 . 11 :

 511 Figure 5.11: Compensated spectrum and Reynolds number of Case 5.4 from time 310 to 400.

( a )Figure 5 . 12 :

 a512 Figure 5.12: The Log-log plot of energy and integral scale of Case 5.4

( a )Figure 5 . 13 :

 a513 Figure 5.13: Energy spectrum of Case 5.5.

Fig 5 .

 5 Fig 5.13 is time evolution of ensemble-averaged energy spectrum of Case 5.5 which is similar to Case 5.3 that spectrum coefficient in infrared range rises due to the reversed energy and has vanished between time 200 and 300 in Fig 5.13c.

Figure 5 . 16 :

 516 Figure 5.16: Energy spectrum of Case 5.6.

4 ,

 4 and further in-depth research will be summarized and published in related paper.Large infrared range research will provide a clear insight into the impact of semilocal interactions of SGC model on basckscatter term and scaling laws. Until now, more cases with large infrared range are limited by the drawbacks of cutoff wave number method, as well as successfully assessing the role of backscatter term in slowing energy decay by Loitsiansky exponent over large time scales and a wide range of k 4 energy spectrum. Besides, helical modes decompositions which is explored for the non-isotropic impacts of solid-body rotation on homogeneous turbulence, appears consistent with a trend toward two-dimensionalityas a result of allowing a departure from isotropy via external rotation effects and the loss of spectral energy from the parallel to the normal wave vectors (relative to the rotation axis). It will provide new perspectives on the dynamics of SGC model, such as the manner in which energy is transferred.

B. 2 Figure B. 2 :

 22 Fig.B.2a. In order to make Eq. 3.1.9 link to gyroscope equation Eq. B.1.1, the M i m and Ω i m has the orthogonal decomposition ( u i m , u 2i m+1 , u 2i+1 m+1 ) Fig.B.2b and (k 2i m+1 , ω 2i+1 m+1 , ω 2i m+1 ) Fig.B.2c. These orthogonality conditions yield

Algorithm 1 1 : 2 : 3 :

 1123 The Euler method Input: The initial velocity U , The wave number K, The time steps ∆t, The viscosity v, The computation time T , The maximum cascade step n, The ratio l Output: The velocities U T for all locations at time 0 -T function computing velocities (U (t + ∆t)) ▷ By the former velocity U (t) Set L as the size of U ▷ Number of eddies at all layers Set the working list Q = [0 . . . L -1] 4:

3 7 :

 7 else if m < n then ▷ Selecting velocity equation at layer m # ## k0 is the wave number at top layer # ## k1 is the wave number at cascade step 1 def Ve lo ci ty_ fu nc ti on 0 ( v0 , v1 , v2 , k1 , k0 ,v , time_step , Force ) : a = 1 -v * np . power ( k0 ,2) * time_step d = np . power ( v1 ,2) -np . power ( v2 ,2) new_velocity = v0 * a +( k1 * d + Force ) * time_step return new_velocity Code Listing B.1: Computing velocity at top layer u 0 0 (t + ∆t) by Euler method The corresponded Python code for the velocities u i n which are located at last layer Eq.3.2.4 is presented:# ## i is the location in ( n ) th step # ## v_n0is the initial velocity # ## a_n is the ancestor eddy # ## k_n is the wave number at layer n def Ve lo ci ty_ fu nc ti on 2 ( v_n0 ,i , a_n , k_n ,v , time_step , Force ) : a =( -1) **( i +1) b = a * k_n * v_n0 * a_n c = v * np . power ( k_n ,2) B =1 -c * time_step new_velocity = v_n0 * B +( b + Force ) * time_step return new_velocity Code Listing B.2: Computing velocity u i n (t + ∆t) by Euler method The Python code for the velocities u i m whose cascade step m is between 1 and n -1 Eq.3.2.5 are simulated by : # ## i is the location at ( s ) th step # ## v_s0 is the initial velocity # ## a_s is it 's ancestor eddy # ## v_s1 , v_s2 is it 's sub -eddies # ## v is the viscosity # ## k_s is the wave number at cascade step s # ## k_s1 is the wave number at cascade step s +1 def Ve lo ci ty_ fu nc ti on 1 ( v_s0 ,i , a_s , v_s1 , v_s2 , k_s1 , k_s ,v , time_step , Force ) : a =( -1) **( i +1) b = np . power ( v_s1 ,2) -np . power ( v_s1 ,2) c = v * k_s **2 C =1 -c * time_step + a * k_s * a_s * time_step new_velocity = v_s0 * C +( k_s1 * b + Force ) * time_step return new_velocity Code Listing B.3: Computing velocity u i m (t + ∆t)(1 ≤ m ≤ n -1) by Euler method Python code for Function 1 of Algorithm 1 computing U (t + ∆t) is presented in the following. # ## n is the maximum cascade step , l is the separating ratio # ## wave_number is the wave number for all cascade steps # ## v is the viscosity # ## Forme r_velo cities is the velocity at time t def Time_function (n , l , wave_number , Former_velocities , v , time_step , Force ) : length = len ( Fo rmer_v elocit ies ) # ## the number of eddies needs to compute New_velocities =[[]]* length # ## Prepare array for new velocities for layer in range ( n +1) : # ## Selecting the velocities equation by cascade step

2 d 7 :

 27 v0 is the initial value of the top layer # ## v is the viscosity # ## v1 , v2 is v0 's sub -eddies # ## k0 is the wave number at top layer # ## k1 is the wave number at cascade step 1 def Ve lo ci ty_ fu nc ti on 0 ( v0 , v1 , v2 , k1 , k0 ,v , time_step , Force ) : a = 1+ v * np . power ( k0 ,2) * time_step d = Square ( v1 ) -Square ( v2 ) b = v0 +( k1 * d + Force ) * time_step new_velocity = b / a return new_velocity Code Listing B.5: Computing velocity at top layer u 0 0 (t + ∆t) by Semi-implicit Euler methodThe corresponded Python code for the velocities u i n which are located at last layer Eq.3.2.18 is presented:# ## i is the location in ( n ) th step # ##v_n0 is the initial velocity # ## a_n is the ancestor eddy # ## k_n is the wave number at layer n def Ve lo ci ty_ fu nc ti on 2 ( v_n0 ,i , a_n , k_n ,v , time_step , Force ) : a =( -1) ** i b = a * k_n * v_n0 * a_n c = v * k_n **2 d = v_n0 +( b + Force ) * time_step B = 1 + c * time_step new_velocity = d / B return new_velocity Code Listing B.6: Computing velocity u i n (t + ∆t) by Semi-implicit Euler method The Python code for the velocities u i m whose cascade step m is between 1 and n -1 Eq.3.2.19 are simulated by : # ## i is the location at ( s ) th step # ## v_s0 is the initial velocity # ## a_s is it 's ancestor eddy # ## v_s1 , v_s2 is it 's sub -eddies # ## v is the viscosity # ## k_s is the wave number at cascade step s # ## k_s1 is the wave number at cascade step s +1 def Ve lo ci ty_ fu nc ti on 1 ( v_s0 ,i , a_s , v_s1 , v_s2 , k_s1 , k_s ,v , time_step , Force ) : a =( -1) ** i b = Square ( v_s1 ) -Square ( v_s2 ) A = k_s1 * b + a * k_s * v_s0 * a_s c = v * k_s **Computing velocity u i m (t + ∆t)(1 ≤ m ≤ n -1) by Semi-implicit Euler method B.3.2 The fourth order Runge-Kutta method

1 : 2 : 3 : 4 : 5 : 6 : 7 : 8 :

 12345678 is the initial value of the top layer # ## v is the viscosity # ## v1 , v2 is v0 's sub -eddies # ## k0 is the wave number at top layer # ## k1 is the wave number at cascade step 1 def Runge_0 ( v0 , v1 , v2 , k1 , k0 ,v , time_step , Force ) : a = -v * np . power ( k0 ,2) d = np . power ( v1 ,2) -np . power ( v2 ,2) k_1 = v0 * a +( k1 * d + Force ) # ## the derivative function y_1 = time_step * k_1 # ## y_1 = dx * f ( x , y ) # ## y2 = dx * f ( x + 0.5 * dx , y + 0.5 * y1 ) c = v0 +0.5* y_1 k_2 = c * a +( k1 * d + Force ) y_2 = time_step * k_2 # ## y3 = dx * f ( x + 0.5 * dx , y + 0.5 * y2 ) f = v0 +0.5* y_2 k_3 = f * a +( k1 * d + Force ) y_3 = time_step * k_3 # ## y4 = dx * f ( x + dx , y + y3 ) g = v0 + y_3 k_4 = g * a +( k1 * d + Force ) y_4 = time_step * k_4 # ## the new Velocity y + ( y1 + 2 * y2 + 2 * y3 + y4 ) / 6. V = v0 +( y_1 + 2 * y_2 + 2 * y_3 + y_4 ) / 6. return V Code Listing B.8: Computing velocity u 0 0 (t + ∆t) simulated by RK4 methodThe Python code for the velocity u i n is:# ## i is the location in ( n ) th step # ##v_n0 is the initial velocity # ## a_n is the ancestor eddy # ## k_n is the wave number at layer n def Runge_2 ( v_n0 ,i , a_n , k_n ,v , time_step , Force ) : # ## y1 = dx * f ( x , y ) a =( -1) **( i +1) b = a * k_n * a_n c = -v * k_n **2+ b k_1 = v_n0 * c + Force y_1 = time_step * k_1 # ## y2 = dx * f ( x + 0.5 * dx , y + 0.5 * y1 ) C = v_n0 +0.5* y_1 k_2 = C * c + Force y_2 = time_step * k_2 # ## y3 = dx * f ( x + 0.5 * dx , y + 0.5 * y2 ) f = v_n0 +0.5* y_2 k_3 = f * c + Force y_3 = time_step * k_3 # ## y4 = dx * f ( x + dx , y + y3 ) g = v_n0 + y_3 k_4 = g * c + Force y_4 = time_step * k_4 # ## y + ( y1 + 2 * y2 + 2 * y3 + y4 ) / 6. V = v_n0 +( y_1 + 2 * y_2 + 2 * y_3 + y_4 ) / 6. return V Code Listing B.9: Computing velocities u i n by RK4 method The Python code for the velocities u i m from layer 1 to n -1 is: # ## i is the location at ( s ) th step # ## v_s0 is the initial velocity # ## a_s is it 's ancestor eddy # ## v_s1 , v_s2 is it 's sub -eddies # ## v is the viscosity # ## k_s is the wave number at cascade step s # ## k_s1 is the wave number at cascade step s +1 def Runge_1 ( v_s0 ,i , a_s , v_s1 , v_s2 , k_s1 , k_s ,v , time_step , Force ) : # ## y1 = dx * f (x , y ) a =( -1) **( i +1) b = np . power ( v_s1 ,2) -np . power ( v_s2 ,2) A = k_s1 * b c = -v * k_s **2+ a * k_s * a_s k_1 = v_s0 * c + A + Force y_1 = time_step * k_1 # ## y2 = dx * f ( x + 0.5 * dx ,y + 0.5 * y1 ) C = v_s0 +0.5* y_1 k_2 = C * c + A + Force y_2 = time_step * k_2 # ## y3 = dx * f ( x + 0.5 * dx , y + 0.5 * y2 ) f = v_s0 +0.5* y_2 k_3 = f * c + A + Force y_3 = time_step * k_3 Input The initial velocity U , The wave number K, The time steps ∆t, The viscosity v, The computation time T , The maximum cascade step n, The ratio l Output The velocities U T for all locations at time 0 -T function Computing velocities(U (t + ∆t)) ▷ By the former velocities U (t), U (t -∆t) Set L as the size of U ▷ Number of eddies for all layers Set the working list Q = [0 . . . L -1] ▷ Velocities U (t + ∆t) for m = 0, . . . , n do ▷ Selecting velocities equations by cascade step m if m = 0 then ▷ Using velocity equation at top layer Compute Q[0] = U 0 0 (t + ∆t) according to Eq.3.2.13 ▷ The velocity at layer 0 else if m < n then ▷ Selecting velocity equation at layer m for i = 0, . . . , K[m] -1 do ▷ i is the location at layer m 9:

19 :Figure B. 6 :

 196 Figure B.5: Slaved ABF method

  Slaved ABF method: Velocities U i m for all layers (0 ≤ m ≤ n) at each time step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207 Cascade Step in most cases. The word layer is used as a synonym in explaining SGC model. Besdies, it is shell in GOY model (Gledzer 1973, Ohkitani and Yamada 1989) . The fact is that all of these terms refer to the same thing n.One of UM parameters, the codimension of singularity of the mean field, which

	ℓ(t) A rather arbitrary scale of energy-containing eddies in the energy spectrum F (x) Dimensionless functions of the dimensionless wavevector x in the energy spec-r(N ) Similarity ratio or the inverse resolution, 27 δ i,j The Kronecker's δ, 58 t f Maximum time, 82 DIA Direct Interaction Approximation
	function, 151 trum function, 151 R α D The modified Reynolds number, 152 η * Artificial dissipation cascade step , 69 equals its singularity. , 50 U Velocity, 6 EDQNM Eddy Damped Quasi-Normal Markovian
	ℓ d,α D (t) The dissipation scale related to R α D , 152 G(x) Dimensionless functions of the dimensionless wavevector x in the energy spec-R ij (l) Correlation function , 9 γ + The positive singularity , 45 C 1e Mean codimension obtained from ensemble analysis , 88 v Viscosity, 6 UM Universal Multifractals Key words-turbulence, intermittency; multifractality; the SGC model; free energy decay law; backscatter term, nonlocal triad interactions; ∞ Infinite, 14 trum function, 151 s The exponent of stationary of big eddies, 151 γ -The negative singularity , 45 c ij Sub-codimensions , 46 Operations SGC Scaling Gyroscope Cascade Chapter 1
	λ I(t) Loitsiansky integral, 15 Separate ratio , 30 S q (l) The moments function, 10 λ k The arbitrary nonlocalness parameter , 59 D The diameter of the pipe, 6 Integral, 14 DNS Direct Numerical Simulation
	µ K v turb (k) k-independent viscosity, 21 Increment of the energy flux at cascade step n, 30 kinetic energy, 14 λ k Parameter to separate the local and nonlocal interactions , 21 E(t) System energy, 82 Product , 30 GSI Generalized Scale Invariance Introduction
	I would like to dedicate this thesis to my family, especially my grandmother who cared about me at her last moment. . 200 B.11 RK4 method. Velocities U i m for all layers (0 ≤ m ≤ n) at each time step201 B.12 Computing velocity u 0 0 (t + ∆t) by slaved ABF method . . . . . . . . 205 B.13 Computing velocities u i n (t + ∆t) by slaved ABF method . . . . . . . 206 α e Exponent of energy decay law, 15 α l Exponent of integral scale, 16 f Body force, 7 k Wave number, 9 p Pressure, 7 r Location vector, 7 t Time, 7 u E(k) Energy spectrum, 9 ∆t The time step, 62 S(n) The storage space, 72 NS Navier-Stokes Velocity which has vector, 7 ρ Density, 7 σ 2 The variance, 22 ε Energy flux rate or energy flux density, 8 u Velocity u in Fourier space, 17 A n ε n = λ nc , 34 A c n ε n = 0, 34 c Co-dimension, 28 C s The prefactor of stationary of big eddies, 151 D Euclidean dimension, 28 D F Fractal dimension, 26 L(t) Integral scale, 16 l 0 Initial scale, 28 β * The power spectrum exponent, 13 Ω The rotation , 58 ω The vorticity , 58 e , 18 M The angular momentum, 58 P ij , 18 C 1 Averaged mean codimension obtained from all sub-series , 88 q s The mean critical moment obtained from all sub-series , 88 τ Eddy turnover time, 13 τ e Initial large-eddy turnover time, 81 θ k,p,q The correlation relaxation time of the triad k, p, q, 155 ε n The energy flux density at step n , 30 L(t) Integral scale, 81 n + The number of more active steps , 45 n -The number of less active steps , 45 P (∇) The projector of nabla operator, 58 q se Critical moment obtained from ensemble analysis, 88 Re Reynolds number, 7 ∝ Proportional to, 13 curl The curl , 58 grad The gradient , 58 Abbreviations RTE Gestionnaire du Réseau de Transport d'Electricité ABF Adams-Bashforth (ABF) Method l , 28 α One of UM parameters, multifractality index , 50 α * kinematic pressure, 58 α e Multifractality index obtained from ensemble analysis , 88 α v The exponent of the velocity v(t), 151 ∇ The nabla operator , 58 α Averaged multifractality index obtained from all sub-series , 88 ϵ The average of energy transfer rate , 13 u Conjugate of u, 17 Scale of motion, 7 Ω The divergence , 58 • Scalar product , 7 ∇ Gradient, 7 TKE Turbulence Kinetic Energy DTM Double Trace Moment RK4 Fourth-order Runge-Kutta Method E L ∆ The Laplace Operator , 58 ODE Ordinary Differential Equation 1.1 Background
	q δ r u(l) The velocity increment, 10 The moments order, 10 F Fractal set, 27 δ i,j The Kronecker's δ, 18 T (n) Statement frequency, 75 QNA Quasi-Normal Approximation
	viii xxix xxx xxxi xxxii xxxiv xxxv

ix Code Listings B.1 Computing velocity at top layer u 0 0 (t + ∆t) by Euler method . . . . . 191 B.2 Computing velocity u i n (t + ∆t) by Euler method . . . . . . . . . . . . 192 B.3 Computing velocity u i m (t + ∆t)(1 ≤ m ≤ n -1) by Euler method . . 192 B.4 Euler method. The velocities U i m for all layers (0 ≤ m ≤ n) at each time step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 B.5 Computing velocity at top layer u 0 0 (t+∆t) by Semi-implicit Euler method196 B.6 Computing velocity u i n (t + ∆t) by Semi-implicit Euler method . . . . 196 B.7 Computing velocity u i m (t + ∆t)(1 ≤ m ≤ n -1) by Semi-implicit Euler method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 B.8 Computing velocity u 0 0 (t + ∆t) simulated by RK4 method . . . . . . . 198 B.9 Computing velocities u i n by RK4 method . . . . . . . . . . . . . . . . 199 B.10 Computing velocities u i m (t + ∆t)(1 ≤ m ≤ n -1) simulated by RK4 method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.14 Computing velocities u i m (t + ∆t) (1 ≤ m ≤ n -1) by slaved ABF method206 xxv B.15 xxvi List of Algorithms 1 The Euler method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 2 The slaved Adams-Bashforth Method . . . . . . . . . . . . . . . . . . 203 xxvii Notation and Abbreviations xxviii Notations Notations α ℓ The exponent of scaling law of scale ℓ(t), 151 α D The dissipativity is the power of negative Laplacian when v ∂ 2 ∂x 2 is replaced by -v(-∂ 2 ∂x 2 ) α D , 152 I The demarcation scale between the anisotropic large eddy and the isotropic small eddies , 12 l n Scale at steps n , 28 L s (t) Saffman invariant, 16 l η Kolmogorov scale/ dissipation scale, 12 N The number of parts that the object exactly decomposed into, 27 n ζ(q) Exponent of the qth-order moment structure function, 10 c(γ) Co-dimension function , 47 xxxiii C 1 m (t) Energy at wave number k m , 82 H * A measure of correlation, 13 K(q) The moment scaling function, 48

  Now, initial values for Case 3.1, which is used for the comparison, are shown in Table3.1. Time step ∆t should be small enough to guarantee the stability of numerical simulation methods. RK4 method has a broader absolute stability zone than Euler method, whereas slaved ABF method is strongly stable. As a result, time step of these approaches is uniformly assumed by the stability of Euler method. It is also limited by cascade steps. Hence, the restriction on time step ∆t is that it should approximately obey |1vk 2 n ∆t| < 1, since the main dissipation term is viscous term and the convection term can be neglected. The running time Tmax is set as sufficient enough time to investigate the difference in results caused by numerical methods.

.3.7) Related artificial dissipation cascade step η * is generally assumed as two steps before the maximum cascade step n.

Table 3 .

 3

		1: Case 3.1
	Case n η *	∆t	Tmax f
	3.1 12 10 0.0002 30000 0
	Energy flux Π(k 6 ) calculated from these different numerical methods is displayed
	in Fig 3.

Table 3 .

 3 2: The space complexity for numerical methods.

	Lines	S(n)	Space Complexity
	Line 2,6,9,13	1	O(1)

Table

.

3.2. There are a total of 2 n+1 -1 spatial eddies in SGC model, all of whose velocities must be calculated. Storage space S(n) of Function 1 in Algorithm 1 becomes 2 n+2 -1 after accounting for extra memory occupied in computation as well as new velocities needed to be written out. Except for velocities equations, Euler method including classic Euler method and semi-implicit Euler method share the same algorithm with RK4 method, hence space complexity of these one-step

methods is O(2 n ). Storage space S(n) of Function 1 in Algorithm 2 is 3 • 2 n+1 -2

and complexity is O(2n) since slaved ABF method, in comparison to those one-step approaches, requires an additional set of inputs U (t -∆t). Table

.

3.2 indicates that these four numerical methods have the same space complexity O(2 n ) due to spatial structure of SGC model.

Table 3 .

 3 .3 displaying time complexity of these algorithms. Simple operations, such as addition, have an execution number of 1 and time complexity of O(1). 3: The time complexity for algorithms.SGC model that preserves the most fundamental properties of NS equation, including triad interactions, is applied to characterize the properties of natural turbulent flows. In this chapter, four numerical methods-classic Euler approach, RK4 method, slaved ABF method (considering fast dissipation at high wave numbers), and semiimplicit Euler method (widely applied to the motion of rigid body) -are evaluated to determine the effective numerical simulation for SGC model and to provide simulated transient motion of eddies on the full-scale range for next chapters. Algorithms analysis consists of UM analysis, complexity analysis and empirical analysis. However, slaved ABF method per time step requires more memory than the other two methods. A more thorough empirical posterior analysis reveals that maximum memory for Case 3.1 which has a large enough time scale, is almost the same for all three methods because the output velocity field is preserved at all time steps. Euler method takes the least amount of time in terms of user time.

	Executions Lines	T (n)	Time Complexity
	Line 2,3,6,9,13	1	O(1)
	Lines 8-10	2 m	O(2 m )
	Lines 12-14	2 n	O(2 n )
	Function 1	2 n+1 + 1	O(2 n )
	Both space complexity and time complexity defined in terms of cascade step n are
	O(2 n ) for these algorithms, suggesting that spatial structure of SGC model determines
	its complexity. Time and storage space required to simulate SGC model increases
	exponentially when cascade step n grows, and numerical simulation approaches are
	incapable of improving SGC model's computation. Complexity analysis, however,

number of computing velocities for each cascade step m is 2 m . Subsequently, all eddies velocities U (t + ∆t) takes T (n) = 2 n+1 + 1 iterations, which is the sum of all steps executions, and time complexity is O(2 n ). is a rather rough evaluation of algorithms. A more thorough posterior analysis is presented to enhance the overall evaluation.

3.5 Summary of Chapter 3

UM parameters of Case 3.1 demonstrate that multifractality index of intermittency using classic Euler scheme is close to that of the fourth-order Runge Kutta scheme, while multifractality index of semi-implicit Euler method is close to the result obtained by Slaved ABH method. The definition of multifractality index makes it common for numerical approaches to result in relatively small difference. Priori analysis reveals that spatial-temporal complexity of these algorithms computing velocities for each time step is O(2 n ) because of spatial structure of SGC model. Execution time T (n) and storage space S(n) of SGC model grow exponentially when cascade step n increases.

Table 4 .

 4

1 lists key parameters of forced SGC model with a cascade step of 12 ensuring a sufficient inertial range in cascade model. External forcing f is limited to the range of 0.1 to 1.0 as a result of assumptions of C k and ε in energy spectrum. Energy flux whose step is m has to be in inertial range. The time in simulation is displayed in terms of initial eddy turnover time τ e with a large time scale t f reaching 10 4 .

Table 4 .

 4 1: Parameters of Case 4.1 whose cascade step is n = 12.

	Case	Case 4.1 n η * 12 10 0.0002 6 ∆t m f t f
	4.1.1	0.1	23810
	4.1.2	0.2	23810
	4.1.3	0.3	23810
	4.1.4	0.4	23810
	4.1.5	0.5	23810
	4.1.6	0.6	23810
	4.1.7	0.7	23810
	4.1.8	0.8	18000
	4.1.9	0.9	12000
	4.1.10	1	11000

Table 4 .

 4 2: Frequency of UM parameters of Case 4.1 closing to atmospheric turbulence.

				Case 4.1
		n	η *	∆t	m resolution
		12	10	0.0002	6	512
	Case	f	t f	N α	N C 1 α	C 1	α e	C 1e
	4.1.1 0.1 23810	17	3	1.13 0.12 1.44 0.08
	4.1.2 0.2 23810	12	5	1.57 0.25 1.27 0.09
	4.1.3 0.3 23810	11	7	1.47 0.23 1.28 0.09
	4.1.4 0.4 23810	11	10 1.59 0.29 1.51 0.17
	4.1.5 0.5 23810	16	10 1.51 0.3 1.26 0.15
	4.1.6 0.6 23810	17	14 1.54 0.24 1.35 0.18
	4.1.7 0.7 23810	14	9	1.54 0.23 1.13 0.12
	4.1.8 0.8 18000	14	9	1.5 0.22 1.42 0.19
	4.1.9 0.9 10000	11	1	1.53 0.1 1.58 0.11
	4.1.10 1 10000	6	0	1.59 0.1 1.66 0.1

No sub-series consistent with atmospheric turbulence exists due to the limitation of N C 1 when external forcing is large, as in the range 0.9-1, or low, as in the range f ≤ 0.3. Averaged UM parameters of Cases 4.1.4-4.1.8, along with N α and N C 1 , imply that the middle range of forcing is appropriate to describe intermittency. UM parameters become stable once initial instability period is finished and injected forcing have maintained the system. Except for Case 4.1, averaged multifractality index α is 1.55 ± 0.05, on the other hand, most of C 1 is 0.25 ± 0.05. However, UM parameters obtained by ensemble analysis are inaccurate, since it is performed on all sub-series including initial unstable time, averaged UM parameters instead of ensemble values are sufficient to study the entire time series.

Table 4 .

 4 1 to 1.1 because of larger cascade step. Time step and dissipation step have both been reset. Initial values are presented in Table 4.3 as well as extended cascade step that was used to study energy flux. The 3: Case 4.2 whose cascade step is n = 14

		Case 4.2 n η * 14 12 0.0001 8 ∆t m
	Case 4.2.1 4.2.2 4.2.3 4.2.4 4.2.5 4.2.6 4.2.7 4.2.8 4.2.9 4.2.10 4.2.11	f 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1	τ 25001 25001 25001 25001 25001 25001 21371 14000 20974 19000 7000

study is divided into three parts: intermittency, non-Gaussian distribution, and time evolution of multifractality.

4.4.1 Intermittency of energy flux

Energy flux Π(k 8 ) of Case 4.2 is displayed in Fig 4.26. (a) Case 4.2.1 (b) Case 4.2.2 (c) Case 4.2.3 (d) Case 4.2.4 (e) Case 4.2.5 (f) Case 4.4.2.6 (g) Case 4.2.7 (h) Case 4.2.8 (i) Case 4.2.9 (j) Case 4.2.10 (k) Case 4.2.11

Table 4 .

 4 4. Case 4.2.7 -Case 4.2.11, whose forcing is large, doesn't reach expected time because of enormous energy flux. And few subseries satisfy C 1 ≈ 0.25 with the exception of Case 4.2.9 resulting in mean values below 0.25. Case 4.2.1 and Case 4.2.2, on the other hand, have relatively few cases satisfying C 1 ≈ 0.25 due to small forcings, as indicated by mean values. As for Case 4.2.3-Case 4.2.6, there are many subseries whose UM parameters are close to α ≈ 1.5 and C 1 ≈ 0.25. Case 4.2 with sample size 512 has many subseries α > 2, which is caused by extreme values dominating scaling properties and leads to α being larger than 1.5.Since averaged values are obtained over all subseries, including unstable state at the beginning or the end, it cannot be used to determine whether or not the entire time series is close to atmospheric turbulence. Time period of Case 4.2 that is consistent with atmospheric turbulence is displayed in the following.Due to UM analysis for time evolution of fluctuations in the last section, this part presents ensemble analysis for the time exhibiting atmospheric turbulence-like properties. It starts with cases with small forcing. Here it has to point out that time scale for ensemble analysis doesn't necessarily have to be slot subseries from time evolution UM analysis. Due to the limitation of C 1 , cases whose time scales reflect 4.2 is α ≈ 1.66 ± 0.06, as there are several values α > 2 caused by the extreme values in the field dominating scaling properties. Frequency plots, as well as Table4.4, reveal that there are more subseries with UM parameters close to α ≈ 1.5, C 1 ≈ 0.25 for middle-range forcings. Mean codimension C 1 is 0.25 ± 0.05 with the exception of small forcing Case 4.2.1-Case 4.2.2, large forcing Case 4.2.10-Case 4.2.11, and Case 4.2.8. Temporal evolution provides time periods and ensemble analysis confirms that part of time evolution is in agreement with atmospheric turbulence.

	Table 4.4: Frequency of UM parameters of Case 4.2 closing to atmospheric
				turbulence.
				Case 4.2
		n	η *	∆t	m resolution
		14	12	0.0001	8	512
	Case	f	t f	N α	N C 1 α	C 1	α e	C 1e
	4.2.1 0.1 25001	16	1	1.58 0.12 1.62 0.09
	4.2.2 0.2 25001	16	2	1.62 0.14 1.63 0.09
	4.2.3 0.3 25001	14	17 1.68 0.25 1.3 0.12
	4.2.4 0.4 25001	11	10 1.7 0.31 1.38 0.1
	4.2.5 0.5 25001	11	9	1.61 0.23 1.37 0.1
	4.2.6 0.6 25001	17	12 1.54 0.19 1.29 0.11
	4.2.7 0.7 21000	19	4	1.65 0.23 1.49 0.17
	4.2.8 0.8 13500	9	3	1.6 0.16 1.35 0.12
	4.2.9 0.9 20000	10	12 1.58 0.19 1.29 0.13
	4.2.10 1 19000	14	0	1.66 0.08 1.72 0.08
	4.2.11 1.1 7000	2	0	1.62 0.1 1.63 0.09
	4.4.4 Time scale closing to atmospheric turbulence

Table 4 .

 4 5: Case 4.3 whose cascade step is n = 15.

		Case 4.3
		n η *	∆t	m
		15 13 0.00007 9
	Case	f	t f
	4.3.1	0.5	25001
	4.3.2	0.6	25001
	4.3.3	0.7	13233
	4.3.4	0.8	25001
	4.3.5	0.9	19000
	4.3.6	1	18000
	4.3.7	1.1	11203

  .6. N C 1 of Case 4.3.7 whose forcing is extremely large causing SGC model stops at time 11000 is 2, while it is 3 for Case 4.3.1 whose fluctuations vanish. Because the energy flux of spatial eddies is so large that results in an early stopping, Case 4.3.3 only has two subseries satisfying C 1 ≈ 0.25. These cases demonstrate the absence of a large time scale with UM parameters that are consistent with empirical UM parameters. For the rest, N α and N C 1 in the table suggest that energy flux reflecting properties of atmospheric turbulence has sufficient time scale, and ensemble analysis is presented in the following part. The mean of multifractality index α is 1.63 ± 0.04, which is greater than the empirical value 1.5 as many values α > 2 appears in Case 4.3. C 1 is 0.25 ± 0.07 excluding Case 4.3.1, Case 4.3.3, and Case 4.3.7. Ensemble analysis for Case 4.3 is presented in this part. Case 4.3.1, Case 4.3.3, and Case 4.3.7 have inadequate time scales whose UM parameters agree with empirical UM parameters, according to N C 1 in Table 4.6. Overlapping time scale is not enough for Case 4.3.2, since Case 4.3.2 exhibiting gradually vanishing intermittency has abnormal C 1 after subseries N = 20 and α greater than 2 in the first half of time series. Ensemble analysis of Case 4.3.4-Case 4.3.6 is therefore carried out to demonstrate partial time series reflecting the characteristics of empirical atmospheric turbulence.

	4.5.4 Time scale closing to atmospheric turbulence
	Table 4.6: Frequency of UM parameters of Case 4.3 closing to atmospheric
			turbulence.
			Case 4.3
	n	η *	∆t	m resolution
	15	13	0.00007	9	512
	Case f	t f	N α	N C 1 α	C 1	α e	C 1e
	4.3.1 0.5 25001	16	3	1.63 0.32 1.5 0.1
	4.3.2 0.6 25001	11	7	1.65 0.32 1.47 0.1
	4.3.3 0.7 12000	10	2	1.64 0.13 1.53 0.11
	4.3.4 0.8 25001	15	17 1.66 0.21 1.35 0.11
	4.3.5 0.9 19000	17	10 1.63 0.18 1.33 0.13
	4.3.6 1 19000	11	8	1.66 0.23 1.46 0.13
	4.3.7 1.1 11000	14	2	1.59 0.12 1.56 0.1

Ensemble analysis of Case 4.3.4 is shown in

Table 4 .

 4 7: Frequency of UM parameters of Case 4.1 simulated by semi-implicit Euler method. Mean codimension C 1 is small at the beginning as SGC model is quitting from an ad-hoc initial state and forcing is not yet functioning, followed by stable and rising at the end. The occurrence of abnormal values C 1 > 1 when forcing is small suggests that these subseries are not multifractal fields, but this does not happen when forcing is extremely small. The analysis of these three cascade steps reveals that the multifractality index is significantly smaller than 2. Averaged multifractality index of cascade step 12 except the minimum forcing 0.1 is α = 1.53 ± 0.06 and is 1.66 ± 0.06 for cascade step 14.It is also confirmed by averaged multifractality index 1.63 ± 0.04 for cascade step 15. There are many cases for these three cascade steps whose energy flux is consistent with the empirical atmospheric turbulence α ≈ 1.5, C 1 = 0.25. Time scale for ensemble analysis is limited when forcing is either large or small because of few subseries satisfying C 1 ≈ 0.25 As for forcings of medium value, many subseries satisfying α ≈ 1.5, C 1 ≈ 0.25, which is supported by ensemble analysis, demonstrates that

	Case 4.1 simulated by semi-implicit Euler method.
		n	η *	∆t	m resolution
		12	10	0.0002	6	512
	Case	f	t f	N α	N C 1 α	C 1	α e	C 1e
	4.1.1 0.1 23810	14	5	1.61 0.2 1.35 0.09
	4.1.2 0.2 23810	14	5	1.61 0.2 1.35 0.09
	4.1.3 0.3 23810	10	14 1.62 0.28 1.21 0.13
	4.1.4 0.4 23810	12	8	1.61 0.3 1.32 0.15
	4.1.5 0.5 23810	19	6	1.57 0.28 1.36 0.17
	4.1.6 0.6 18000	17	7	1.58 0.2	1.5 0.17
	4.1.7 0.7 23810	16	11 1.48 0.24 1.21 0.16
	4.1.8 0.8 18000	16	12 1.48 0.23 1.5 0.2
	4.1.9 0.9 23810	17	8	1.52 0.24 1.37 0.15
	4.1.10 1	6000	7	4	1.62 0.12 1.34 0.08
	4.7 Summary of Chapter 4
	Intermittency of forced SGC model with three large cascade steps n = 12, 14, 15,
	which is simulated by classic Euler method, is investigated in this chapter by injecting
	forcing at top layer to achieve quasi-equilibrium state. Forcing is 0.1 ≤ f ≤ 1 for
	cascade step n = 12, and ranges from f = 0.1 to f = 1.1. Cascade step 15 with forcing
	from 0.5 to 1.1 is studied. The proof of intermittency is given out first, followed by
	its multifractality, and then the investigation about computational instability.
	With the exception of one special case, energy flux in the inertial range for these
	three cascade steps yields extraordinarily high values as forcing increases, even bring-
	ing an early ending to SGC model. All these cases confirm the existence of intermit-
	tency. Strong intermittency is clearly observed in cases injecting the medium value

  .5.13) To emphasize that the exponent of energy decay is obtained through a sufficiently high local Reynolds number Re L , the local Reynolds number Re L of SGC model is computed by:

Table 5 .

 5

	1: Cases without power spectrum
	Case n n peak η *	∆t	Tmax
	5.1 11	0	9 0.0002 6000
	5.2 12	0	10 0.0002 6000

Since there is no independent case to perform an ensemble average, statistical descriptions for time slots are carried out with resolution 10.

Table 5 .

 5 3: The exponent of energy decay and integral scale of cases with power spectrum k4 Due to the fact that SGC model preserves the triad interaction generating intermittency, the stationarity of large eddies and scaling laws taking intermittency into consideration are studied by such a model. Theoretical backscatter term is initially tried to be derived from SGC equations, but the semi-local interaction might prevents direct verification. The study thus is carried out with the classic Euler method numerical simulation of SGC model whose two-range spectrum contains power spectrum k 4 and inertial spectrum k -5/3 . Reynolds number of all cases is sufficient to guarantee that SGC model is full developed. It's typical that infrared range decreases and vanishes with time, as evidenced by four cases whose two-range spectrum includes power spectrum k 4 . But phenomena including the conservation of large eddies and the increased exponent of infrared spectrum are observed although with intermittency effects. Compensated spectrum confirms the change of infrared spectrum and inertial spectrum k -5/3 . With the exception of Case 5.3, exponents of energy decay law for case containing infrared spectrum are larger than those for cases lacking infrared spectrum, indicating that backscatter term slows down the energy decay. The buffs happened in the integral scale of these four cases, which suggests the energy is stored at large scales in a discontinuous way. Besides, k 4 at large eddies is observed over a short time scale and with a small wave number length, making it impossible to estimate the significance of backscattering term by Loitsiansky index that slowing down the energy decay. As for numerical simulation, stable semi-implicit Euler method is worth considering for SGC model, as numerical results in this chapter are achieved by the classical Euler method. Due to time restrictions, results of the semi-implicit Euler simulation have not yet been sorted out and included in this chapter. It's also important to look into the energy decay using only k 4 infrared spectrum to enlarge infrared range, revealing the effect brought on by semi-local interaction of SGC model.

	Case	t	α e	γ *
	5.3 100-300 -1.56 0.77
	5.4 100-400 -1.04 0.19
	5.5 100-300 -1.08 0.35
	5.6 100-400 -0.84 0.12
	5.6 Conclusions of Chapter 5 and prospects

More work needs to be done on energy backscatter term to provide both theoretical and numerical simulation evidence. It's challenging to obtain theoretical backscatter term from semi-local triad interactions of SGC model. Besides, it will be interesting to study SGC model by helical modes decompositions, which is typically used to analyze the non-isotropic effects of solid-body rotation on homogeneous turbulence, to provide detailed information about the dynamics of SGC model, such as the manner in which energy is transferred.

  Aiming at filling part of this gap, SGC model which is more close to NS equation than the well known shell model, is simulated in this study, since it preserves the most important energy transfer structure-triad interactions, which is able to generate intermittency and backscatter term. To begin, different numerical simulation approaches are investigated thorough complexity analysis, UM analysis and empirical analysis to determine the effective one. Two of them are the classic Euler method and classical RK4 method. The third one is slaved ABF method, which takes into account the rapid damping of high cases with large forcing, enormous fluctuations emerge and even stop the simulation, resulting in computational instability. Strong fluctuations are obviously observed in cases injecting medium value forcing. Therefore, part of cases supports empirical UM parameters for atmospheric turbulence. Averaged multifractality index for all cases as well as ensemble value is significantly smaller than 2 casting doubt on the log-normal model. Besides, computation instability displayed by these three cascade steps using classic Euler method is investigated by cascade step n = 12 simulated by semi-implicit Euler method, becasue semi-implicit Euler method preserves many symmetries and is more stable than the Euler method, widely utilizing in the motion of rigid body. The outcomes demonstrate that computation instability still exists in the simulation of semi-implicit Euler method. Strong fluctuations are exhibited in the majority of cases n = 12 and UM analysis confirms the conclusion carried out by classic Euler method that multifractality index is significantly less than 2.

	evidence of backscatter term is presented.
	Last but not least, the role of intermittency in backscattering energy ans scaling
	laws are investigated by energy spectrum containing power spectrum k 4 and the
	inertial spectrum, since energy is backscattering for power spectrum k s (s ≤ 4) with
	the help of EDQNM, a closure approach without considering intermittency. Time
	evolution of four cases exhibits the conserved large eddies and increased spectrum
	exponent in infrared range, supporting that phenomenology resulting from power
	spectrum remains the same, although with intermittency effects. It also confirmed
	by observed slower energy decay, when energy decay process is compared to the
	energy decay of case with same inertial range but no infrared spectrum. Besides,
	integrals scale demonstrates the energy is stored at large scales by buffs, no longer
	in a continuous manner. Unfortunately, neither theoretical evidence nor numerical

This figure is from JohnnyOneSpeed

The data is from Hydrology Meteorology and Complexity laboratory of École des Ponts Paris-Tech, which is measured by one of 3D sonic anemometers (manufactured by Thies). It is installed at a wind farm located the southeast of Paris in France.

n-z -g (z) h (z)).And g (z) is the number of supposed dead eddies at (nz)th step which is caused by former steps. Thereforeg (z) = m ′ =z+1 h (m ′ ) • 2 m ′ -z corresponds to the supposed

n-z -g (z) h (z)

The first two decimal places of UM parameters are retained for the sake of labeling concisely in figures, but q s is calculated from the original UM parameters. q s in figures and the one estimated from the labeled UM parameters, thus, differ slightly.

If not specifically noted, theoretical UM parameters remain the same.

Thesis -X. LI; Ecole des Ponts ParisTech -Hydrology Meteorology et Complexity

Acknowledgements

larger than 1.

Based on the β-model, the multifractal is introduced to provide a remarkable accuracy description of irregular fields. Meanwhile, UM framework, which is inspired by the multifractal process and the central limit theorem, is presented as a powerful framework for characterizing the multifractality of weather and climate on a large spatial and temporal scale. A detailed analysis is carried out in the next chapter.

Chapter 3

The Scaling Gyroscope Cascade model and numerical methods

The intermittency of natural turbulence is considered one of the most important phenomena in classical physics, nevertheless, for which neither a thorough theoretical foundation nor phenomenological knowledge is sufficient. Scaling Gyroscope Cascade (SGC) model [START_REF] Chigirinskaya | Cascade of scaling gyroscopes: Lie structure, universal multifractals and self-organized criticality in turbulence[END_REF][START_REF] Chigirinskaya | Scaling gyroscopes cascade: universal multifractal features of 2-d and 3-d turbulence[END_REF][START_REF] Schertzer | Multifractal cascade dynamics and turbulent intermittency[END_REF] was motivated by Arnold's insight [4] that vorticity equation arising from NS equation is analogous to Euler's gyroscope equation. On the basis of this similarity, SGC model preserves the majority of fundamental properties of NS equation, such as the triad interaction and energy conservation. With its capacity to generate the intermittency and quantitative description of natural turbulence, SGC model is regarded as being extremely close to turbulence in nature and boosting the possibility to precisely study intermittency and its role in energy decay. Various numerical techniques are used to provide the velocity time evolution of SGC model,

Algorithm analysis

Algorithm analysis [START_REF] Shaffer | A practical introduction to data structures and algorithm analysis[END_REF] is a quantitative evaluation of resources needed to run a task, which primarily takes execution time and memory usage into account. Due to the fact that it can be carried out in two stages: pre-implementation and postimplementation, priori analysis of these algorithms is performed first, and then the posterior analysis follows.

A priori analysis is an approach that does not rely on specific testing or data, so the pseudo-code of velocity function for previously presented algorithms is provided along with complexity analysis. Space complexity characterizes required storage space S(n) for full execution, including temporary auxiliary space, input space, and output space, as a function of problem size n. It is defined as an asymptotic behavior by following generic formula :

where f (n) is a statement function with regard to the storage space. Space complexity is denoted as S(n) = O(f (n)) merely preserving the highest order term without coefficient. The order of algorithm complexity from simplest to most complicated is

However, the empirical efficiency of algorithms not only depends on the complexity, but also problem size n, as shown in Fig. 3.4.

In order to conduct empirical investigation of algorithms, posterior analysis utilizes data to execute a coded algorithm on target computer and collect practical statistics such as running time and required space. The accuracy of this approach is reliant on testing environment and compiler languages. Indeed, computations for velocity equations distinguishing classical Euler algorithm from semi-implicit Euler technique are basic operations so that posterior analysis for these two approaches is almost identical. Posterior analysis of classical Euler algorithm representing semi-implicit Euler method is used in comparison with the other two algorithms.

Posterior analysis of Case 3.4, the order from smallest to largest is classic Euler method, RK4 method and slaved ABF method.

Max RSS is the maximum physical memory utilised by execution process at any given time. Max RSS of three approaches is extremely close as the maximum memory relies on velocity field that has to be written out and velocities for all eddies at each time step are saved in case 3.1. In terms of efficiency, the classic Euler method, which takes the least amount of time, is an good option, based on posterior analysis. However, time scale used to compare algorithms in Case 3.1 isn't enough to perform properties analysis and time differences will be far more significant when time scale is extended. It also likely that the difference will increase as the cascade step increases, but that is not the focus of this study. Euler methods including classic Euler method and semi-implicit Euler method are effective numerical methods and are employed in the following chapters. In addition, simulation of SGC model is complex and can not be simply improved by numerical methods, as the complexity analysis demonstrates.

It has been demonstrated that, in accordance with the central limit theorem, the sum of random variables from "heavy" tail eventually converges to an α-stable distribution [START_REF] Franzke | The structure of climate variability across scales[END_REF]. UM framework based on the central limit theorem and multiplicative cascade processes is a powerful tool [START_REF] Schertzer | Structures in turbulence and multifractal universality[END_REF] to describe this nonlinear high variability.

The entire time series of energy flux is divided into many sub-series to carry out the temporal evolution of UM parameters as well as corresponding frequency. Ensemble analysis is then conducted on all sub-series. Multifractality index obtained from ensemble analysis of all subsequences is represented as α e , while mean codimension is written as C 1e . α is the average of multifractality index of all sub-series, and C 1 represents the averaged mean codimension. q se is critical moment obtained by ensemble analysis, and q s is the mean of critical moment of all sub-series.

UM analysis of Case 4.1 starts with cases in which a middle range of forcings is imposing and high intermittent energy flux is directly observed. Cases injecting large forcing are analyzed next, followed by those in which small forcing is injected. UM parameters are estimated as the mean value of three scaling moment orders q = 0.9, q = 1.5, and q = 2. Absolute energy flux |Π(k 6 )| has to be normalized to maintain the conservation and perform the UM analysis.

Sample size and bias generated by fixed range estimation

To begin studying intermittency using UM analysis, a sample size has to be selected. The frequency of UM parameters that is close to empirical values is summarised in As statistically uniform estimating range may not be well fitted for every subseries, adjusted range estimation is utilized to estimate UM parameters and critical moments q s . The number of subseries with UM parameters close to α = 1.5 and C 1 = 0.25 is summarized in Table 4 Non-Gaussian distribution of these fluctuations and heavy tails are not given in this section. UM analysis of all cases will be provided in Appendix B. where it remains fully developed, i.e. with a high Reynolds number. The exponent of this decay law can be obtained with the help of only two hypotheses:

• global self-similarity of the energy spectrum and transfer, with the possible exception of the dissipation range [START_REF] Schertzer | Comportements auto-similaires en turbulence homogene et isotrope[END_REF] • stationarity of the energy spectrum of big eddies, i.e. those corresponding to vanishing wave numbers (k → 0) Indeed, these hypotheses yield the two equations that follow:

where F (x) and G(x) are dimensionless functions of the dimensionless wavevector

x, v 2 (t)/2 is the average energy K(t) (Eq.1.5.1); ℓ(t) is a rather arbitrary scale of the energy-containing eddies (e.g. the integral scale); s is the exponent of stationary of big eddies and C s is its prefactor. The latter, via Eq.5.1.1b, yields a first scaling relation between v(t) ∝ t αv and ℓ(t) ∝ t α ℓ :

A second scaling relation is obtained by considering the phenomenological relation of energy dissipation to the velocity ε (Eq.1.4.1):

Combining both scaling relations (Eqs.5.1.2, 5.1.3) yields:

A more rigorous and more general derivation was obtained by Schertzer in 1980 [START_REF] Schertzer | Comportements auto-similaires en turbulence homogene et isotrope[END_REF] by firstly introducing a modified dissipativity [START_REF] Lions | Quelques méthodes de résolution de problemes aux limites non linéaires[END_REF] α D ̸ = 1 to enable to obtain a strict global self-similarity, i.e., including in the dissipation range. The spectral equation is then generalized to:

The scales of dissipation ℓ d,α D and of energy containing eddies ℓ are linked through the modified Reynolds number R α D :

Then it is assumed that the modified dissipativity α D ̸ = 1 can be chosen to obtain a strict global self-similarity. Spectral equation Eq.5.1.5 can be expressed with the help of functions F and G:

Deriving this equation only with respect to time, one faces then the alternative:

• the three functions F, xF ′ , x 2α D F (x) are linearly dependent

• they are linearly independent

The first hypothesis yields solutions of type:

that was rejected by Schertzer [START_REF] Schertzer | Comportements auto-similaires en turbulence homogene et isotrope[END_REF] as non-physically interesting, despite they were investigated in the physical space by [START_REF] Sedov | Methods of similitude and dimensional analysis in mechanics[END_REF].

The second hypothesis yields the following form of spectral equation:

with the series of constant positive prefactors:

The two first prefactors are equivalent to the aforementioned phenomenological relation and the resulting relation between scaling exponents of v and ℓ (Eq.5.1.3).

Together with the assumed stationarity of big eddies (Eq.5.1.2), they yield the expressions of scaling exponents α v and α ℓ with respect to big eddies stationary exponent s. The third prefactor, which is the (modified) Reynolds number, provides the expression of adequate modified dissipativity α D with respect to the exponent s:

In this respect, the case s = 1 is the unique case that does not require any modification of the dissipativity to ensure strict global self-similarity. It confirms a remark made by [START_REF] Lesieur | Self-similar damping of turbulence at large reynolds number[END_REF] on the fact that in this case the analysis of free decay can be performed with

Energy decay law of spectrum injecting power spectrum

Now, free energy decay law of energy spectrum Eq.5.5.11 is investigated and cases are given in Table 5 The existence of intermittency and "heavy" tail of its probability distribution are confirmed by all cases, which suggests UM framework is a powerful tool to describe this nonlinear extreme events. The intensity of fluctuations is related to injected forcing, which is normal considering forcing has to prevent energy from losing to prolong running time. Fluctuations in cases with small forcing vanish over time, however in other

Appendix A

A.1 Legendre Transform

The connection between the scaling moment function K(q) and the codimension function c(γ) through the Legendre transform is explained.

From Eq 2.4.5, the probability of

where the field has many variables γ. The qth moment of ε has

Since Eq.2.5.4 suggests ⟨ε q λ ⟩ can be expressed by the resolution λ, it yields

where ξ = log(λ). Due to Eq A.1.2, the function f (γ) is

Considering the mathematical technique of "steepest descents", the dominant term

. Therefore, the corresponded moment equation is

Eq A.1.5 describing the moment equation K(q) and the codimsion function c(γ) is a Legendre transform. Following same method, another Legendre transform is obtained

For a given q, γ q is defined by the maximum value that increase qγc(γ). q γ is defined as same way.

A.2 Universal Multifractal

The general form of scaling moment function K(q) and codimension function C(γ)

which relates to the UM parameters α, C 1 are explained.

Based on the Taylor expansion, the moment function K(q) is obtained corresponding to a pure log-normal cascade

where A 2 is coefficient. In order to make sure the small scale cascade limit is well behaved, the field has to be normalized by ε → ε/⟨ε⟩, so that K(q) → K(q) -qK [START_REF] Aitchison | the lognormal distribution[END_REF].

And the moment function K(q) is

where A 2 is replaced by C 1 due to Eq 2.6.

However, this is not sufficient to describe the cascade process. A random generator following an "extremely asymmetric" Lévy distribution is added. And it yields

For α = 1, K(q) = C 1 qLn(q) is obtained from L'Hopital's rule for the limit α → 1.

Taking the Legendre transformation of the moment function K(q), the codimen-

When α = 0, it's the monofractal β-model. And α = 2 is the Log-normal model which is the upper limit . The illustration that shows c(γ) could be locally characterized by the singularity C 1 and the local radius of curvature R c (C 1 ) at C 1 . K(q) can be characterized through same method.

Appendix B B.1 Gyroscope Equation

The Euler's equation for a rigid body is known as the gyroscope equation. There are two invariant, the square of the angular momentum and the kinetic energy.

The gyroscope equation is

where M is the angular momentum of a rigid body; Ω is its rotation; ∧ is the vector product. The linear relationship between angular momentum M and rotation Ω is symmetric

where I is the moment of inertial tensor.

Considering Eq. B.1.1 in so(3), the components of the fundamental asymmetric tensor I is the structure constants and the vector product ∧ corresponds to the Lie bracket [ , ]. Now Eq. B.1.1 equals to an adjoint action on Lie algebra a

The velocity field and the wave vector has

The last orthogonality triad is Fig. B.2d

for i = 0, . . . , K[m] -1 do ▷ i is the location at layer m 9:

Compute Set the final velocity list U T := U ▷ The final velocity includes the initial value 21:

for t = 0, . . . , T do ▷ Computing velocities for each location at time t

22:

Compute U (t + ∆t) by Function 1

23:

Store U (t + ∆t) in U T

24:

end for 25:

return U T

26: end function

The Python codes for velocities equations using classic Euler's scheme are given.

First of all, the associated Python code for the velocity at top layer u 0 0 Eq.3.2.3 is:

B.3.3 The slaved Adams Bashforth method

Since the input required by the slaved ABF method is different from that of the one-step methods, the pseudo-code for the slaved ABF method is given.

Algorithm 2

The slaved Adams-Bashforth Method for t = 0, . . . , T do ▷ Computing velocities for each location at time t

22:

Compute U (t + ∆t) by Function 3 ▷ When t = 0, the former velocities

B.3.4 Results

The computing time and maximum memory of Euler method:

Declaration of Academic

Achievement