Marcello Meldi
email: marcello.meldi@ensma.fr

Pierre Sagaut

Marcello Orcid

Meldi

Investigation of anomalous very fast decay regimes in homogeneous isotropic turbulence

Keywords: Isotropic turbulence, EDQNM, turbulence theory

published or not. The documents may come

Investigation of anomalous very fast decay regimes in homogeneous isotropic turbulence

Introduction

The investigation of homogeneous isotropic turbulence (HIT) dynamics is among the most important research topics in the field of turbulence research, since it lies at the heart of most existing turbulence theories (see [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF] for a recent survey). A comprehensive analysis of nonlinear mechanisms which drive the global features of the integral quantities such as turbulent kinetic energy K, the integral scale L and energy dissipation rate ε are still open challenges despite the advances in this research area since the beginning of the twentieth century. Many efforts have been devoted to the case of statistically steady isotropic turbulence, i.e. forced turbulence with a constant turbulent kinetic energy production rate, leading to a significant amount of theoretical, experimental and numerical results. An important point is that such a configuration is associated to an equilibrium state, at which production rate, net energy transfer rate via the non linear kinetic energy cascade and the dissipation rate ε are equal. This state is represented in the spectral space by the Lin equation, which describes the time evolution of the three-dimensional kinetic energy spectrum:

∂E(k, t ) ∂t + 2 ν k 2 E(k, t ) = T (E; k, t ) + F (k, t ) (1)
where ν, E(k, t) and T(E; k, t) denote the kinematic viscosity of the fluid, the energy spectrum and the non-linear energy transfer associated to E(k, t), respectively. F(k, t) represents an arbitrary isotropic forcing term. Statistical equilibrium is reached for

+∞ 0 F (k, t )dk = 2 ν +∞ 0 k 2 E(k, t )dk = ε, dK dt = 0
and, in the case of asymptotic high Reynolds numbers and a forcing term restricted to a finite range of large scales, an inertial range such that ∂E(k,t ) ∂t = T (E; k, t ) = 0 is expected complying with the Kolmogorov solution E(k)Ýε 2/3 k -5/3 . For a wavenumber k located in such a range, the spectral energy flux is constant since (k) = +∞ k T (E; p)dp = ε. One important consequence is that the non-dimensional dissipation parameter C ε = √ 27/8 ε L/K 3/2 [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF][START_REF] Djenidi | On the normalized dissipation parameter c ϵ in decaying turbulence[END_REF] is constant for statistically steady HIT. These last two features are classically associated with turbulence equilibrium. Recent investigation by J.C. Vassilicos [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF] proves that while C ε = const. is actually a true indicator of turbulence equilibrium, an inertial range complying with the relation E(k)Ýε 2/3 k -5/3 is observed in a larger class of flows.

The existence of anomalous fast decay regimes downstream fractal grids and fractal bluff bodies has been recently identified by Vassilicos at Imperial College, followed by a growing number of research groups (see [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF] for a survey).

The origin of this anomalous decay rate might be attributed to each of the two following reasons: [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF] an effect of the peculiar spectrum shape at energetic scales produced by fractal objects or [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF] an effect of long-lasting production mechanisms due to the vanishing shear effects in the wake of such bodies. These effects appear to be associated with a bump of the energy spectrum in the energy production region. This feature is also observed in direct numerical simulation (DNS) solutions with anomalous fast decay reported in [START_REF] Goto | Energy dissipation and flux laws for unsteady turbulence[END_REF][START_REF] Goto | Unsteady turbulence cascades[END_REF] and in data-assimilated initial solutions obtained in [START_REF] Mons | Is isotropic turbulence decay governed by asymptotic behavior of large scales? An EDQNM-based data simulation study[END_REF], showing that such non-classical decay regimes are all tied to a non-classical initial spectrum shape. Additional studies highlighted how these anomalous features are observed as well in classical grid turbulence and that multiscale/fractal designs simply magnify the visibility of these phenomena. This picture is made more complex by the difficulties arising when bridging theoretical analysis with experiments/DNS. In fact while theoretical results are obtained for asymptotic high Reynolds number regimes, experimental and DNS data are related for lower Reynolds numbers at which finite Reynolds number (FRN) effects are not negligible.

The present study aims at providing new results dealing with this question, based on both theoretical analysis and the eddy-damped quasi-normal Markovian (EDQNM) model. The recent version developed by the research group [START_REF] Meldi | On the emergence of non-classical decay regimes in multiscale/fractal generated isotropic turbulence[END_REF] will be used, which includes an ad hoc spectral isotropic forcing F. This term has been rigorously derived from proposals employed for DNS-forced turbulence [START_REF] Mazzi | Fractal-generated turbulence[END_REF]. Using EDQNM, one can ensure that pure isotropy is preserved and that finite Reynolds effects are under control. Detailed features of the initial solution are also prescribed in a fully arbitrary way. Therefore, a clear distinction can be made between the different possible effects. One of the main goals is to assess previous observations dealing with the scaling of C ε and underlying physical mechanisms associated.

The paper is organised as follows. First, some theoretical considerations dealing with the possible increase of the dissipation rate by adding a disturbance to the initial energy spectrum are given in Section 2. Numerical details and EDQNM simulation set-up are addressed in Section 3. EDQNM results are then analysed in Sections 4-6. Several topics are addressed: (1) the influence of the ratio of the mean advection time versus the turbulent eddy turnover time, (2) the effect of a continuous fractal/multiscale forcing over a range of large scales and then (3) the effect of a forcing restricted to a discrete set of scales within the same range. Concluding remarks are given in Section 7. Details of the EDQNM model used in the present study are given in the Appendix for the sake of clarity.

Decay acceleration via local energy injection: a look at Lin equation

We first analyse the influence on the dissipation rate of the local injection of kinetic energy into a classical solution. First, it is reminded that the exact evolution equations for kinetic energy and dissipation obtained by integrating Lin's Equation [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF] are

dK dt = -2 ν +∞ 0 k 2 E(k, t )dk =ε + +∞ 0 T (E; k, t )dk =0 + +∞ 0 F (k, t )dk =P (2)
dε dt = -4 ν 2 +∞ 0 k 4 E(k, t )dk + 2ν +∞ 0 k 2 T (E; k, t )dk + 2ν +∞ 0 k 2 F (k, t )dk =P (3)
Now let us consider a base solution E 0 (k, t) obeying the Lin equation to which a disturbance E(k, t) is added. Assuming that (E 0 + E) is also a solution of the Lin equation, one obtains by subtraction

