# Equilibrium and non-equilibrium statistics in inhomogeneous and unsteady turbulence

Wouter J T Bos, Ryo Araki

► To cite this version:

Wouter J T Bos, Ryo Araki. Equilibrium and non-equilibrium statistics in inhomogeneous and unsteady turbulence. Physical Review Fluids, 2025, 10 (4), pp.044603. 10.1103/PhysRevFluids.10.044603. hal-04967602

HAL Id: hal-04967602

https://hal.science/hal-04967602v1

Submitted on 26 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Equilibrium and non-equilibrium statistics in inhomogeneous and unsteady
turbulence

Wouter J. T. Bos¹ and Ryo Araki¹,²

1) Univ Lyon, École Centrale de Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, Laboratoire de Mécanique des Fluides et d'Acoustique, UMR 5509, 36 Avenue Guy de Collongue, F-69134, ECULLY, France<sup>a)</sup>

2) Current address: Department of Mechanical and Aerospace Engineering, Faculty of Science and Technology, Tokyo University of Science, Yamazaki 2641, Noda-shi 278-8510, Japan

(Dated: 26 February 2025)

We model the turbulent kinetic energy of an arbitrary turbulent flow as an equi-
librium plus a perturbation. Subsequently, perturbation theory is applied to un-
derstand the nature of both temporal and spatial fluctuations in either unsteady or
inhomogeneous flows. The models are assessed using numerical simulations of the
three-dimensional Kolmogorov flow. Implications for the formulation of turbulence
models are discussed.

<footnote>a) Electronic mail: wouter.bos@ec-lyon.fr</footnote>
# I. INTRODUCTION

The simplest model for the dissipation rate of turbulent kinetic energy $\epsilon$ in a turbulent flow is Taylor's relation$^{1,2}$

$$ \epsilon(\mathbf{x}, t) \sim \frac{k(\mathbf{x}, t)^{3/2}}{L(\mathbf{x}, t)}, \qquad (1) $$

where $k$ is the turbulent kinetic energy and $L$ a characteristic lengthscale set by the turbulence generating mechanism, or the size of the energetic eddies. Expression (1) gives an order of magnitude estimate of the dissipation rate, once $k, L$ are known.

To assess the validity of expression (1), it is helpful to focus on the quantity

$$ C_{\epsilon}(\mathbf{x}, t) \equiv \frac{\epsilon(\mathbf{x}, t)L(\mathbf{x}, t)}{k(\mathbf{x}, t)^{3/2}}, \qquad (2) $$

also called normalized dissipation rate. This quantity might behave differently in distinct flows$^3$. But for a given flow, if suitable averages for $k, \epsilon$ and $L$ are used, its value tends, in general, to a constant$^{4-6}$ for high Reynolds numbers. However, the instantaneous and/or local values of the dissipation rate deviate in general from Taylor's estimate.

Nevertheless, a certain universality is observed for the temporal fluctuations of ensemble averages or volume averages of the temporal variations of the normalized dissipation rate$^{7-9}$ (see the review$^{10}$), and similar universality is observed for time-averages of the normalized dissipation rate in spatially inhomogeneous flows$^{11}$. Understanding this quasi-universal behavior of the quantity $C_{\epsilon}(\mathbf{x}, t)$ and its implications for turbulence modeling are the main subjects of the present investigation.

A related subject is the wavenumber scaling of the kinetic energy spectrum. Indeed, one of the most influential contributions to turbulence theory is Kolmogorov's inertial range concept$^{12,13}$, predicting that the energy spectrum scales proportional to $\kappa^{-5/3}$ in an intermediate range of wavenumbers $\kappa$. Recently, interest has revived in the characterization of perturbations to Kolmogorov's inertial range scaling. Yoshizawa$^{14}$ predicted that unsteady perturbations of the kinetic energy distribution satisfy, in the inertial range, a scaling proportional to $\kappa^{-7/3}$. Theoretical support was presented in$^{15,16}$. The prediction was assessed in direct numerical simulations$^{17-19}$ and closure$^{20}$. The case of spatial inhomogeneity was recently addressed and it was shown that a similar approach leads also to a $\kappa^{-7/3}$ correction for this case$^{21}$. A recent review is dedicated to this spectral perturbation approach$^{22}$ and an extension to surface quasi-geostrophic turbulence is presented in Ref.$^{23}$.
Interestingly, a direct relation was pointed out between the unsteady scaling of the dissipation rate and these perturbations of the energy spectrum$^{24}$. Integrating the inertial range energy spectrum associated to equilibrium, and the one associated to its perturbations, it was shown that the fluctuations of the normalized dissipation rate can be quantitatively explained by spectral perturbation theory. The implications for the dissipation rate and integral lengthscale were also verified in grid-turbulence, wake-flow simulations and atmospheric measurements$^{25-28}$. Although results were supportive, this approach obviously raises the question whether we can effectively extend the inertial range approach to the largest scales of the flow. The present investigation attempts to get around this question by directly applying the perturbation approach to the integral quantities.

In the following section, we introduce the perturbation approach and apply it to the kinetic energy balance. Both temporal and inhomogeneous fluctuations are considered. The predictions are assessed in Sec. III. After this validation we show in Sec. IV how the analytical expressions allow us to understand non-equilibrium fluctuations of the normalized dissipation rate. Sec. V discusses the implications of our results for turbulence modeling before concluding in Sec. VI. In the appendices we discuss the influence of a decomposition into large and small scales as well as the influence of the choice of the integral lengthscale in defining the normalized dissipation rate.

## II. DERIVATION OF THE PREDICTIONS

We focus in the present investigation on the integral values of turbulent statistics. In a flow $u(x, t)$ with an (ensemble) average velocity $U(x, t) = \langle u(x, t) \rangle$ and fluctuations $u' \equiv u - U$, the turbulent kinetic energy is defined

$$k(x, t) = \frac{1}{2} \langle |u'(x, t)|^2 \rangle. \qquad (3)$$

In the numerical experiments in Sec. III we will, instead of ensemble averages, use both temporal and spatial averages. Where needed we specify which average is considered.

The evolution equation of the turbulent kinetic energy is given by

$$\frac{\partial k(x, t)}{\partial t} + U_m(x, t) \frac{\partial k(x, t)}{\partial x_m} = p(x, t) - \epsilon(x, t) - \nabla \cdot \phi(x, t), \qquad (4)$$

where $p, \epsilon, \phi$ are, respectively, the production, dissipation and spatial flux of kinetic energy. When we consider spatial inhomogeneity in the present investigation, we restrict ourselves
<figure><img src="image_2.png" /><figcaption>FIG. 1. Sketch of the decomposition into large and small scales. Represented is the energy spectrum E as a function of the wavenumber κ. The flux Π of kinetic energy transfers the energy between scales associated with wavenumbers smaller and larger than $κ_{\Pi}^{-1}$, respectively.</figcaption></figure>

to the case where statistical inhomogeneity is present in one direction only, perpendicular to the mean-flow direction. Thanks to this simplification, the second term on the left-hand side can be neglected (see for instance the appendix of [29] for definitions of p, ϵ, φ for this specific case). Typical examples are channel flow, planar Couette flow and the case we will consider in Sec. III, Kolmogorov flow. In such flows, where mean-flow and statistics only vary in a single direction, which we call z, the energy balance simplifies to

$$ \frac{\partial k(z, t)}{\partial t} = p(z, t) - \epsilon(z, t) - \frac{\partial}{\partial z} \phi_z(z, t). \qquad (5) $$

In this equation we have

$$ p(z, t) = -\langle u'w' \rangle \frac{\partial U}{\partial z} \qquad (6) $$

$$ \phi_z(z, t) = \langle w'\pi' \rangle + \frac{1}{2} \langle u'_i u'_i w' \rangle - \frac{\nu}{2} \frac{\partial}{\partial z} \langle u'_i u'_i \rangle \qquad (7) $$

$$ \epsilon(z, t) = \nu \langle \frac{\partial u'_i}{\partial x_j} \frac{\partial u'_i}{\partial x_j} \rangle. \qquad (8) $$