∂ E(k, t ) ∂t + 2 ν k 2 E(k, t ) = T (E 0 , E; k, t ) + F (k, t ) (4)
where

T(E 0 , E; k, t) = T(E 0 + E; k, t) -T(E 0 ; k, t)
is the disturbance in the transfer term and F(k, t) is the possible change in the forcing term in the case it is solutionsensitive. It is important noting that Equation ( 4) is general and does not come from a linearisation process. Therefore, high-energy disturbances can be considered. Associated exact equations for disturbances in the kinetic energy

K = +∞ 0 E(k)dk and dissipa- tion rate ε = 2ν +∞ 0 k 2 E(k)dk are d K dt = -2 ν +∞ 0 k 2 E(k, t )dk = ε + +∞ 0 T (E 0 , E; k, t )dk =0 + +∞ 0 F (k, t )dk = P (5) d ε dt = -4 ν 2 +∞ 0 k 4 E(k, t )dk + 2ν +∞ 0 k 2 T (E 0 , E; k, t )dk + 2ν +∞ 0 k 2 F (k, t )dk = P (6)
It is worth noting that T(E 0 , E; k, t) is a true non-linear energy transfer term, since +∞ 0

T (E 0 , E; k, t )dk = 0. The emphasis is put here on disturbances that are local in the wavenumber space, i.e. they are restricted to a narrow band in the Fourier space. Two cases can be considered: (i) a disturbance localised at energetic scales near the peak of E 0 (k, t) or (ii) a disturbance localised at very small scales, in the neighbourhood of a wavenumber k δ such that k δ L 1. In this last case, L is the integral scale associated with the base solution E 0 . Fractal grid experiments and previous EDQNM results dealing with fractal forcing do all correspond to case (i), which is justified by the consideration that the same turbulence production mechanisms are driving both the base flow and the disturbance.

Let us now neglect forcing effects during the decay phase, so that the analysis will be focused on the free decay of an initially perturbed base flow. This corresponds to case in which the fractal forcing vanishes quickly and its impact is limited to the initial flow configuration. Equation ( 5) simplifies as

d K dt = -ε (7)
while the situation is more complex for the dissipation rate disturbance. In the case of a large scale forcing, it is reasonable to neglect the first term in the right-hand side of (6) because of the k 4 weight in the integrand, leading to

d ε dt = 2ν +∞ 0 k 2 T (E 0 , E; k, t )dk (8)
Since T(E 0 , E; k, t) is a true transfer term associated to a net direct energy cascade (this point is a posteriori assessed looking at EDQNM solution), the k 2 weight leads to the fact that the right-hand side of ( 8) is positive at least for the very first stages of decay, leading to a growth of ε and therefore an increase in the total dissipation ε 0 + ε.

The case of small-scale forcing is more complex, since no term can be a priori neglected. The destruction of dissipation rate -4 ν 2 +∞ 0 k 4 E(k, t )dk is strictly negative, and the evolution of ε is governed by the balance between a strictly negative term and a positive one, whose exact amplitudes are solution-dependent. One could a priori consider cases with positive (forcing in the far inertial range of E 0 (k) with dominant transfer effects) or negative (forcing in the dissipative region of E 0 (k) with negligible non-linear transfers) right-hand side.

Now restricting the analysis to the large-scale forcing case, the transfer disturbance T(E 0 , E; k, t) can be further analysed. To this end, one must use a closed expression for the transfer term. The EDQNM model used here for the triad k, p, q so that p + q = k (see Appendix) yields

T (E(k), E(p), E(q); k, t ) = p+q=k kpq g E(q, t ) E(p, t ) pk 2 -E(k, t ) p 3 dp dq p q (9)
where the characteristic time kpq is also a function of E(k), E(p) and E(q) and g is a geometric flow-independent factor. Using that expression, one obtains for the transfer disturbance:

T (E 0 , E; k, t ) = T ( E(k), E 0 (p), E 0 (q); k, t ) = T 1,1 + T (E 0 (k), E(p), E 0 (q); k, t ) = T 1,2 + T (E 0 (k), E 0 (p), E(q); k, t ) = T 1,3 + T ( E(k), E(p), E 0 (q); k, t ) = T 2,1 + T ( E 0 (k), E 0 (p), E(q); k, t ) = T 2,2 + T (E 0 (k), E(p), E 0 (q); k, t ) = T 2,3 + T ( E(k), E(p), E(q); k, t ) = T 3 ( 10
)
Term T 3 is related to transfer disturbance self-induced by the disturbances, while all other terms are associated to cross-interactions between the base flow and the disturbance. This expression is fully general and must be refined for the case of a localised disturbance around kL = 1, i.e. E is zero outside this narrow spectral region. Since we are interested in describing the evolution of K, the emphasis is put on the transfer disturbance for wavenumbers k close to k L ∼ 1/L, and therefore E(k) ࣔ 0. Expression [START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of grid-generated turbulence[END_REF] shows that T 1, 3 , T 2, 2 , T 2, 3 and T 3 vanish if q is not close to k L , i.e. if E(q) = 0. In the non-vanishing case, p is also related to large scales, since p + q = k. Therefore, these terms correspond to transfers induced by interactions between large scales. T 3 represents local triadic interactions among modes located in the forced spectral region. Since the forcing is restricted to a narrow region, one has for non-vanishing values E(k) ∼ E(p) ∼ E(q) ∼ L K, along with dominant cross-interactions such that E(k) ∼ LK. In the same way, dominant interactions are associated to kpq ∼ L/ √ K 0 + K. Arguably, this energy transfer is magnified with an increasing size of the forcing range.

Using these evaluations, and considering a small disturbance case, i.e. now assuming that K is small with respect to K 0 , dimensional analysis yields

T 1,3 ∝ L √ K 0 + K (LK 0 )(L K)L -3 ∼ K 0 K T 1,i (i = 1, 2) ∝ L √ K 0 + K (LK 0 )(LK 0 )L -3 ∼ K 3/2 0 T 2,i (i = 1, 2, 3) ∝ L √ K 0 + K (LK 0 )(L K)L -3 ∼ K 0 K T 3 ∝ L √ K 0 + K (L K)(L K)L -3 ∼ K 2 / K 0
The dominant terms are T 1, 1 and T 1, 2 , showing that the growth of ε and therefore the acceleration of the decay of K is mainly governed by cross-interactions between the base flow and the disturbance field associated to local energy transfers.

In the opposite case of large disturbances, i.e. K K 0 , which is the most relevant to analyse existing data that exhibit an anomalous fast decay, one has

T 1,3 ∝ L √ K 0 + K (LK 0 )(L K)L -3 ∼ √ KK 0 T 1,i (i = 1, 2) ∝ L √ K 0 + K (LK 0 )(LK)L -3 ∼ K 2 0 / √ K T 2,i (i = 1, 2, 3) ∝ L √ K 0 + K (LK 0 )(L K)L -3 ∼ √ KK 0 T 3 ∝ L √ K 0 + K (L K)(L K)L -3 ∼ K 3/2
in which the dominant term is T 3 . Dimensional analysis then yields for short-time evolution during which K K 0

d dt K ∝ T 3 /L ∼ K 3/2 /L (11)
In the same way, one obtains

d dt ε ∝ 2ν T 3 /L 3 ∼ 2ν K 3/2 /L 3 (12)

Comte-Bellot/Corrsin type theoretical analysis

In order to understand how previously discussed modal interaction affects HIT decay, we now develop a Comte-Bellot-Corrsin-type simplified model [START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of grid-generated turbulence[END_REF] based on dimensional analysis and simplified spectrum shapes. Let us consider a global energy spectrum E = E 0 + E defined as

E 0 (k) = ⎧ ⎨ ⎩ Ak σ , k < k L Ak σ +5/3 L k -5/3 , k L ≤ k ≤ k η 0, k > k η (13) E(k) = D F E 0 (k L ), k f ≤ k ≤ k F 0, elsewhere (14)
The spectrum E 0 in Equation ( 13) is a classical two-range asymptotic spectrum exhibiting an inertial range for wavenumbers larger than k L = 1/L. A feature that is here introduced is that the spectrum is filtered for wavenumbers larger than

k η = 1/η, where η = L Re -3/4 T = (20/3) 3/4 L Re -3/2 λ
is the Kolmogorov scale. The integral of this spectrum provides the base value for the turbulent kinetic energy K 0 = k η 0 E 0 dk. On the other hand, the spectrum E is considered as a perturbation to the base spectrum E 0 which can be physically associated with production mechanisms. In this case, this perturbation is constant and proportional to the value of the peak of the base spectrum E 0 (k L ) through a strictly positive constant D F . The perturbation is non-zero in a limited range of the spectral space, which is included between the two wavenumbers k f and k F . The integral of the spectrum E provides the turbulent kinetic energy perturbation

K = k F k f E(k)
dk and, as a consequence, the integral of the total energy spectrum E is equal to

+∞ 0 E(k)dk = K = K 0 + K.
A representation of the energy spectrum is provided in Figure 1.

Two different families of configurations can be identified, depending on the values provided for D F and Re λ . The resulting initial features of HIT decay can exhibit very different behaviours. The reason why can be understood via the analysis of the spectrum k 2 E(k), which is proportional to the total dissipation spectrum 2νk 2 E(k). A visualisation of two asymptotic cases representing the features of the two groups of configuration is provided in Figure 2. In Figure 2(a), a value of the parameters D F = 3 and Re λ = 10 6 is chosen. In this scenario, the energy perturbation associated with the spectrum E does not significantly modify the energy dissipation spectrum and the peak (i.e. maximum dissipation) is observed at the very small scales. Results for this configuration are expected to be very close to classical asymptotic analyses where D F = 0 and Re λ → +Ý. On the other hand, the configuration in Figure 2(b) for D F = 10 2 and Re λ = 10 3 represents a much more plausible configuration which can be analysed in experiments, because of the value of Re λ and the clear energy bump at the large scales that is usually observed in grid turbulence [START_REF] Melina | Vortex shedding effects in grid-generated turbulence[END_REF]. In this case, it is clear that the energy dissipation is most active at the large scales. This relevant difference sets the stage for important differences in the initial evolution of the turbulent kinetic energy [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF].

Let us neglect active production effects so that P = 0 and let us consider a suitable time scale t F ࣔ 0 which is in magnitude very similar to the base turnover time t 0 = K 0 /ε 0 . In addition, an evolution law for the physical quantities must be provided. Consistently with most of the theoretical works in the literature, we will consider instantaneous power law evolution of the terms in Equation ( 2) so that

∂K(t )/∂t = n K /t F (K 0 + K)(1 + t/t F ) n K -1 and equivalently ε(t ) = (ε 0 + ε)(1 + t/t F ) n ε , with n ε = n K -1.
This implies that an instantaneous estimation of the power law coefficient n K can be derived as

n K = -t F ε 0 + ε K 0 + K = -t F ε 0 K 0 -t F εK 0 -ε 0 K K 0 (K 0 + K) ( 15
)
where

-t F ε 0 / K 0 = n * K = -2(σ + 1)/(σ + 3
) is connected with the features of the initial spectrum and in particular with the large scale slope coefficient σ . Here, n * K is the power law exponent derived via Comte-Bellot/Corrsin analysis for K 0 and is negative in value. In the case D F = 0, Re λ → +Ý the second right hand term of Equation ( 15) is zero, and the Comte-Bellot asymptotic analysis is recovered. On the other hand, the second right hand term of Equation ( 15) becomes progressively more important for increasing values of D F . Further manipulation yields

n K = n * K 1 + ε/ε 0 1 + K/K 0 (16)
Thus, the decay exponent n K can exhibit significant variations when compared with n * K . If k F if significantly larger than k f , the following approximation holds:

ε = k F k f k 2 D F E 0 (k L )dk ≈ 1 3 D F E 0 (k L )k 3 F (17) K = k F k f D F E 0 (k L )dk ≈ D F E 0 (k L )k F ( 18
)
this approximation indicates that, if k F is sufficiently larger, the augmented dissipation rate due to the spectral imbalance is initially much stronger than the energy injection. This aspect results in a much faster turbulent kinetic energy dissipation. However, because of triadic interactions previously discussed, this augmented dissipation can vanish rapidly. Eventually, a classical dynamic behaviour exhibiting a peak for k 2 E at the small scales is obtained, if the Reynolds number initially imposed is sufficiently high.

EDQNM: set-up and numerical details

The EDQNM model is based on the numerical resolution of the Lin equation in the spectral space, which is obtained via the application of two fundamental hypotheses: (1) the deviation of the fluctuating velocity pdf is accounted via a linear eddy-damping term and (2) the characteristic time associated with the eddy-damping rate is neglected with respect to the characteristic evolution time. Elements describing the closures for the non-linear energy transfer term T and the multiscale forcing F are provided in the Appendix and extended discussion is elaborated in [START_REF] Meldi | On the emergence of non-classical decay regimes in multiscale/fractal generated isotropic turbulence[END_REF][START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF][START_REF] Meldi | The signature of initial production mechanisms in isotropic turbulence decay[END_REF]. In particular, the multiscale forcing term F is governed by the three parameters 2fγ R , β and α which control the intensity of the forcing, its spectral distribution and its time evolution, respectively. Each EDQNM calculation goes through three different global steps in order to mimic the generation of turbulence in grid wind-tunnel experiments:

(1) Initialisation. The initial energy spectrum is imposed using the functional form proposed by Meyers and Meneveau [START_REF] Meyers | A functional form for the energy spectrum parametrizing bottleneck and intermittency effects[END_REF]:

E(k) = C k ε 2/3 k -5/3 f L (kL) f η (kη) (19)
with:

f L (kL) = kL [(kL) 3/2 + c L )] 2/3 5/3+σ , f η (kη) = exp(-ξ ([(kη) 4 + c 4 η ] 1/4 -c η )) (20)
where c η = 0.4, ξ = 5.3 and c L is chosen to obtain L(0) = 1. In addition, the initial Reynolds number is set equal to Re λ = 2000 to prevent FRN effects. This implies that an approximation for the initial Kolmogorov length scale is η(0) = L(0) Re -3/4 T = (20/3) 3/4 L(0) Re -3/2 λ . The largest wavenumber of the numerical mesh, k max , is selected so that k max η(0) = 10. On the other hand, the smallest wavenumber k min is set to comply with the relation k min L(0) = 10 -9 . This choice allows to exclude large scales saturation effects, which eventually affect the numerical prediction when k min L(0) > 10 -2 [START_REF] Meldi | Turbulence in a box: quantification of large-scale resolution effects in isotropic turbulence free decay[END_REF]. The resulting mesh elements are distributed over 14 decades in the spectral space. Their distribution follows a geometrical progression of ratio 1.12 and the spectral space investigated is discretised in 287 elements.

(2) Steady forced state. The EDQNM model is advanced in time-solving Lin's Equation [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF], with the addition of the fractal forcing term (A6). This phase corresponds to evolution times t < τ, for which the forcing F is assumed to be steady. τ is a time virtual origin which has to be set. Preliminary analysis indicated a best value of τ = 10 t 0 (0). Because of the choice of the value of the forcing constant 2 f γ R discussed in the Appendix, the solution evolves towards a new, more energetic state, until an equilibrium between the production term and the energy dissipation rate is reached. The equilibrium is reached at t 0.8τ , which implies that a completely steady turbulent state is observed for 0.8τ < t < τ . This new flow condition is characterised by the quantities L(τ ) 0.9 L(0), Re λ (τ ) 4000 and η(τ ) 2η(0) for a forcing range of one decade (i.e. r f = k F /k f = 10). (3) Time-evolving forced state. For times larger than τ , the forcing F exhibits a time evolution, which is governed by the parameter α. This regime mimics the evolution of the wake of fractal grids downstream the location at which the peak of turbulent kinetic energy is observed in wind-tunnel experiments (the wake interaction length in the parlance of J.C. Vassilicos) [START_REF] Staicu | Turbulent wakes of fractal objects[END_REF][START_REF] Seoud | Dissipation and decay of fractal-generated turbulence[END_REF][START_REF] Valente | The decay of turbulence generated by a class of multiscale grids[END_REF][START_REF] Valente | Universal dissipation scaling for nonequilibrium turbulence[END_REF]. Depending on the value of this parameter, various regimes can emerge as investigated by Meldi et al. [START_REF] Meldi | On the emergence of non-classical decay regimes in multiscale/fractal generated isotropic turbulence[END_REF]. Two specific features of these regimes are (1) both classical algebraic and anomalous fast evolution of the physical quantities can be observed and (2) the turbulence coefficient

C ε = √ 27/8 ε L/K 3/2 can exhibit time evolution.
A qualitative visualisation of the energy spectra E(k, t) obtained via EDQNM is shown in Figure 3. It is important to note that the fractal/multiscale forcing results in a non-classical shape of the peak kinetic energy spectrum, with a significant bump near the peak. Here, a classical initial spectrum shape is a shape with a smooth peak, as found in existing models for E(k), e.g. Pope's model and the Meneveau-Meyers model. Such an unusual spectrum can be adequately modelled using the framework developed in Section 2.1, the disturbance E(k, t) being associated with the anomalous shape of the peak of the computed energy spectrum. Meldi et al. [START_REF] Meldi | On the emergence of non-classical decay regimes in multiscale/fractal generated isotropic turbulence[END_REF] provided a comprehensive investigation of long time evolution of several regimes for α ࢠ [10 -5 , 10 2 ], while Meldi [START_REF] Meldi | The signature of initial production mechanisms in isotropic turbulence decay[END_REF] focused the attention of the features of early time development for very low values of α. In the present manuscript, HIT decay is investigated in the first characteristic times of evolution in the range of parameters α ࢠ [10 -1 , 10 1 ] and [START_REF] Djenidi | On the normalized dissipation parameter c ϵ in decaying turbulence[END_REF][START_REF] Melina | Vortex shedding effects in grid-generated turbulence[END_REF] where r f is the total forced range in the spectral space. Such values of α correspond to a fast decay of the forcing term after τ for the values of α investigated, as shown in Figure 4. The choice for this range of investigation stems from observation of several experimental work reported in the literature dealing with grid turbulence produced by unconventional grids [START_REF] Valente | The decay of turbulence generated by a class of multiscale grids[END_REF][START_REF] Valente | Universal dissipation scaling for nonequilibrium turbulence[END_REF][START_REF] Krogstad | Freely-decaying, homogeneous turbulence generated by multiscale grids[END_REF][START_REF] Krogstad | Near-field investigation of turbulence produced by multi-scale grids[END_REF][START_REF] Thormann | Decay of homogeneous, nearly isotropic turbulence behind active fractal grids[END_REF][START_REF] Hearst | Decay of turbulence generated by a square-fractal-element grid[END_REF][START_REF] Hearst | Velocity derivative skewness in fractal-generated, non-equilibrium grid turbulence[END_REF][START_REF] Hearst | Effects of multi-scale and regular grid geometries on decaying turbulence[END_REF], and it may adequately be analysed using the theoretical framework previously discussed, in which residual effects of the forcing term F were neglected. The observation region of these experiments spans moderate decrease in Re λ , usually by a factor about two. One of the main advantages of the present analysis is that FRN effects are completely excluded. In fact, because of the very high Reynolds number at the virtual origin τ , the condition Re λ > 250 in the observation window is satisfied for all EDQNM simulations.

k 2 F = k -2 f = r f ∈ [
One must make a distinction between the analysis here proposed and the classical initialisation of numerical simulations, where the bulk of the turbulent kinetic energy is initially distributed over the infrared modes. In the present work, the exotic spectrum obtained at the beginning of the time-evolving state exhibits fundamental similarities with spectra obtained in grid experiments [START_REF] Melina | Vortex shedding effects in grid-generated turbulence[END_REF]. Thus, the initial phase of the time evolution is not here a numerical transient, but a fast evolution decay regime starting from a realistic configuration.

Advection time vs. eddy turnover time in grid turbulence

Among the various elements of discussion in recent works dealing with grid turbulence, the most common conclusions drawn are that (1) classical scaling laws [START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of grid-generated turbulence[END_REF][START_REF] Davidson | The minimum energy decay rate in quasi-isotropic grid turbulence[END_REF][START_REF] Meldi | On non-self-similar regimes in homogeneous isotropic turbulence decay[END_REF] are eventually observed sufficiently far downstream and that (2) a higher degree of complexity in the shape of the grid usually results in a faster decay of turbulent kinetic energy in the wake. Hearts and Lavoie [START_REF] Hearst | Velocity derivative skewness in fractal-generated, non-equilibrium grid turbulence[END_REF] suggest that this latter observation could be connected with lingering effect of energy production mechanisms. However, one aspect that would deserve more investigation is the connection between the initial turbulent intensity and the average advection time/characteristic time of the mean shear induced by the bar. Equivalently, this aspect represents the relation between advection time and turbulence turnover time. One of the main goals of the present work is to investigate this aspect via EDQNM model which, thanks to the friendly features of the model, will be isolated from FRN effects, confinement effects, anisotropy effects and intermittency effects [START_REF] Meldi | Turbulence in a box: quantification of large-scale resolution effects in isotropic turbulence free decay[END_REF]. Also, the analysis of this feature is a simple parametric analysis in this case. As a matter of fact, the parameter α is determined by the ratio of the initial eddy turnover time and the advection time, αÝt 0 /t A = C M t 0 /t A . Thus, introducing the grid Reynolds number Re M associated with the production mechanisms, the following relation can be derived:

α = C M t 0 t A = 3 L(τ ) 2 C ε (τ ) u (τ ) U ∞ d → α = C M 3 2 C ε L 2 τ d 2 Re M Re L (τ ) (21)
where U Ý and d are the characteristic velocity/length of the production mechanism. Thus, a direct link between EDQNM calculations and experiments can be established, once a value for the constant C M is provided. A qualitative comparison is shown in Figure 5 for the case α = 0.5, r f = 3, where a very good agreement with experimental data for bi-planar grids is obtained via tuning of the constant C M . The constant has been set so that Re M /Re L (t = τ + 10 -3 t 0 ) = 35. The quantity investigated is the non-dimensional parameter C ε . This coefficient has been object of recent studies [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF][START_REF] Goto | Energy dissipation and flux laws for unsteady turbulence[END_REF][START_REF] Goto | Unsteady turbulence cascades[END_REF][START_REF] Valente | Universal dissipation scaling for nonequilibrium turbulence[END_REF][START_REF] Hearst | Velocity derivative skewness in fractal-generated, non-equilibrium grid turbulence[END_REF][START_REF] Bos | Spectral imbalance and the normalized dissipation rate of turbulence[END_REF][START_REF] Bos | Dissipation in unsteady turbulence[END_REF] because it is a clear indicator of the turbulence equilibrium state. The results indicate that for higher values of the parameter α the curves tend to slides towards higher initial values of Re M /Re L . Assuming that the turbulence turnover time is the same, this implies that higher values of α are associated with a reduction of the characteristic advection time t A . This aspect is tied with different slopes for the initial evolution of C ε , which justifies the difficulties in the comparison of different experimental set-ups by the research teams working on this topic. In addition, this behaviour for the coefficient C ε has been observed for every value of the parameter α investigated and, more in general, for every forced EDQNM simulation. Actually, the only EDQNM calculation that does not exhibit this behaviour is the reference non-forced one, for which the decay phase starts from the Meyers-Meneveau spectrum of Equations ( 19) and [START_REF] Krogstad | Freely-decaying, homogeneous turbulence generated by multiscale grids[END_REF]. This observation suggests that the bump in the spectrum is probably the governing aspect of the anomalous decay observed. This feature is clearly identified as well in experiments for classical and fractal grids [START_REF] Melina | Vortex shedding effects in grid-generated turbulence[END_REF] and DNS [START_REF] Laizet | The spatial origin of -5/3 spectra in grid-generated turbulence[END_REF]. These results indicate that re-normalisation of experimental data with similar values of the parameter α (or equivalently, ratio t A /t 0 ) could highlight similarities in the physical evolution of turbulent features.

In addition, the parameter α could be used to control the energy distribution in the wake in industrial applications dealing with multiscale turbulent energy production [START_REF] Gehlert | Effect of fractal endplates on the wingtip vortex[END_REF].

Early decay stages: continuous forcing

EDQNM results for continuous forcing, i.e. forcing over all wavenumbers in the range

k ∈ [ r -1 f , √ r f ] = [k f , k F ]
, are now investigated. This is an asymptotic approximation of realistic configurations, where a discrete set of modes are associated with production mechanisms.

Results for α = 0.5 along with r f = 3 and r f = 10 are first investigated. In order to perform a relevant comparison, an ad hoc value of the constant C M is calculated for the case r f = 10, so that the condition Re M /Re L (t = τ + 10 -3 t 0 ) = 35 is valid for both calculations. Energy spectra E(k, t) and the non-linear transfer energy budgets k T(E; k, t) are reported in Figure 6 at three different times corresponding to Re M /Re L = 35, 45, 60. Here, the Reynolds number Re L decreases approximately by a factor about two in the investigation range, similarly with the observation windows of grid turbulence experiments. In addition, these values encompass the transient in which the coefficient C ε exhibits non-classical time evolution. If we consider the physical time normalised over t 0 , the sampled spectra are taken for t -τ = 10 -3 , 0.5, 2.4 for r f = 10 and t -τ = 10 -3 , 1.1, 2.5 for r f = 3. It is possible to observe that while the last sample is performed approximately at the same time, the second one occurs much earlier for r f = 10. This effect is due to a much faster energy dissipation for the case r f = 10, which is associated with the energetic local interaction that occurs in the forced range. In both cases one observes that the energy spectrum exhibits a non-classical shape characterised by a range of very energetic modes at the peak when compared to classical spectrum shape obtained in the unforced case shown for the sake of comparison. This bump is associated with the disturbance E(k, t) investigated in the previous theoretical analysis and it is progressively smoothed and dissipated with the time evolution of turbulence. This unusual form of the peak shape is also seen to be associated to a dramatic increase of T(E; k, t) for modes located in the energy peak, showing that the transfer disturbance T(E, E; k, t) may be very strong, leading to a transient growth of ε. Because of the very nature of the EDQNM model, we can thus exclude that this result, which is also reported in grid turbulence experiments [START_REF] Hearst | Velocity derivative skewness in fractal-generated, non-equilibrium grid turbulence[END_REF][START_REF] Hearst | Effects of multi-scale and regular grid geometries on decaying turbulence[END_REF], is associated with anisotropy or intermittency phenomena.

As previously mentioned, the configurations for r f = 10 appear to decay faster than the simulations for r f = 3, because of the stronger interactions between excited modes, still in agreement with the analysis of triadic interactions presented above. This result complies with the statement that a non-zero T(E, E; k, t) yields a non-zero initially growing ε and therefore an increase of the total dissipation rate. This effect is well observed in the local prediction of the power law coefficients of the turbulent kinetic energy K and the energy dissipation rate ε, which are displayed in Figure 7. The values of the power law exponents here reported are much higher than the ones obtained using classical decay theory and the ones reported in fractal grid experiments, which are usually such that n K > -3. The reason for that difference may be attributed to two distinct aspects: (1) the time sampling in EDQNM is extremely fine, allowing for capturing more features than experimental sampling and (2) the forcing here investigated is continuous over the spectral mesh, which results in 20 and 10 excited modes for r f = 10 and r f = 3, respectively. In wind-tunnel grid turbulence usually 3-4 discrete modes are excited, so that the very strong energy transfer that is observed for this configurations is not as magnified as for the continuous forcing applied in EDQNM simulations, even if the parameter r f is the same. This aspect will be further discussed in Section 6. Similar conclusions can be drawn from the time evolution of the power law exponents, which is displayed in Figure 8. It appears that the higher the value of the parameters α and r f , the faster the turbulent kinetic energy decays before approaching the asymptotic limit of classical free decay predicted by Comte-Bellot and Corrsin theory. Again, these findings are consistent with observations by Hearst and Lavoie [START_REF] Hearst | Velocity derivative skewness in fractal-generated, non-equilibrium grid turbulence[END_REF], which observed that fractal grid turbulence exhibits a much lower turbulent kinetic energy downstream, when compared to standard grid turbulence. This can be justified if one assumes that a higher fraction of the total turbulent kinetic energy is associated with the peak disturbance E(k, t) in fractal grid turbulence when compared with classical grid turbulence. The fast decay of K related to ε yields the conclusion. However, a completely different time evolution is observed in Figure 8 for the classical EDQNM calculation. While the final asymptotic limit is the same, the present results show how strong the impact of production mechanisms is in the initial phases of turbulence decay. Arguably, the integration of realistic features of the energy spectrum is a necessary step to capture the emergence of anomalous regimes via numerical simulation, at least for turbulence closures. In addition, the extent of this initial range of decay is driven by t 0 and t A . This implies that, depending of the initial set-up, this region can easily cover a significant part if not the whole test region in grid turbulence experiments.

0 1 2 3 -16 -12 -8 -4 0 0 1 2 -24 -20 -16 -12 -8 -4 0 (a)
The analysis of Figures 6 and9 indicates that the anomalous behaviour of C ε is associated with the behaviour of E(k) and T(E; k) at the upper end of the inertial range. As a matter of fact, intense triadic interactions are observed for Re M /Re L = 35, which are responsible for a faster evolution of global physical quantities such as turbulent kinetic energy and dissipation. Clearly, the larger the parameter r f , the more extreme the non-linear energy transfer.

It appears that the spectrum and the non-linear energy transfer become almost smooth for Re M /Re L = 60, when the time evolution of C ε becomes significantly less intense. This is well understood comparing the forced simulation data with a classical simulation without forcing initialised using a functional form as in Equations ( 19) and [START_REF] Krogstad | Freely-decaying, homogeneous turbulence generated by multiscale grids[END_REF]. In this case, no mode is excited and the initial spectrum is smooth. The spectra, which are sampled for t = 0.3, 1, 2.5 (see Figure 6(e,f)), indicate that in this case no peculiar feature is observed. An important point that should be stressed is that the results in Figure 6 exhibit a sensitivity to two different aspects of the production mechanisms: (1) production is responsible for the shape of the energy spectrum at the virtual origin and (2) lingering production effects influence turbulence evolution after the virtual origin.

Finally, the time evolution of the coefficient C ε is reported in Figure 9. Early times evolution is shown in Figure 9(a). It appears that for the classical EDQNM simulation without forcing the condition C ε const. emerges at very early evolution times when compared to forced EDQNM simulations. The convergence toward an equilibrium state in classical freely decaying turbulence was studied in [START_REF] Bos | Spectral imbalance and the normalized dissipation rate of turbulence[END_REF], in which EDQNM was observed to very accurately match theoretical analysis. We could argue that this observation is associated with the precision of the initial energy spectrum imposed, where minimal variations are observed when compared with EDQNM calculated spectra. On the other hand, the higher the value of the parameter α and r f , the faster the evolution of C ε towards a constant value. This observation is again associated with the intense non-linear interactions due to the effect of production mechanisms. In Figure 9 shown. Here, the coefficient C ε goes through two different phases of time evolution. The first, observed for higher Re λ values, is associated with previous discussion about the initial shape of the energy spectrum. The second is driven by interactions between the lasting effects of production mechanisms and turbulent dissipation, as investigated by Meldi et al. [START_REF] Meldi | On the emergence of non-classical decay regimes in multiscale/fractal generated isotropic turbulence[END_REF]. It is interesting that both regimes exhibit similar Re a λ , a ≈ -1 behaviour and they are consistent with DNS results by Goto and Vassilicos [START_REF] Goto | Energy dissipation and flux laws for unsteady turbulence[END_REF], also in good agreement with the value -15/14 predicted by Bos and Rubinstein [START_REF] Bos | Dissipation in unsteady turbulence[END_REF]. This observation was obtained for different initial values of the parameter Re λ (0), which allows to exclude any FRN effect.

At 10 for three different calculations. The energy flux has been calculated for the three wavenumbers k = L(τ ) -1 , k = λ(τ ) -1 and k = (L(τ ) • λ(τ )) -0.5 . In addition, a fourth wavenumber has been investigated for the forced simulations. It has been identified as the first non-forced wavenumber in the inertial range, so that k = r 0.5 f . The results indicate that, for t = τ , the energy flux for the classical simulation is equal to zero, while a perfect balance in the inertial range is obtained for the forced simulation. This observation is justified if one considers that, for t ࢠ [0, τ ], a steady-forced turbulence state is artificially produced. While the normalised energy flux seems approximately the same for k = λ(τ ) -1 and k = (L(τ ) • λ(τ )) -0.5 , it appears that a return towards equilibrium is accelerated increasing the range of the initial forcing. This observation is linked with the higher energy dissipation at early stages for higher r f values, which implies a faster decrease of Re λ and progressively lower spectral unbalance [START_REF] Bos | Spectral imbalance and the normalized dissipation rate of turbulence[END_REF]. Similar considerations can be performed for k = L(τ ) -1 and, in particular, for k = r 0.5 f . A very important energy transfer is observed for this wavenumber in early stages of decay, which increases with the extension of the forced range r f .

Case of a discrete forcing

The forcing term in the EDQNM model is now modified to achieve a more realistic representation of the physical phenomena acting in grid turbulence experiments. In that case, production effects are not continuous in the excited spectral width, but they are associated with wavenumbers linked to the thickness of grid bars. The EDQNM cannot represent physical instantaneous features of the velocity field. Thus, the objective is to obtain realistic features of the spectrum [START_REF] Laizet | The spatial origin of -5/3 spectra in grid-generated turbulence[END_REF] at the beginning of phase 3 (time-evolving forced state) of the calculation as described in Section 3. To this aim, the EDQNM forcing term is modified in order to excite four modes only in the range r f = 3 and r f = 10, respectively. The discrete-forced range will be referred to as r fD . The wavenumbers are distributed following a geometric progression of constant ratio, which is reminiscent of the process of the geometric pattern of fractal grids. In addition, an EDQNM simulation with a single excited wavenumber has been produced. This last case should be a closer representation of the production mechanisms of classical grids, where only one mode is excited. For the analysis performed in this section, the value for the parameter α has been fixed to 0.5. The evolution of the spectra E(k, t) and the non-linear energy transfer budget kT(E; k, t) are shown in Figure 11 for the case r fD = 3. The analysis here is similar to the study of the continuous forcing case, and results are reported for Re M /Re L = 35, 45, 60. The maximum energy variation is obtained for the excited mode closer to the inertial range, which arguably sets the stage for the non-classical evolution observed in the physical quantities. However, the qualitative evolution of the energy spectra is very similar to the cases presented in Figure 6, which supports the idea that the continuous forcing is an asymptotic limit for the representation of production mechanisms.

Results for the time evolution of the global dissipation coefficient C ε are shown in Figure 12. Generally speaking, the evolution observed for the continuous forcing and the discrete forcing are similar. In particular, the comparison of the continuous forcing/discrete forcing (r f = 3 and r fD = 3) indicates that the return towards a constant behaviour is slightly slower. This aspect is associated with less intense action of the production mechanisms, which result in a smaller energy dissipation as shown in Figure 13. In addition, it is worth noting that the evolution for r fD = 1 is the most distant from the classical simulation. Assuming that the results from this EDQNM calculation represent regular grid turbulence, this observation suggests that classical simulation of HIT free decay should be not considered as the closest numerical representation of grid turbulence unless a realistic initial energy spectrum conserving the features of the basic production mechanisms is imposed.

The analysis of the initial decay of (K 0 + K) and (ε 0 + ε) confirms initial qualitative guesses. The results, which are shown in Figure 14, indicate that the intensity of the power law exponents decreases when a discrete forcing is considered. This result emerges from the weaker local interaction among very energetic modes associated with E(k, t), which is responsible of a higher turbulent kinetic energy dissipation. At last, the time evolution of the normalised energy flux (k, t)/ε > (k, t) is investigated for the case r fD = 3 as shown in Figure 15. The wavenumbers are chosen following the same strategy previously proposed. The observation of the results shows that similar distribution with respect to the continuous forcing case is obtained for every wavenumber. However, ε > (k, t) is lower for the discrete forcing, so that the resulting energy flux will be less intense as well.

Conclusions

The effects of turbulent kinetic energy production on isotropic turbulence decay have been investigated via EDQNM simulations. Observations reported in wind-tunnel experiments dealing with turbulence decay in the wake of fractal grids and fractal bluff bodies are obtained, including scaling laws for the dissipation parameter C ε . Detailed analysis of kinetic energy spectrum E(k) and energy transfer T(k) shows that anomalous decay regimes are associated with the relaxation of initial energy spectra with a bump at energetic scales toward solutions with a smooth shape at the spectrum peak. The existence of a bump in E(k) induces an increase in the energy cascade rate, in agreement with previously existing theoretical and numerical results dealing with triadic interactions in isotropic turbulence. Fractal/multiscale objects are associated to a injection of turbulent kinetic energy over a spectral band due to the multiscale velocity shear that occurs in their wake, which magnifies the observability of anomalous decay regimes. The relaxation toward classical solutions is due to a double relaxation phenomena, namely (1) the vanishing of the multiscale shear and the multiscale energy production and (2) the smoothing of the spectrum peak due to the increase of the energy cascade associated with the peak itself. The fact that these phenomena, including the correct scaling laws, are captured by the EDQNM model shows that arguably they are not due to breakdown of isotropy or homogeneity, but only to peculiar features of kinetic energy production and induced energy spectrum shapes.

Figure  .Figure  .

  Figure . Perturbed energy spectrum E. The vertical grey dashed line represents the position of the Kolmogorov scale wavenumber k η . In this representation D F =   and Re λ =   .

Figure  .

  Figure . Evolution of the energy spectrum E(k, t) via EDQNM calculations. The case here represented is for α =   , β = . and r f = . Times are normalised over the initial turnover time t  .

Figure  .

  Figure . Evolution of the production (via forcing term F) to dissipation ratio in early stages of decay.

Figure  .

  Figure . Evolution of the coefficient C ε , when a multiscale forcing is applied. EDQNM data are compared with grid turbulence data (bi-planar grid) by Valente and Vassilicos. Data are reported for values of the model parameters of α = ., , . and r f = .

Figure  .

  Figure . (left column) Energy spectra E(k, t) and (right column) energy transfer budget k • T(k, t). The results are shown for (first row) forced simulation for α = ., r f = , (second row) forced simulation for α = ., r f =  and (third row) classical EDQNM solution with F = .

Figure  .

  Figure . Evolution of the power law exponents characterising HIT decay with respect to Re M /Re L . Exponents for (a) the turbulent kinetic energy n K and (b) the energy dissipation rate n ε are shown, respectively.

Figure  .

  Figure . Time evolution of the power law exponents characterising HIT decay. Exponents for (a) the turbulent kinetic energy n K and (b) the energy dissipation rate n ε are shown, respectively.

Figure  .

  Figure . evolution of the coefficient C ε calculated via EDQNM model. The initial transient and (b) long time decay are shown, respectively.

Figure  .

  Figure . Time evolution of the normalised energy flux (k, t)/ε > (k, t) calculated via EDQNM. (a) Classical free decay and forced decay with (b) r f =  and (c) r f =  are shown, respectively.

Figure

  Figure . (a) Energy spectra E(k, t) and (b) non-linear energy transfer budget k • T(k, t) calculated via EDQNM simulations, in the case of discrete forcing.

Figure

  Figure . Early time evolution of the coefficient C ε . Comparison between classical, continuous forcing and discrete forcing EDQNM data is provided.

Figure

  Figure . Early time evolution of the forcing term F. Discrete forcing and continuous forcing are compared for α = ..

Figure  .Figure

  Figure . Time evolution of the power law exponents characterising HIT decay, in the case of discrete forcing. Exponents for (a) the turbulent kinetic energy n K and (b) the energy dissipation rate n ε are shown, respectively.

  last, the energy flux defined as (k, t ) =

			+∞ k	[T (E; p, t ) + T (E, E; p, t )]dp is
	investigated. It can also be split into a base and a disturbance component:
	(k, t ) = (E; k, t ) +	(E, E; k, t )
	+∞		+∞
	=	T (E; p, t )dp +	T (E, E; p, t )dp	(22)
	k		k

Results dealing with the normalised flux across mode k, i.e. (k, t)/ε > (k, t), with ε > (k, t ) = +∞ k 2ν p 2 E(p, t )dp are reported in Figure

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendix. The EDQNM model

In this appendix, elements of the derivation of the eddy-damped quasi-normal Markovian model are provided. Detailed discussion and analysis is reported in the literature [START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF]. This model allows for a very accurate description of the non-linear energy transfers in isotropic turbulence via direct investigation of the statistical evolution of the energy spectrum. The instantaneous velocity field is not simulated, so that the computational resources required are negligible when compared to DNS. In addition, convergence of the statistical moments is granted. These favourable characteristics indicate that the EDQNM model is a very powerful complementary tool to DNS, which allows for the screening of a very large number of configurations and providing reliable information of HIT statistical behaviour. In recent years, this model has been used to investigate HIT evolution [START_REF] Meldi | On non-self-similar regimes in homogeneous isotropic turbulence decay[END_REF][START_REF] Meldi | Further insights into self-similarity and self-preservation in freely decaying isotropic turbulence[END_REF] as well as features of the pressure spectrum [START_REF] Meldi | The signature of initial production mechanisms in isotropic turbulence decay[END_REF][START_REF] Meldi | Pressure statistics in self-similar freely decaying isotropic turbulence[END_REF] and the transport of scalars [START_REF] Briard | Mixed-derivative skewness for high prandtl and reynolds numberd in homogeneous isotropic turbulence[END_REF].

The time evolution of E(k, t) and T(k, t) is investigated solving Equation (1) over a spectral mesh discretisation (wavenumbers k). Because of the nonlinear nature of turbulence, the physical problem represented by Equation ( 1) is not closed. Thus, the model relies on two fundamental approximations for the calculation of the term T(k, t):

(1) The deviation of the fluctuating velocity probability density function (pdf) from a Gaussian distribution, which is usually referred to as fourth-order cumulant, is accounted via the addition of a linear eddy-damping term. (2) The characteristic time associated with the eddy-damping rate is assumed to be negligible with respect to the characteristic evolution time (Markovianisation).

The use of the previous hypotheses provides a dramatic simplification of the expression of the non-linear energy transfer T(k, t), which can be calculated by the algebraic relation:

g is a geometric function and the wavenumbers [p, q, k] describe a triad in the spectral space.

We usually refer to local triadic interactions when the magnitude of p, q and k is similar. Conversely, non-local interactions are observed when the wavenumbers exhibit important differences, up to several orders of magnitude in size. The term kpq is the relaxation time derived via the Markovianisation hypothesis:

where a 1 is a model constant. The geometric coefficient g in Equation (A1) is calculated as

The isotropic forcing term F(k, t) in Equation ( 1) represents the lingering effect of production mechanisms in HIT decay. It is designed to mimic the effect of turbulence production mechanisms in the wake of fractal grids. This term has been modelled for inclusion in the EDQNM calculation starting from the proposal by Mazzi and Vassilicos [START_REF] Mazzi | Fractal-generated turbulence[END_REF], which has been developed for DNS of statistically steady isotropic turbulence. In addition, a time evolution of the model is included via a parameter α to account for the evolution of the wake of the grid used to generate turbulence in wind-tunnel experiments, as proposed by the authors in previous studies [START_REF] Meldi | On the emergence of non-classical decay regimes in multiscale/fractal generated isotropic turbulence[END_REF][START_REF] Meldi | The signature of initial production mechanisms in isotropic turbulence decay[END_REF]. The algebraic expression for the forcing, which is applied to every mesh element in a limited range of the spectral space k ࢠ [k f , k F ] (continuous forcing) or just to selected modes in the same range (discrete forcing), is

where β is a shape parameter driving the distribution of the forcing while α drives the time evolution of the forcing for times larger than a virtual origin τ . In addition, t 0 = K(0)/ε(0) is the initial turnover time and L τ = L γ = k -1 f are the characteristic lengths of the forcing. The forcing is applied in a limited range of the spectral mesh, which is controlled by the coefficient a k . In fact, a k = 1 for k ࢠ [k f , k F ], a k = 0 elsewhere. For the discrete forcing, the coefficient a k = 1 only for four selected modes, which are distributed following a geometric progression in the range of investigation [k f , k F ]. In this work, the forcing has been applied for every time step of the EDQNM calculation, i.e. the forcing is continuous in time. Finally, the constant 2 f γ R allows to set the initial intensity of the forcing. In this case, this parameter is chosen such that

Fdk is 10 times larger than the energy dissipation rate ε(0).