where ν is the kinematic viscosity, π' a pressure fluctuation and $u' = (u', v', w')$ denote x, y, and z components of the velocity fluctuations, respectively.

In statistically steady and statistically homogeneous flows the average energy production $\langle p \rangle$ is equal to the average dissipation $\langle \epsilon \rangle$. In general, the production and dissipation act at different lengthscales and positions, and are therefore not equal if we evaluate them locally
and instantaneously. We are interested in the present study exactly in these deviations from the equilibrium $p = \epsilon$.

To take into account that different scales do not react simultaneously, we introduce a two-scale description, $k = k^< + k^>$ using an appropriate filter. When considering volume averaged periodic domains, one can conveniently use a sharp spectral low-pass filter (see Fig. 1). In such case, the filter is determined by a cut-off wavenumber $\kappa_{\Pi}$ and its associated lengthscale is $l \sim \kappa_{\Pi}^{-1}$.

For inhomogeneous flows, the choice of a filter is more delicate, but spectral filters can still be used, either invoking scale separation$^{30,31}$, or by applying the Fourier-description in the homogeneous directions only$^{21}$. Alternatives would consist in using filters more adapted to inhomogeneous flows as discussed in ref.$^{32,33}$. In the present investigation the exact type of filter is not important, as long as we can distinguish small and large scales which exchange energy by means of an interscale flux $\Pi$.

We consider large Reynolds number flow, with scale separation between production and dissipation ranges. Considering that the filter scale $l$ is chosen in between the production and dissipation ranges, allows us to write separate equations for the large- and small-scale energy balance,

$$ \frac{\partial k^<(z, t)}{\partial t} = p(z, t) - \Pi(z, t) - \left( \frac{\partial}{\partial z} \phi_z(z, t) \right)^<, \quad (9) $$

$$ \frac{\partial k^>(z, t)}{\partial t} = \Pi(z, t) - \epsilon(z, t) - \left( \frac{\partial}{\partial z} \phi_z(z, t) \right)^> . \quad (10) $$

The notion of equilibrium in turbulence is on this level represented by the balance between the flux, production and dissipation. Large-scale equilibrium can be written

$$ p(z, t) = \Pi(z, t), \quad (11) $$

and small scale equilibrium

$$ \Pi(z, t) = \epsilon(z, t). \quad (12) $$

It is a matter of taste, or considered phenomenon, which of these two equilibria is chosen as the reference case for a perturbation analysis and we will discuss both. Indeed, not only the large scales, but also the small scales can be out of equilibrium as was illustrated in numerical experiments$^{34}$. We introduce as the notation for equilibrium kinetic energy the subscript zero and for the perturbations the subscript one, $k = k_0 + k_1$. In the present
investigation the first (unsteady) and last (inhomogeneous) terms of Eq. (5) are considered to give rise to perturbations $k_1$ around the equilibrium $k_0$. Note that one should not confuse the large vs. small scale decomposition (using superscripts <, >) and the equilibrium vs. non-equilibrium decomposition (using subscripts 0, 1). Indeed, even without separating the large and the small scales, one can apply a perturbation approach to the system as is shown in the appendix.

To get definite predictions for the behavior of the non-equilibrium behavior of the kinetic energy, we make several assumptions,

* The small scales adapt rapidly to the large-scales so that

$$k^>(z, t) = \gamma k^<(z, t), \qquad (13)$$

where $\gamma$ is a positive constant for a given flow. The precise value of $\gamma$ depends on the type of flow and the ratio of the filter-scale $l$ to the injection scale $L$. Taking a value of $\gamma$ close to zero corresponds to a value $l/L$ tending to zero.

* The inter-scale energy flux is locally (in time and space) determined by the large-scale kinetic energy

$$\Pi(z, t) \sim \frac{(k(z, t))^3}{L} \qquad (14)$$

This expression is the flux-equivalent of Eq. (1). Since the energy-input scale is fixed in the present investigation, we consider that $L$ is constant. Because the flux (for large filter-scales) is mainly determined by the large scales, the validity of expression (14) is in general better satisfied than Eq. (1) (see e.g. [9]).

The introduction of a constant value of $\gamma$ is a common assumption in the turbulence modeling literature to reduce the model complexity. See for instance Ref. [35]. It can and will be dependent on flow-type and filter size. In appendix A we elaborate on the validity of the assumption of constant $\gamma$. It is shown that it is a valid assumption for the case of a fixed lengthscale $L$, which is the case we consider.

A direct consequence of relation (13) is that $k = (1 + \gamma)k^<$, so that we can write the energy balance as

$$\frac{\partial k(z, t)}{\partial t} = (1 + \gamma) \left[ p(z, t) - \Pi(z, t) + \left( \frac{\partial}{\partial z} \phi_z(z, t) \right)^< \right]. \qquad (15)$$
and the flux as

$$ \Pi(z, t) = C_{\Pi} \frac{k(z, t)^{3/2}}{L}. \qquad (16) $$

In this expression $C_{\Pi}$ is a constant, in contrast to the normalized dissipation rate $C_{\epsilon}(z, t) \equiv \epsilon(z, t)L/k(z, t)^{3/2}$, which can depend on space and time.

At equilibrium, where $p = \Pi$, (16) gives,

$$ k_0^{(p)}(z, t) = C_{\Pi}^{-2/3} (Lp(z, t))^{2/3}. \qquad (17) $$

An alternative equilibrium is obtained by starting from $\Pi = \epsilon$,

$$ k_0^{(\epsilon)}(z, t) = C_{\Pi}^{-2/3} (L\epsilon(z, t))^{2/3}. \qquad (18) $$

The superscripts $(p, \epsilon)$ indicate whether we use the production or dissipation of kinetic energy to determine equilibrium. For the sake of clarity, these superscripts are omitted in the derivations but we will keep them in the final expressions.

We focus first on unsteady perturbations to these equilibrium states.

## A. Unsteady perturbations

We consider volume averages so that the space dependence and the inhomogeneous transport $\partial_z \phi_z$ vanish. As stated in the foregoing, we write the kinetic energy as an equilibrium plus a perturbation, $k = k_0 + k_1$,

$$ \frac{d(k_0(t) + k_1(t))}{dt} = (1 + \gamma) \left[ p(t) - C_{\Pi} \frac{(k_0(t) + k_1(t))^{3/2}}{L} \right]. \qquad (19) $$

Retaining only leading order terms on both sides and using that $p - \Pi[k_0] = 0$, we obtain

$$ \frac{dk_0(t)}{dt} = -(1 + \gamma) \frac{C_{\Pi}}{L} \left[ \frac{3}{2} k_0(t)^{3/2} \frac{k_1(t)}{k_0(t)} \right], \qquad (20) $$

where we have linearized the term $(1 + k_1/k_0)^{3/2}$ assuming that $k_1 \ll k_0$. We rewrite this as

$$ k_1(t) = - \frac{2L}{3(1 + \gamma)C_{\Pi}} k_0(t)^{-1/2} \frac{dk_0(t)}{dt}. \qquad (21) $$

Substituting Eq. (17), this yields

$$ k_1^{(p)}(t) = -\frac{4}{9} \frac{1}{(1+\gamma)} \left(\frac{L}{C_{\Pi}}\right)^{4/3} \frac{\dot{p}(t)}{p(t)^{2/3}}, \qquad (22) $$

where $\dot{p} = dp/dt$. This expression assumes, as we mentioned, a fixed lengthscale. In the flow we will investigate, $L$ is the forcing scale of the simulation. The use of a different definition of the lengthscale, will be discussed in the Appendix.
## B. Inhomogeneous perturbations

In determining the inhomogeneous correction there is one supplemental difficulty. Whereas in the unsteady case, the temporal perturbations are determined by the time-derivative $dk/dt$, which is directly expressed as a function of $k(t)$, the inhomogeneous transfer $\partial_z\phi_z$ is not an explicit function of the kinetic energy. We therefore need to model $\phi_z$ as a function of $k$. We take the simplest possible model, where the flux is modeled as a gradient of the considered quantity (see e.g. [36]),

$$ \phi_z(z,t) = -\mu_T(z,t) \frac{\partial k(z,t)}{\partial z}. \qquad (23) $$

We assume that the turbulent diffusion $\mu_T$ is dominantly determined by the large scales³¹ and that its value is thereby insensitive to the filtering operation.

Now that we have a closed expression, the derivation of the inhomogeneous correction is analogous to that of the unsteady correction. We consider the balance (15) in the statistically stationary case

$$ \frac{d}{dz} \left( \mu_T(z) \frac{d(k_0(z) + k_1(z))}{dz} \right) = (1 + \gamma) \left[ p(z) - C_{\Pi} \frac{(k_0(z) + k_1(z))^{3/2}}{L} \right]. \qquad (24) $$

We follow the same procedure as in the previous section. The only difference with the unsteady perturbation is that we have replaced

$$ \frac{dk}{dt} \rightarrow \frac{d}{dz} \left( \mu_T \frac{dk}{dz} \right). \qquad (25) $$

Then, substituting again the equilibrium value, retaining leading order contributions and linearizing, the expression for $k_1(z)$ is

$$ k_1(z) = - \frac{2}{3C_{\Pi}(1+\gamma)} \frac{L}{k_0(z)^{1/2}} \frac{d}{dz} \left( \mu_T(z) \frac{dk_0(z)}{dz} \right). \qquad (26) $$

For the turbulent diffusion we use a classical model $\mu_T \sim \mathcal{UL}$ based on typical velocity $\mathcal{U} \sim k^{1/2}$ and mixing-length $\mathcal{L} \sim L$,

$$ \mu_T(z) = \frac{C_\mu}{C_{\Pi}} k(z)^{1/2} L, \qquad (27) $$

where the constants are chosen such that $C_\mu$ coincides with its definition in the modeling literature³⁶. Substituting again the equilibrium value, retaining leading order contributions and linearizing, the expression for $k_1(z)$ becomes

$$ k_1^{(p)}(z) = -\frac{4 C_\mu}{9(1+\gamma)} \left(\frac{L}{C_{\Pi}}\right)^{8/3} \frac{p_{zz}(z)}{p(z)^{1/3}}, \qquad (28) $$
where we used that

$$ \frac{d}{dz} \left( \mu_T [k_0(z)] \frac{dk_0(z)}{dz} \right) = \frac{2}{3} \left( \frac{C_\mu L}{C_\Pi} \right) k_0(z)^{3/2} \frac{p_{zz}(z)}{p(z)}, \quad (29) $$

and where we introduced the short-hand notation $p_{zz} \equiv d^2p/dz^2$.

### C. Expressions for large-scale and small-scale perturbation theory

To avoid repetitions of the expressions for the non-stationary and inhomogeneous cases we introduce the quantity $\Omega[p(z, t)]$, where

$$ \Omega[p(t)] = \frac{\dot{p}(t)}{p(t)} (1 + \gamma)^{-1} \quad (30) $$

and

$$ \Omega[p(z)] = C_{\mu} \left( \frac{L}{C_{\Pi}} \right)^{4/3} \frac{p_{zz}(z)}{p(z)^{2/3}} (1 + \gamma)^{-1}. \quad (31) $$

We can then write

$$ k^{(p)}(z, t) = C_{\Pi}^{-2/3} p(z, t)^{2/3} L^{2/3} - \frac{4}{9} C_{\Pi}^{-4/3} \Omega[p(z, t)] p(z, t)^{1/3} L^{4/3} \quad (32) $$

Up to here, the kinetic energy was expressed as a function of the energy input $p(z, t)$ and its derivatives. This corresponds, as we mentioned, to perturbations to large scale equilibrium. Exactly the same derivations can be carried out by perturbing about small scale equilibrium $\Pi(z, t) = \epsilon(z, t)$. Introducing $\epsilon_{zz} = d^2\epsilon/dz^2$ and

$$ \Omega[\epsilon(t)] = \frac{\dot{\epsilon}(t)}{\epsilon(t)} \frac{\gamma}{(1 + \gamma)} \quad (33) $$

and

$$ \Omega[\epsilon(z)] = C_{\mu} \left( \frac{L}{C_{\Pi}} \right)^{4/3} \frac{\epsilon_{zz}(z)}{\epsilon(z)^{2/3}} \frac{\gamma}{(1 + \gamma)}, \quad (34) $$

the expression for the kinetic energy is

$$ k^{(\epsilon)}(z, t) = C_{\Pi}^{-2/3} (\epsilon(z, t)L)^{2/3} + \frac{4}{9} C_{\Pi}^{-4/3} \Omega[\epsilon(z, t)] \epsilon(z, t)^{1/3} L^{4/3}. \quad (35) $$

This is the integral version of Yoshizawa's unsteady correction to inertial range scaling¹⁴ and the equivalent correction in inhomogeneous turbulence²¹. The non-equilibrium corrections obtained starting from $\Pi = \epsilon$ (Eq. (35)) are similar in shape to those obtained by perturbing about a $p = \Pi$ equilibrium (Eq. (32)). Note however that the sign of the correction is
<table><thead><tr><th>N</th><th>ν</th><th>u'</th><th>λ</th><th>R<sub>λ</sub></th><th>T<sub>total</sub>/T</th></tr></thead><tbody><tr><td>512</td><td>0.01</td><td>1.33</td><td>0.138</td><td>102</td><td>170</td></tr></tbody></table>

TABLE I. Control Parameters of the direct numerical simulations of three dimensional Kolmogorov flow: resolution N, kinematic viscosity ν. Key statistical quantities: fluctuating isotropic RMS velocity $u' \equiv \sqrt{2k/3}$ where energy of the temporal fluctuating velocity $k \equiv \langle u'_i u'_i \rangle_{x,t}/2$ and $u'_i(x,t) \equiv u_i(x,t) - \langle u_i \rangle_t(x)$; the Taylor microscale $\lambda \equiv u'\sqrt{15\nu/\epsilon}$ where the energy dissipation rate is evaluated by $\epsilon = \nu\langle\omega_i\omega_i\rangle_{x,t}$; the Taylor-length Reynolds number $R_\lambda \equiv u'\lambda/\nu$; the integral time scale $T \equiv L/u'$ with $L = k_f^{-1} = 1$, with $k_f$ the forcing wavenumber; the simulation time in the statistically steady state $T_{\text{total}}$ as a function of $T$.

opposite, reflecting the fact that the dissipation lags the production of kinetic energy. Note further that the dependence on $\gamma$ is slightly different, since

$$k = (1 + \gamma)k^{<} = \frac{1 + \gamma}{\gamma}k^{>}. \qquad (36)$$

For convenience, the different expressions are summarized in table II in the conclusions of this article. We conclude this section, by stressing that the introduction of the large-small description is not essential. Even in the limit $\gamma$ tending to zero, the perturbation analysis can be carried out. This is explained in the appendix.

### III. NUMERICAL VERIFICATION OF THE NON-EQUILIBRIUM CORRECTIONS

We use a data-base of results of numerical simulations of three-dimensional Kolmogorov flow. The numerical code is pseudo-spectral$^{37}$ and we consider periodic boundary conditions. The resolution is 512$^3$ grid-points, viscosity is $\nu = 0.01$, corresponding to an average Reynolds number $R_\lambda = \sqrt{20k^2/(3\epsilon\nu)} = 102$. Note that if we add the contribution of the mean-flow $U$ to the kinetic energy, the value of the Reynolds number increases to $R_\lambda = 184$. The parameters of the flow are summarized in table 1. Further details of the flow, such as energy spectra and visualizations can be found in [21]. All statistics are reported during a statistically steady state.

To determine $C_\Pi$, we proceed as follows: since in a steady state the time- and space-averaged flux should equal the dissipation $\langle p \rangle_{t,z} = \langle \Pi \rangle_{t,z} = \langle \epsilon \rangle_{t,z}$, the average dissipation
<figure><img src="image_3.png" /><figcaption>FIG. 2. Comparison of the DNS results for the volume averaged kinetic energy in Kolmogorov flow to equilibrium and non-equilibrium estimates, using the equilibrium $p(t) = \Pi(t)$. (a) kinetic energy compared to the equilibrium estimate $k_0^{(p)}(t) = (p(t)L/C_\Pi)^{2/3}$. (b) deviation from equilibrium $k(t) - k_0^{(p)}(t)$, compared to $k_1^{(p)}(t) = -(4/9)(L/C_\Pi)^{4/3}(1+\gamma)^{-1}p^{-2/3}(t)\dot{p}(t)$, the non-equilibrium estimate (last term of Eq. (32)).</figcaption></figure>

rate relates to the flux constant by

$$C_{\Pi} = \frac{\langle \epsilon \rangle_{t,z} L}{\langle k^{3/2} \rangle_{t,z}} = \frac{\langle p \rangle_{t,z} L}{\langle k^{3/2} \rangle_{t,z}}. \qquad (37)$$

where the brackets denote in the present case an average over both space and time.

For the lengthscale, we use $L = 1/\kappa_f = 1$, corresponding to the forcing scale. Clearly, this choice is not unique. In Taylor's original work the lengthscale $L$ was chosen to represent
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. As in Figure 2, temporal signals of the kinetic energy and estimates $k_0^{(\epsilon)}(t)$, $k_1^{(\epsilon)}(t)$ obtained from perturbation theory using the equilibrium $\epsilon(t) = \Pi(t)$.</figcaption></figure>

the size of the system, which is consistent with the present choice. In literature, the choice of the lengthscale has been debated³⁸⁻⁴⁰ and there is no definite answer on a correct choice of $L$. In any case, in our description of turbulence, the lengthscale $L$ appears only in linear combination with the constant $C_{\Pi}$, so that any constant factor in the definition of $L$ can be absorbed in the constant $C_{\Pi}$. Using the current definition, we find $C_{\Pi} = 0.14$. In order to compare this with reference values from literature, we note that the average value of the longitudinal correlation lengthscale is $\langle\Lambda(t)\rangle_t = 1.68\kappa_f^{-1}$ and the velocity fluctuations can be related to the kinetic energy by $u'^2 = \frac{2}{3}k$ so that $\langle\epsilon\rangle\langle\Lambda\rangle/\langle u'\rangle^3 = 0.43$, in agreement with
values from literature⁵.

A. Results for temporal fluctuations

In Fig. 2(a) we show $k(t)$ compared to $k_0^{(p)}(t)$, (expression (17)) as a function of time. Clearly, the order of magnitude of $k_0$ is the same as that of $k$, but the instantaneous fluctuations are not perfectly described. Indeed, there is a time-lag between the two signals. Note however that this estimate, which is simply based on Taylor's relation (1) and the assumption $p = \epsilon$ does reproduce the data already considerably well. This is an indication that the turbulence is not far from equilibrium.

The difference between the two signals is shown in Fig. 2(b). The fluctuations around the equilibrium estimate are compared to the non-equilibrium estimate $k_1^{(p)}(t)$, (expression (22)). It is observed that this quantity describes the data fairly well for a value $\gamma = 0.38$. It is this value which gives the best agreement with the data minimizing the mean-square error between data and estimate. Once $\gamma$ and $C_{\Pi}$ are known, we can thus predict for a given $L$, not only the average, or the equilibrium kinetic energy, but also the instantaneous fluctuations around equilibrium from the time-series of the energy input $p$ into the system.

We can also try to obtain the non-equilibrium estimates from the temporal signal of the dissipation rate $\epsilon(t)$, and perturbing about the small scale equilibrium $\Pi = \epsilon$, i.e., using expression (18). It is observed in Fig. 3(a) that, similarly, the equilibrium estimate follows the data, but the time-lag is now of opposite sign, compared to the result in Fig. 2(a). In Fig. 3(b) we compare the difference between the kinetic energy and $k_0^{(\epsilon)}(t)$ to the non-equilibrium correction $k_1^{(\epsilon)}(t)$, (expression (35)). It is observed that the agreement is less satisfactory than for $k_1^{(p)}(t)$. The reason for this is that the temporal variation of the dissipation rate is less smooth than that of the energy injection. Indeed, considering the ratio of the equilibrium and non-equilibrium estimates, the perturbation approach is valid for $\epsilon L^{2/3} \epsilon^{-4/3} \ll 1$ so that for rapid fluctuations the validity of the approach breaks down. For this reason, using the time-signal of the dissipation rate to estimate the non-equilibrium contribution to the kinetic energy yields less satisfactory agreement with the data than the same procedure using $p$.

If one wants to improve the agreement with data using the dissipation rate signal, one could consider, instead of the instantaneous values, phase-averaged statistics of the dissipation rate, for a well defined forcing protocol$^{20,41}$. Indeed, using a controlled periodic forcing
allows one to study the response of the turbulence for a given large scale frequency. In such an approach the phase-averaged dissipation rate behaves in a smooth fashion$^{42,43}$. We will use this approach, using harmonic forcing, in Section IV to obtain insights on fluctuations of the normalized dissipation rate, but we will now first focus on inhomogeneous non-equilibrium.

## B. Results for inhomogeneous perturbations

The Kolmogorov flow is inhomogeneous in the z-direction. The mean flow and visualizations are given in [21]. Whereas the meanflow is sine-shaped with a period $2\pi/L$ ($k_z = 1$), where $L = 2\pi$ is the domain-size (and the forcing scale), the time-averaged kinetic energy, dissipation and production can all be described, to a good approximation, by harmonic functions with $k_z = 2$ (see also$^{44-46}$ for a discussion of the flow properties). These time-averaged profiles of turbulent kinetic energy, its production and dissipation are shown in Fig. 4(a).

In principle both small scale and large scale equilibria can be used in our perturbation approach, but there is an issue in the present case with the large scale equilibrium: the time-averaged production of kinetic energy is never negative in the present flow, but it can be zero. This is the case in the Kolmogorov flow, where the mean velocity gradient vanishes on two planes perpendicular to the z-axis (see Fig. 4(a) at $z = 0, \pi$). Around these planes, the equilibrium estimate $k_0^{(p)}(z)$ is therefore close, or equal, to zero. A zero equilibrium kinetic energy is not convenient for a perturbation approach since perturbations should be small compared to the equilibrium. Since the perturbation approach breaks down for $p$, we will use the small-scale equilibrium (involving $\epsilon$) to perturb about. Indeed, the dissipation profile is closely described by a smoothly varying function which is not close to zero.

In Fig. 4(a) we show the fit of a cosine-function to the dissipation profile

$$ \bar{\epsilon}(z) = 1.40(1 + 0.125 \cos(2z)). \qquad (38) $$

Since the fit is very close to the actual profile, we will use this analytical profile as an estimate for $\epsilon$, which simplifies the determination of the second spatial derivative $\epsilon_{zz}$ needed to estimate the non-equilibrium part of the kinetic energy.

In Fig. 4(b) we show the equilibrium estimate $k_0^{(\epsilon)}(z)$ (Eq. (18)) compared to the actual kinetic energy profile. The z-averaged value $\langle k \rangle_z$ is correctly reproduced, but the spatial
<figure><img src="image_6.png" /><figcaption>FIG. 4. (a) Averaged profiles of the production $p(z)$, kinetic energy $k(z)$ and dissipation $\epsilon(z)$ in 3D Kolmogorov flow. Standard deviations are indicated by black lines. A dashed line indicates the fit (38) to the dissipation profile. (b) Spatial equilibrium $k_0^{(\epsilon)}(z)$ and non-equilibrium estimate $k_0^{(\epsilon)}(z) + k_1^{(\epsilon)}(z)$ of the kinetic energy, compared to its measured value.</figcaption></figure>

variation of the profile around this value is completely reversed. The reason for this is that $k_0^{(\epsilon)}(z) \sim \epsilon(z)^{2/3}$ and the fluctuation contributions to $k_0^{(\epsilon)}(z)$ and $\epsilon(z)$ have opposite phases.

Adding the non-equilibrium estimate, with $C_\mu = 0.12$ gives excellent agreement between the quantity $k_0(z) + k_1(z)$ and the time-averaged profile $k(z)$ as illustrated in Fig. 4(b). It seems that the combination of our perturbation approach plus the turbulent viscosity
assumption for the inhomogeneous flux allows us to closely reproduce the observations.

## IV. ON THE NORMALIZED DISSIPATION RATE

In the foregoing sections we expressed, quite successfully, the kinetic energy fluctuations as a function of the energy injection, or as a function of the dissipation rate. In this section we will use these relations to compute the normalized dissipation rate and the Taylor-scale Reynolds number for a given temporal or spatial dissipation profile.

To make the theoretical expressions as clear as possible, we will prescribe the temporal or spatial dependence of the dissipation or production. This line of research, also called periodically forced turbulence, was initiated in Ref.⁴⁷ and further explored theoretically, experimentally and numerically in ⁴¹⁻⁴³,⁴⁸⁻⁵⁰, recently an investigation has used periodically forced turbulence to study temporal intermittency of turbulence⁵¹. We will show how this approach allows us to recover analytical expressions for the parametric plots of the normalized dissipation rate.

The main result that we obtain in this section is that we can write the parametric expression of $C_\epsilon(z, t)$ as a function of $R_\lambda(z, t)$ in the form,

$$ \frac{C_{\epsilon}(t)}{\langle C_{\epsilon} \rangle} = 1 - (3a/2) \cos(\omega t), \qquad (39) $$

$$ \frac{R_{\lambda}(t)}{\langle R_{\lambda} \rangle} = 1 + a \cos(\omega t) + b \sin(\omega t), \qquad (40) $$

and

$$ \frac{C_{\epsilon}(z)}{\langle C_{\epsilon} \rangle} = 1 + (3a'/2) \sin(qz), \qquad (41) $$

$$ \frac{R_{\lambda}(z)}{\langle R_{\lambda} \rangle} = 1 - (a' - b') \sin(qz), \qquad (42) $$

where in this section the averages $\langle C_\epsilon \rangle$, $\langle R_\lambda \rangle$ correspond to averages over both space and time.

One can observe that the inhomogeneous and unsteady perturbations do differently alter the parametric relation between $C_\epsilon(z, t)$ and $R_\lambda(z, t)$. Indeed, the unsteady correction, being proportional to the first time-derivative introduces a $\cos(\omega t)$, whereas the second space derivative of the inhomogeneous correction introduces a $-\sin(qz)$ dependence. This difference explains why the inhomogeneous perturbations are approximately captured by a
functional relation $C_\epsilon(z) = F[(R_\lambda(z))]$ (see¹¹ and the present results), whereas the unsteady perturbations describe curved-elliptic orbits in the $C_\epsilon(t) - R_\lambda(t)$ plane²⁰.

Using expression (35), the spatio-temporal fluctuations of the normalized dissipation rate are given by

$$C_{\epsilon}(z, t) = \frac{\epsilon(z, t)L}{(k_{0}^{(\epsilon)}(z, t) + k_{1}^{(\epsilon)}(z, t))^{3/2}} \quad (43)$$

$$= C_{\Pi} \left( 1 - \frac{4 L^{2/3}}{9 C_{\Pi}^{2/3}} \Omega[\epsilon(z, t)] \epsilon(z, t)^{-1/3} \right)^{-3/2} \quad (44)$$

and for the Reynolds number

$$R_{\lambda}(z, t) = \sqrt{\frac{20}{3} \frac{k_{0}^{(\epsilon)}(z, t) + k_{1}^{(\epsilon)}(z, t)}{(\epsilon(z, t)\nu)^{1/2}}} \quad (45)$$

$$= \sqrt{\frac{20}{3} \frac{\epsilon(z, t)^{1/6} L^{2/3}}{\nu^{1/2} C_{\Pi}^{2/3}} \left( 1 - \frac{4 L^{2/3}}{9 C_{\Pi}^{2/3}} \Omega[\epsilon(z, t)] \epsilon(z, t)^{-1/3} \right)}. \quad (46)$$

Note that in the expression for $R_\lambda$ a factor $\epsilon(t, z)^{1/6}$ appears outside the brackets. Without this term, the parametric dependence would be a straight-forward relation $C_\epsilon(z, t) \sim (R_\lambda(z, t))^{-3/2}$. This is thus not the case and the parametric relation is more involved and depends (weakly) on the dissipation profile.

To get a better understanding, we use the dissipation profiles

$$\epsilon(z) = \langle \epsilon \rangle (1 + \alpha_{\epsilon} \sin(qz)); \quad \epsilon(t) = \langle \epsilon \rangle (1 + \alpha_{\epsilon} \sin(\omega t)). \quad (47)$$

We substitute these relations into the expressions for $C_\epsilon(z, t)$, $R_\lambda(z, t)$ yielding

$$\frac{C_{\epsilon}(t)}{C_{\Pi}} = \left( 1 + \frac{4}{9} \frac{\gamma}{(1+\gamma)} \left( \frac{L}{C_{\Pi}} \right)^{2/3} \langle \epsilon \rangle^{-1/3} \alpha_{\epsilon} \omega \cos(\omega t) \right)^{-3/2}, \quad (48)$$

$$\frac{R_{\lambda}(t)}{\langle R_{\lambda} \rangle} = (1 + \alpha_{\epsilon} \sin(\omega t))^{1/6} \left( 1 + \frac{4}{9} \frac{\gamma}{(1+\gamma)} \left( \frac{L}{C_{\Pi}} \right)^{2/3} \langle \epsilon \rangle^{-1/3} \alpha_{\epsilon} \omega \cos(\omega t) \right), \quad (49)$$

and

$$\frac{C_{\epsilon}(z)}{C_{\Pi}} = \left( 1 - \frac{4}{9} C_{\mu} \frac{\gamma}{(1+\gamma)} \left( \frac{qL}{C_{\Pi}} \right)^{2} \alpha_{\epsilon} \sin(qz) \right)^{-3/2}, \quad (50)$$

$$\frac{R_{\lambda}(z)}{\langle R_{\lambda} \rangle} = (1 + \alpha_{\epsilon} \sin(qz))^{1/6} \left( 1 - \frac{4}{9} C_{\mu} \frac{\gamma}{(1+\gamma)} \left( \frac{qL}{C_{\Pi}} \right)^{2} \alpha_{\epsilon} \sin(qz) \right). \quad (51)$$

with

$$\langle R_{\lambda} \rangle = \sqrt{\frac{20}{3} \frac{1}{C_{\Pi}^{2/3}} \frac{\langle \epsilon \rangle^{1/6} L^{2/3}}{\nu^{1/2}}}. \quad (52)$$
<figure><img src="image_7.png" /><figcaption>FIG. 5. Parametric representation of the equilibrium coefficients $C_{\epsilon} = \epsilon L/k^{3/2}$ as a function of the Taylor-Reynolds number $R_{\lambda}$. Comparison of analytical results and numerical observations. (a) Results for temporal fluctuations. (b) Results for spatial variations.</figcaption></figure>

Linearizing the expressions, we find expressions (39), (41) with

$$a = \frac{4}{9} \frac{\gamma}{(1 + \gamma)} \left( \frac{L}{C_{\Pi}} \right)^{2/3} \langle \epsilon \rangle^{-1/3} \alpha_{\epsilon} \omega \quad (53)$$

$$b = b' = \alpha_{\epsilon}/6 \quad (54)$$

$$a' = \frac{4}{9} C_{\mu} \frac{\gamma}{(1 + \gamma)} \left( \frac{qL}{C_{\Pi}} \right)^2 \alpha_{\epsilon}. \quad (55)$$

We plot in Fig. 5 expressions (48)-(51), together with the numerical data. The parameters
used in these plots are $C_Π = 0.14$, $γ = 0.38$, $C_μ = 0.12$, $q = ω = 2$ and $α_ε = 0.35, 0.125$ for the unsteady and inhomogeneous case, respectively.

In Fig. 5(a) we see that the analytical expressions describe the temporal data rather well, in particular given the fact that the sine-function is not the actual functional shape of the data. Fig. 5(b) illustrates that the agreement in the inhomogeneous case is even more successful.

## V. TURBULENCE MODELING

In general in studies of the normalized dissipation rate it is stressed that the precise understanding of $C_ε$ is of major importance in turbulence modeling. For that reason we will here highlight the implications of the present study in the context of classical turbulence models. What we want to stress below is that $C_ε(z, t)$ is an output of a well-designed model, where the model constants can be chosen fixed.

We focus first on statistically homogeneous, but unsteady flow. The investigation of the accuracy of two-equation models for this test-case is also reported in$^{52–54}$. Let us assess the most well known model, the $k-ε$ model$^{55}$ for our case,

$$ 
\begin{aligned}
\frac{dk(t)}{dt} &= p(t) - \epsilon(t), \\
\frac{d\epsilon(t)}{dt} &= \frac{\epsilon(t)}{k(t)}(C_{\epsilon1}p(t) - C_{\epsilon2}\epsilon(t)).
\end{aligned}
\qquad (56) $$

In a homogeneous and stationary flow the constants $C_{ε1}, C_{ε2}$ need to be equal to allow for a non-trivial steady-state $p = ε$. We consider the case where $L$ is constant and assume $p(t)$ given.

The question is now what the role is of Taylor’s relation in this type of model. Assuming directly Taylor’s relation $\epsilon(t) = Ck(t)^{3/2}/L$, we can simplify the system to a one-equation model (since $L$ is fixed),

$$ \frac{dk(t)}{dt} = p(t) - C \frac{k(t)^{3/2}}{L}. \qquad (57) $$

We note that this is equivalent to setting $C_{ε1} = C_{ε2} = 3/2$, as can be seen by elimination of $ε$ from (56). We have thus imposed Taylor’s relation, or, in other words, a constant value of $C_ε = \epsilon(t)L/k^{3/2}(t)$.
If we want to allow for variations of $C_\epsilon$, to more closely reproduce observations, we use relations (15) and (16) from Section II. We have then,

$$ \frac{dk(t)}{dt} = (1 + \gamma) \left( p(t) - C_{\Pi} \frac{k(t)^{3/2}}{L} \right). \qquad (58) $$

The difference with the standard $k - \epsilon$ model (57) for this very simple case is the presence of the factor $(1 + \gamma)$ in the expression. Since the balance $dk/dt = p - \epsilon$ is exact, we see that,

$$ \epsilon(t) = (1 + \gamma)\Pi(t) - \gamma p(t). \qquad (59) $$

Alternatively, we see that this model predicts for the normalized dissipation rate,

$$ C_{\epsilon}(t) = (1 + \gamma)C_{\Pi} - \gamma \frac{p(t)L}{k(t)^{3/2}}. \qquad (60) $$

This shows that a simple one-equation model can qualitatively represent temporal non-equilibrium, as reflected by the temporal dependence of $C_\epsilon(t)$ allowed by expression (60).

Let us conclude this discussion by considering spatial inhomogeneity, considering steady flow while inhomogeneous in the $z$-direction. From our results we can conclude that the flux-gradient assumption was sophisticated enough to reproduce the spatial profile of the kinetic energy. Our flow, with harmonic space dependence of the statistics might not be the most challenging test case and further tests using different profiles, as in$^{46}$ could be considered in future investigations. But let us consider this simple case.

We have the energy balance,

$$ 0 = p(z) - \epsilon(z) + \frac{d}{dz} \mu_T(z) \frac{dk(z)}{dz}. \qquad (61) $$

Filtering this and retaining the small scales, we obtain

$$ 0 = \frac{(1 + \gamma)}{\gamma} (\Pi - \epsilon(z)) + \left( \frac{d}{dz} \left( \mu_T(z) \frac{dk(z)}{dz} \right) \right)^{>} . \qquad (62) $$

Exactly the same reasoning, leads then to an expression similar to (60) for the spatial dependence of the normalized dissipation rate in inhomogeneous flows,

$$ C_{\epsilon}(z) = (1 + \gamma)C_{\Pi} - \gamma \frac{p(z)L}{k(z)^{3/2}}. \qquad (63) $$

Again, we observe that if $L$ is fixed, a one-equation turbulence model allows one to represent non-equilibrium. If the injection lengthscale is not fixed, a two-equation model is needed.
<table><caption>TABLE II. Table of main theoretical results. In our simulations, the lengthscale is fixed at $L = 1/k_f = 1$ and the model constants are $C_π = 0.14$; $\gamma = 0.38$; $C_μ = 0.12$.</caption><thead><tr><th>Equilibrium</th><th>k<sub>0</sub></th><th>k<sub>1</sub></th></tr></thead><tbody><tr><td>p(t) = Π(t)</td><td>C<sub>Π</sub><sup>-2/3</sup>(Lp(t))<sup>2/3</sup></td><td>-<sup>4</sup>&frasl;<sub>9</sub> (L/C<sub>Π</sub>)<sup>4/3</sup> <sup>1</sup>&frasl;<sub>(1+γ)</sub> (p&#775;(t)/p(t))p(t)<sup>1/3</sup></td></tr><tr><td>p(z) = Π(z)</td><td>C<sub>Π</sub><sup>-2/3</sup>(Lp(z))<sup>2/3</sup></td><td>-<sup>4</sup>&frasl;<sub>9</sub> (L/C<sub>Π</sub>)<sup>8/3</sup> C<sub>μ</sub> <sup>1</sup>&frasl;<sub>(1+γ)</sub> (p<sub>zz</sub>(z)/p(z)<sup>2/3</sup>)p(z)<sup>1/3</sup></td></tr><tr><td>ε(t) = Π(t)</td><td>C<sub>Π</sub><sup>-2/3</sup>(Lε(t))<sup>2/3</sup></td><td>+<sup>4</sup>&frasl;<sub>9</sub> (L/C<sub>Π</sub>)<sup>4/3</sup> <sup>γ</sup>&frasl;<sub>(1+γ)</sub> (ε&#775;(t)/ε(t))ε(t)<sup>1/3</sup>.</td></tr><tr><td>ε(z) = Π(z)</td><td>C<sub>Π</sub><sup>-2/3</sup>(Lε(z))<sup>2/3</sup></td><td>+<sup>4</sup>&frasl;<sub>9</sub> (L/C<sub>Π</sub>)<sup>8/3</sup> C<sub>μ</sub> <sup>γ</sup>&frasl;<sub>(1+γ)</sub> (ε<sub>zz</sub>(z)/ε(z)<sup>2/3</sup>)ε(z)<sup>1/3</sup>.</td></tr></tbody></table>

## VI. CONCLUSION

The statement that Taylor's relation is the cornerstone of turbulence modeling is obviously correct. While it is true that blindly using it to model the dissipation rate will miss certain non-equilibrium features, since it imposes a fixed value of the normalized dissipation rate $C_ε(x, t)$, it does not yield catastrophic results, as is illustrated in Figs. 2(a) and 4(a). However, these equilibrium estimates will show positive or negative time or space-lags when compared to the actual data, associated with the spatio-temporal dynamics of the energy flux.

First-order perturbation estimates allow us to improve the agreement. One should be careful, however, that the equilibrium does not approach zero, or that the spatio-temporal variations are not too extreme. This is obviously a limitation of any perturbative approach. For smooth large-scale statistics, the perturbation approach appears to be quite successful to improve the agreement between model and observations. In table II we summarize the analytical expressions derived using our perturbation approach. In future studies it would be particularly valuable to assess the inhomogeneous corrections using different forcing profiles in physical space. A relevant test-case could be a Gaussian mean velocity profile as in Ref. 46. Also, the present results are only assessed for a single Reynolds number. At low Reynolds numbers viscous effects could alter the predictions.

The implication for turbulence models might not be enormous. Indeed, Taylor's relation is quite robust if it is formulated for the energy flux, instead of for the dissipation rate. Using the flux-relation in a turbulence model does change the formulation, but not radically.

For instance, the turbulence production term in Eq. (58) contains a supplemental prefactor $(1+\gamma)$. However, the production term itself is in general modeled using an eddy-viscosity
closure, so that this factor could be absorbed into the constant which defines the value of the eddy-viscosity. We could think that the inhomogeneous corrections are important, considering that the results in Fig. 4(b) show that the equilibrium estimate is out of phase with the total kinetic energy. However, this effect is correctly modeled by an elementary gradient-diffusion model based on a mixing length. The essential insight is here that $C_\epsilon(x, t)$ is not the constant which is used in the model, but an output of the model with a simple inhomogeneous transport term and fixed model constants $C_\Pi$ and $C_\mu$: one should not mix up these model constants with the normalized dissipation rate.

The performance of a model might be less sensitive to the physical improvements than the influence of walls, anisotropy and the optimization of fitting constants. Nevertheless, it might be worth to assess whether the non-equilibrium corrections allow one to change the behavior in canonical test-cases such as homogeneous shear and channel flow.

We note that the current investigation does not explain the origin of the fluctuations, neither the amplitude or period of them. The fluctuations are either determined by an externally defined forcing or are inherent to the dynamics of the Navier-Stokes equations. In the present study we have not considered any detailed physics of the system. To advance on the dynamical origin of turbulent fluctuations, one needs an understanding of the large-scale instability mechanisms of the flow and a promising model to advance in this understanding was recently formulated$^{56}$ for flows in the presence of steady forcing. The questions: what is the origin of the fluctuations and how to characterize their statistics seem to be, at the current understanding of turbulent dynamics, disconnected issues.

As a perspective we want to mention that it can be of major interest to investigate the implication of stationarity on the statistics of the variance of the fluctuations. Indeed, recently$^{57}$ a fundamental relation between correlation times and fluctuations$^{58}$ was used to shed some light on intermittent scaling of turbulent fluctuations. Even more recently it was found that this same intermittent scaling is in inhomogeneous flows correlated to some extent to the scaling of the dissipation rate constant $C_\epsilon^{59}$. It would be interesting to test whether these observations can also be explained using perturbation analysis. However, it is not clear how such an approach can be extended to higher-order velocity moments.
# ACKNOWLEDGMENTS

Personal acknowledgment of WB: The present work builds on the insights I gained during a collaboration with Robert Rubinstein. The influence Bob has had on my current understanding of modeling turbulent flows cannot be overstated. thank him as both a friend and a mentor and wish him a long, happy, and musical retirement.

We further thank Andrea Maffioli, Faouzi Laadhari, Felipe Alves Portela and Susumu Goto for insightful comments and stimulating discussions.

All DNS simulations were carried out using the facilities of the PMCS2I (École Centrale de Lyon). R.A. is supported by the Takenaka Scholarship Foundation.

For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

# APPENDIX

## A. On the validity of the large-small decomposition

To get an intuition on the meaning of $\gamma$, consider a spectral distribution,

$$E(\kappa, t) = C\epsilon(t)^{2/3}\kappa^{-5/3}, \qquad (64)$$

ranging from $1/L$ to $1/\eta$, with $\eta \ll L$. Using a sharp cut-off, the energy contained for scales smaller than $l$, compared to those larger, for $L^{-1} \ll \kappa \ll \eta^{-1}$, is then

$$\gamma = \frac{\int_{l^{-1}}^{\infty} E(\kappa, t) d\kappa}{\int_{l^{-1}}^{l^{-1}} E(\kappa, t) d\kappa} = \frac{l^{2/3} - \eta^{2/3}}{L^{-2/3} - l^{2/3}} \approx \left(\frac{l}{L}\right)^{2/3}. \qquad (65)$$

Adding an infra-red range for $\kappa < 1/L$, so that $E(\kappa, t) = C\epsilon^{2/3}L^{5/3}(\kappa L)^s$, this changes the value of $\gamma$ to

$$\gamma \approx \left(1 + \frac{2}{3(s+1)}\right)^{-1} \left(\frac{l}{L}\right)^{2/3}. \qquad (66)$$

which shows that the scaling remains unchanged. The prefactor for $s = 2$ changes from unity to 0.75, which shows that also quantitatively the value of $\gamma$ is quite insensitive to changes in shape of the large scales.
An important question is how $\gamma$ varies in time if the injection scale $L$ is not fixed in time. Let us use expression (66) for the value of $\gamma$, where we vary $L(t)$ for a fixed filter-scale $l$,

$$ \frac{\dot{\gamma}}{\gamma} = -\frac{2}{3} \frac{\dot{L}}{L}. \qquad (67) $$

This shows that the variations of $\gamma$ are weak for the case of a constant forcing-scale.

We can also directly evaluate the error induced in equation (15) by the assumption of a constant value of $\gamma$,

$$ \frac{dk}{dt} = \frac{d}{dt}[(1 + \gamma)k^{<}] = \underbrace{(1 + \gamma)\frac{dk^{<}}{dt}}_{\text{Eq. (15)}} + k^{<}\frac{d\gamma}{dt}. \qquad (68) $$

Evaluating the ratio of the last two terms we find,

$$ \left| \frac{k^{<} \dot{\gamma}}{(1 + \gamma) \dot{k}^{<}} \right| = \frac{\gamma}{1 + \gamma} \left( \frac{1}{1 + \frac{\dot{\epsilon}}{\epsilon} \frac{L}{L}} \right) \qquad (69) $$

This shows that the approximation of constant $\gamma$ is valid if the fluctuations of the production lengthscale $L$ are slow compared to the fluctuations of the energy production/dissipation rate. In our set-up the injection lengthscale is fixed in time, validating the assumption of constant $\gamma$.

## B. Modeling the system without a large-small decomposition

To illustrate that we can use the perturbation approach directly on the full system, we consider the statistically homogeneous case (or a volume-averaged inhomogeneous case),

$$ \frac{dk(t)}{dt} = p(t) - \epsilon(t). \qquad (70) $$

Using Taylor's estimate directly to model the dissipation rate, we write

$$ \frac{dk(t)}{dt} = p(t) - C \frac{k(t)^{3/2}}{L}, \qquad (71) $$

so that the equilibrium estimate (following from $p = \epsilon$) is

$$ k_0(t) = \left( \frac{p(t)L}{C} \right)^{2/3}. \qquad (72) $$

Now substituting $k = k_0 + k_1$ in (71), we find to leading order,

$$ \frac{dk_0(t)}{dt} = -\frac{3}{2}C \frac{k_0^{3/2}(t) k_1(t)}{L} k_0(t), \qquad (73) $$
where we used $k_1 \ll k_0$. Using expression (72), this yields for the perturbation $k_1$,

$$k_1(t) = -\frac{4}{9} \left(\frac{L}{C}\right)^{4/3} \frac{dp(t)}{dt} p(t)^{-2/3}. \qquad (74)$$

This expression corresponds to the case $\gamma = 0$ in the main text.

Alternatively, we could also assume that it is the energy injection in (70), which is given by Taylor's expression,

$$p(t) = C \frac{k(t)^{3/2}}{L}. \qquad (75)$$

Carrying out the same analysis, we find then,

$$k_0(t) = \left( \frac{\epsilon(t)L}{C} \right)^{2/3}. \qquad (76)$$

and

$$k_1(t) = \frac{4}{9} \left(\frac{L}{C}\right)^{4/3} \frac{d\epsilon(t)}{dt} \epsilon(t)^{-2/3}. \qquad (77)$$

The arbitrariness in this approach to model either the production or the dissipation using Taylor's approach is removed by modeling rather the flux using Taylor's estimate. The disadvantage is that, to model the system using the introduction of the flux $\Pi(t)$, we need information on the decomposition, which in the simplest possible case, leads to the introduction of the additional parameter $\gamma$.

### C. Influence of the choice of the lengthscale

In experiments, the longitudinal correlation length $\Lambda$ is generally chosen as an estimate for $L$ in Taylor's expression (1). Indeed, for experimental measurements using hotwire anemometry this is the most convenient choice. For practical reasons an alternative definition uses the velocity correlations truncated at their first zero-crossing$^{60}$. Another definition involves the peak of the compensated wavenumber spectrum$^{39}$ and this seems to be an experimental equivalent of the choice we used in the present investigation.

In Fig. 6(a) we show how the equilibrium kinetic energy $k_0$ is changed when we use the time-dependent integral lengthscale $\Lambda(t)$, where

$$\Lambda(t) \sim k(t)^{-1} \int \kappa^{-1} E(\kappa, t) d\kappa, \qquad (78)$$
<figure><img src="image_8.png" /><figcaption>FIG. 6. (a) Temporal signals of the kinetic energy and estimates $k_0(t)$, comparing the equilibrium assuming either a constant $L$ indicated by $k_0(t)$ or a time-dependent one $\Lambda(t)$, indicated by $k_0^\Lambda(t)$. (b) Non-equilibrium estimate $k_1(t)$ and $k_1^\Lambda(t)$ using either a constant or a time-dependent integral lengthscale.</figcaption></figure>

instead of the fixed lengthscale $L$ to estimate the equilibrium kinetic energy. Qualitatively, the signals of the equilibrium estimates are close, with a similar timelag compared to the timesignal of the kinetic energy.

Following the derivations in sec. II, the non-equilibrium estimate can be obtained by a first-order perturbation of the equilibrium. The result is

$$k_1^\Lambda(t) = -\frac{4}{9} \frac{1}{(1+\gamma)} \left(\frac{\Lambda(t)}{C_\Pi'}\right)^{4/3} p(t)^{1/3} \left(\frac{\dot{p}(t)}{p(t)} + \frac{\dot{\Lambda}(t)}{\Lambda(t)}\right), \quad (79)$$
where $\dot{\Lambda} = d\Lambda/dt$ and $C'_{\Pi}$ is the flux-constant in

$$ \Pi(z, t) = C'_{\Pi} \frac{k(z, t)^{3/2}}{\Lambda(t)}. \qquad (80) $$

In Fig. 6(b) we show how this influences the estimate $k_1$. The agreement with the data does not substantially improve.

The last point to be discussed is the influence on the normalized dissipation rate. As illustrated in$^{24}$, using the integral-correlation lengthscale $\Lambda$ modifies the slope from a -3/2 powerlaw to a dependence close to a -15/14 dependence (see also$^{20,28}$). The derivation uses, however, the relation between the kinetic energy spectrum and the integral lengthscale to obtain a closed expression. There exists no closed exact expression to relate the integral lengthscale to either the energy injection or energy dissipation, so that we cannot derive a closed expression for $C_{\epsilon}(t)$ using $\Lambda(t)$ by prescribing the dissipation profile $\epsilon(t)$.

In the present investigation, where we tried to get around the use of spectral descriptions, we will not try to go beyond the Expressions (39),(41). Furthermore, in inhomogeneous turbulence the integral correlation lengthscale is not uniquely defined, since the flow is anisotropic so that different choices are possible for the definition of $\Lambda(z)$ (for instance based on correlations parallel or perpendicular to the mean-flow). We will leave the formulation of a satisfactory lengthscale equation and the study of its implication on the spatio-temporal dynamics of the normalized dissipation rate for further research.

## REFERENCES

¹G. I. Taylor. Statistical theory of turbulence. Proc. Roy. Soc. London. Ser. A, Math. Phys. Sci., 151:421-444, 1935.

²H. Tennekes and J.L. Lumley. A first course in turbulence. The MIT Press, 1972.

³Wouter J. T. Bos, Liang Shao, and J-P Bertoglio. Spectral imbalance and the normalized dissipation rate of turbulence. Physics of fluids, 19(4), 2007.

⁴K. R. Sreenivasan. An update on the energy dissipation rate in isotropic turbulence. Phys. Fluids, 10:528, 1998.

⁵Y. Kaneda, T. Ishihara, M. Yokokawa, K. Itakura, and A. Uno. Energy dissipation rate and energy spectrum in high resolution direct numerical simulations of turbulence in a periodic box. Phys. Fluids, 15(L21), 2003.
⁶W. D. McComb, A Berera, S. R. Yoffe, and M. F. Linkmann. Energy transfer and dissipation in forced isotropic turbulence. *Physical Review E*, 91(4):043013, 2015.

⁷P. C. Valente and J. C. Vassilicos. Universal dissipation scaling for nonequilibrium turbulence. *Phys. Rev. Lett.*, 108(21):214503, 2012.

⁸R. J. Hearst and P. Lavoie. Decay of turbulence generated by a square-fractal-element grid. *J. Fluid Mech.*, 741:567–584, 2014.

⁹S. Goto and J. C. Vassilicos. Energy dissipation and flux laws for unsteady turbulence. *Phys. Lett. A*, 379(16):1144–1148, 2015.

¹⁰J. C. Vassilicos. Dissipation in turbulent flows. *Ann. Rev. Fluid Mech.*, 47:95–114, 2015.

¹¹JG Chen, C Cuvier, J-M Foucaut, Y Ostovan, and JC Vassilicos. A turbulence dissipation inhomogeneity scaling in the wake of two side-by-side square prisms. *Journal of Fluid Mechanics*, 924, 2021.

¹²A. N. Kolmogorov. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. *Dokl. Akad. Nauk. SSSR*, 30:301, 1941.

¹³U. Frisch. *Turbulence, the legacy of A.N. Kolmogorov*. Cambridge University Press, 1995.

¹⁴A. Yoshizawa. Nonequilibrium effect of the turbulent-energy-production process on the inertial-range energy spectrum. *Phys. Rev. E*, 49(5):4065, 1994.

¹⁵S.L. Woodruff and R. Rubinstein. Multiple-scale perturbation analysis of slowly evolving turbulence. *J. Fluid Mech.*, 565:95, 2006.

¹⁶Robert Rubinstein and Timothy T Clark. “equilibrium” and “non-equilibrium” turbulence. *Theoretical and applied mechanics letters*, 7(5):301–305, 2017.

¹⁷K. Horiuti and T. Tamaki. Nonequilibrium energy spectrum in the subgrid-scale one-equation model in large-eddy simulation. *Phys. Fluids*, 25(12):125104, 2013.

¹⁸Kiyosi Horiuti and Tetsuya Ozawa. Multimode stretched spiral vortex and nonequilibrium energy spectrum in homogeneous shear flow turbulence. *Physics of Fluids*, 23(3), 2011.

¹⁹S Berti, G Boffetta, and S Musacchio. Mean flow and fluctuations in the three-dimensional turbulent cellular flow. *Physical Review Fluids*, 8(5):054601, 2023.

²⁰Le Fang and Wouter J. T. Bos. An edqnm study of the dissipation rate in isotropic non-equilibrium turbulence. *Journal of Turbulence*, pages 1–18, 2023.

²¹Ryo Araki and Wouter J. T. Bos. Inertial range scaling of inhomogeneous turbulence. *Journal of Fluid Mechanics*, 978:A9, 2024.

²²Wouter JT Bos. Unsteady and inhomogeneous turbulent fluctuations around isotropic
equilibrium. Atmosphere, 15(5):547, 2024.

<footnote n="23"><sup>23</sup>VJ Valadão, T Ceccotti, G Boffetta, and S Musacchio. Nonequilibrium fluctuations of the direct cascade in surface quasi-geostrophic turbulence. Physical Review Fluids, 9(9):094601, 2024.</footnote>

<footnote n="24"><sup>24</sup>Wouter Bos and Robert Rubinstein. Dissipation in unsteady turbulence. Physical Review Fluids, 2(2):022601, 2017.</footnote>

<footnote n="25"><sup>25</sup>Yulin Zheng, Kohtaro Nakamura, Koji Nagata, and Tomoaki Watanabe. Unsteady dissipation scaling in static-and active-grid turbulence. Journal of Fluid Mechanics, 956:A20, 2023.</footnote>

<footnote n="26"><sup>26</sup>Yulin Zheng, Naoya Koto, Koji Nagata, and Tomoaki Watanabe. Unsteady dissipation scaling of grid turbulence in the near-field region. Physics of Fluids, 35(9), 2023.</footnote>

<footnote n="27"><sup>27</sup>X Shao, J Fang, and L Fang. Non-equilibrium dissipation laws in a minimal two-scale wake model. Physics of Fluids, 35(8), 2023.</footnote>

<footnote n="28"><sup>28</sup>Marta Waclawczyk, Jakub L Nowak, Holger Siebert, and Szymon P Malinowski. Detecting nonequilibrium states in atmospheric turbulence. Journal of the Atmospheric Sciences, 79(10):2757-2772, 2022.</footnote>

<footnote n="29"><sup>29</sup>Wouter J. T. Bos. Production and dissipation of kinetic energy in grid turbulence. Phys. Rev. Fluids, 5:104607, 2020.</footnote>

<footnote n="30"><sup>30</sup>A. Laporta. Etude spectrale et modélisation de la turbulence inhomogène. PhD thesis, Ecole Centrale de Lyon, 1995.</footnote>

<footnote n="31"><sup>31</sup>DC Besnard, FH Harlow, RM Rauenzahn, and C Zemach. Spectral transport model for turbulence. Theoretical and computational fluid dynamics, 8(1):1-35, 1996.</footnote>

<footnote n="32"><sup>32</sup>Adrien Cahuzac, Jérôme Boudet, Pierre Borgnat, and Emmanuel Lévêque. Smoothing algorithms for mean-flow extraction in large-eddy simulation of complex turbulent flows. Physics of fluids, 22(12), 2010.</footnote>

<footnote n="33"><sup>33</sup>Yutaro Motoori and Susumu Goto. Hierarchy of coherent structures and real-space energy transfer in turbulent channel flow. Journal of Fluid Mechanics, 911:A27, 2021.</footnote>

<footnote n="34"><sup>34</sup>PC Valente, R Onishi, and CB da Silva. Origin of the imbalance between energy cascade and dissipation in turbulence. Physical Review E, 90(2):023003, 2014.</footnote>

<footnote n="35"><sup>35</sup>A Cadiou, K Hanjalić, and K Stawiarski. A two-scale second-moment turbulence closure based on weighted spectrum integration. Theoretical and Computational Fluid Dynamics, 18(1):1-26, 2004.</footnote>
<footnote n="36"><sup>36</sup>S.B. Pope. *Turbulent Flows*. Cambridge University Press, 2000.</footnote>

<footnote n="37"><sup>37</sup>A. Delache, C. Cambon, and F. Godeferd. Scale by scale anisotropy in freely decaying rotating turbulence. *Phys. Fluids*, 26(2):025104, 2014.</footnote>

<footnote n="38"><sup>38</sup>B. R. Pearson, P.-A. Krogstad, and W. van de Water. Measurements of the turbulent energy dissipation rate. *Phys. Fluids*, 14:1288, 2002.</footnote>

<footnote n="39"><sup>39</sup>P. Burattini, P. Lavoie, and R. A. Antonia. On the normalized turbulent energy dissipation rate. *Phys. Fluids*, 17:098103, 2005.</footnote>

<footnote n="40"><sup>40</sup>Daniel Odens Mora and Martin Obligado. Estimating the integral length scale on turbulent flows from the zero crossings of the longitudinal velocity fluctuation. *Experiments in fluids*, 61(9):199, 2020.</footnote>

<footnote n="41"><sup>41</sup>Wouter J. T. Bos, Timothy T Clark, and Robert Rubinstein. Small scale response and modeling of periodically forced turbulence. *Physics of Fluids*, 19(5), 2007.</footnote>

<footnote n="42"><sup>42</sup>A. K. Kuczaj, B. J. Geurts, and D. Lohse. Response maxima in time-modulated turbulence: Direct numerical simulations. *Europhys. Lett.*, 73:851, 2006.</footnote>

<footnote n="43"><sup>43</sup>Yuyao Yang, Robert Chahine, Robert Rubinstein, and Wouter J. T. Bos. Passive scalar mixing in modulated turbulence. *Fluid Dynamics Research*, 51(4):045501, 2019.</footnote>

<footnote n="44"><sup>44</sup>Vadim Borue and Steven A Orszag. Numerical study of three-dimensional kolmogorov flow at high reynolds numbers. *Journal of Fluid Mechanics*, 306:293-323, 1996.</footnote>

<footnote n="45"><sup>45</sup>Stefano Musacchio and Guido Boffetta. Turbulent channel without boundaries: The periodic kolmogorov flow. *Physical Review E*, 89(2):023004, 2014.</footnote>

<footnote n="46"><sup>46</sup>Wenwei Wu, François G Schmitt, Enrico Calzavarini, and Lipo Wang. A quadratic reynolds stress development for the turbulent kolmogorov flow. *Physics of Fluids*, 33(12):125129, 2021.</footnote>

<footnote n="47"><sup>47</sup>D. Lohse. Periodically kicked turbulence. *Phys. Rev. E*, 62:4946, 2000.</footnote>

<footnote n="48"><sup>48</sup>A. von der Heydt, S. Grossmann, and D. Lohse. Response maxima in modulated turbulence. *Phys. Rev. E*, 67:046308, 2003.</footnote>

<footnote n="49"><sup>49</sup>O. Cadot, J.H. Titon, and D. Bonn. Observation of resonances in modulated turbulence. *J. Fluid Mech.*, 485:161, 2003.</footnote>

<footnote n="50"><sup>50</sup>H.E. Cekli, C. Tipton, and W. van de Water. Resonant enhancement of turbulent energy dissipation. *Phys. Rev. Lett.*, 105(4):044503, 2010.</footnote>

<footnote n="51"><sup>51</sup>Lukas Bentkamp and Michael Wilczek. Temporal large-scale intermittency and its impact on the statistics of turbulence. *arXiv preprint arXiv:2312.08281*, 2023.</footnote>
<footnote n="52">⁵²Robert Rubinstein and Ye Zhou. Schiestel’s derivation of the epsilon equation and two-equation modelling of rotating turbulence. *Computers & Mathematics with Applications*, 46(4):633–638, 2003.</footnote>

<footnote n="53">⁵³R. Rubinstein, T.T. Clark, D. Livescu, and L.-S. Luo. Time-dependent isotropic turbulence. *J. Turbul.*, 5:011, 2004.</footnote>

<footnote n="54">⁵⁴R. Rubinstein and W. J. T. Bos. On the unsteady behavior of turbulence models. *Phys. Fluids*, 21:041701, 2009.</footnote>

<footnote n="55">⁵⁵W Peter Jones and Brian Edward Launder. The prediction of laminarization with a two-equation model of turbulence. *International journal of heat and mass transfer*, 15(2):301–314, 1972.</footnote>

<footnote n="56">⁵⁶Ryo Araki, Wouter JT Bos, and Susumu Goto. Minimal model of quasi-cyclic behaviour in turbulence driven by taylor-green forcing. *Fluid Dynamics Research*, 55(3):035507, 2023.</footnote>

<footnote n="57">⁵⁷Sébastien Aumaître and Stephan Fauve. Intermittency as a consequence of a stationarity constraint on the energy flux. *Physical Review Letters*, 132(11):114002, 2024.</footnote>

<footnote n="58">⁵⁸S Aumaître, J Farago, S Fauve, and S Mc Namara. Energy and power fluctuations in vibrated granular gases. *The European Physical Journal B-Condensed Matter and Complex Systems*, 42(2):255–261, 2004.</footnote>

<footnote n="59">⁵⁹F Schmitt, A Fuchs, J Peinke, and M Obligado. Refined universality of the turbulence cascade. *arXiv preprint arXiv:2407.15953*, 2024.</footnote>

<footnote n="60">⁶⁰F Alves Portela, G Papadakis, and JC Vassilicos. Turbulence dissipation and the role of coherent structures in the near wake of a square prism. *Physical Review Fluids*, 3(12):124609, 2018.</footnote>