# Inertial range scaling of inhomogeneous turbulence

Ryo Araki, Wouter J T Bos

► To cite this version:

Ryo Araki, Wouter J T Bos. Inertial range scaling of inhomogeneous turbulence. Journal of Fluid Mechanics, 2023, 978, 10.1017/jfm.2023.940. hal-04381979

HAL Id: hal-04381979

https://hal.science/hal-04381979

Submitted on 21 Jun 2024

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
Banner appropriate to article type will appear here in typeset article

# 1 Inertial range scaling of inhomogeneous turbulence

## 2 Ryo Araki¹,² and Wouter J. T. Bos¹†

¹Univ Lyon, CNRS, Ecole Centrale de Lyon, INSA Lyon, Univ Claude Bernard Lyon 1, LMFA,

UMR5509, 69340 Ecully, France

²Graduate School of Engineering Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka

560-8531, Japan

7 (Received xx; revised xx; accepted xx)

8 Abstract: We investigate how inhomogeneity influences the $k^{-5/3}$ inertial range scaling of turbulent kinetic energy spectra (with $k$ the wavenumber). For weak statistical inhomogeneity, the energy spectrum can be described as an equilibrium spectrum plus a perturbation.

10 Theoretical arguments suggest that this latter contribution scales as $k^{-7/3}$. This prediction is assessed using direct numerical simulations of three-dimensional Kolmogorov flow.

12

1. Introduction

13 Kolmogorov postulated in 1941 that the small scales of turbulent flows away from boundaries can be considered universal if the Reynolds number is sufficiently large (Kolmogorov 1941). The small scales are then supposed to be in equilibrium, and the energy spectrum satisfies,

$$E(k, \mathbf{x}, t) \sim \epsilon(\mathbf{x}, t)^{2/3} k^{-5/3}, \quad (1.1)$$

17 where $\epsilon$ is the average energy dissipation rate. For this expression to hold, the wavenumber $k$ should be sufficiently large compared to $L(\mathbf{x}, t)^{-1}$, the inverse of the length scale characterising the largest scales of the flow, and sufficiently small compared to the inverse of the Kolmogorov-scale $\eta(\mathbf{x}, t)^{-1}$ (with $\eta = v^{3/4}\epsilon^{-1/4}$), associated with the smallest scale of the flow.

21 In expression (1.1), the time and space dependence of $E$ and $\epsilon$ need some particular attention. Theoretically, the most convenient flow type for investigating inertial range scaling is an infinitely large and statistically stationary flow without boundaries. Since all practical flows are limited in size and lifetime, the dissipation rate will be dependent, even on average, on either position $\mathbf{x}$ or time $t$, or both. Expression (1.1) will therefore hold only locally in subdomains of space and time-intervals large enough compared to the considered length and time scales.

23 Indeed, the assumptions allowing the simple prediction (1.1) are that the scales $k$ can be considered locally isotropic, stationary, and homogeneous. The criterion $k \gg L^{-1}$ represents the implicit assumption that the influence of anisotropy, instationarity and inhomogeneity decreases as a function of scale. The present investigation aims to substantiate this assumption, particularly regarding the influence of inhomogeneity.

35 As an illustration, let us discuss the influence of statistical instationarity on the behaviour

<footnote>† Email address for correspondence: wouter.bos@ec-lyon.fr</footnote>
Banner appropriate to article type will appear here in typeset article

# 1 Inertial range scaling of inhomogeneous turbulence

## 2 Ryo Araki¹,² and Wouter J. T. Bos¹†

¹Univ Lyon, CNRS, Ecole Centrale de Lyon, INSA Lyon, Univ Claude Bernard Lyon 1, LMFA,

UMR5509, 69340 Ecully, France

²Graduate School of Engineering Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka

560-8531, Japan

7 (Received xx; revised xx; accepted xx)

8 Abstract: We investigate how inhomogeneity influences the $k^{-5/3}$ inertial range scaling of turbulent kinetic energy spectra (with $k$ the wavenumber). For weak statistical inhomogeneity, the energy spectrum can be described as an equilibrium spectrum plus a perturbation.

10 Theoretical arguments suggest that this latter contribution scales as $k^{-7/3}$. This prediction is assessed using direct numerical simulations of three-dimensional Kolmogorov flow.

## 1. Introduction

14 Kolmogorov postulated in 1941 that the small scales of turbulent flows away from boundaries can be considered universal if the Reynolds number is sufficiently large (Kolmogorov 1941).

16 The small scales are then supposed to be in equilibrium, and the energy spectrum satisfies,

$$E(k, \mathbf{x}, t) \sim \epsilon(\mathbf{x}, t)^{2/3} k^{-5/3}, \quad (1.1)$$

18 where $\epsilon$ is the average energy dissipation rate. For this expression to hold, the wavenumber $k$ should be sufficiently large compared to $L(\mathbf{x}, t)^{-1}$, the inverse of the length scale characterising the largest scales of the flow, and sufficiently small compared to the inverse of the Kolmogorov-scale $\eta(\mathbf{x}, t)^{-1}$ (with $\eta = v^{3/4} \epsilon^{-1/4}$), associated with the smallest scale of the flow.

20 In expression (1.1), the time and space dependence of $E$ and $\epsilon$ need some particular attention. Theoretically, the most convenient flow type for investigating inertial range scaling is an infinitely large and statistically stationary flow without boundaries. Since all practical flows are limited in size and lifetime, the dissipation rate will be dependent, even on average, on either position $\mathbf{x}$ or time $t$, or both. Expression (1.1) will therefore hold only locally in subdomains of space and time-intervals large enough compared to the considered length and time scales.

23 Indeed, the assumptions allowing the simple prediction (1.1) are that the scales $k$ can be considered locally isotropic, stationary, and homogeneous. The criterion $k \gg L^{-1}$ represents the implicit assumption that the influence of anisotropy, instationarity and inhomogeneity decreases as a function of scale. The present investigation aims to substantiate this assumption, particularly regarding the influence of inhomogeneity.

35 As an illustration, let us discuss the influence of statistical instationarity on the behaviour

<footnote>† Email address for correspondence: wouter.bos@ec-lyon.fr</footnote>
of the small scales. This subject was addressed by Yoshizawa (1994), who proposed that the influence of instationarity at large wavenumbers can be described as a perturbation on the energy spectrum as

$$E(k, \boldsymbol{x}, t) = E_0(k, \boldsymbol{x}, t) + E_1(k, \boldsymbol{x}, t), \qquad (1.2)$$

where the equilibrium part $E_0$ is given by (1.1) and the perturbation scales as

$$E_1^T(k, \boldsymbol{x}, t) = C_Y \frac{d\epsilon(\boldsymbol{x}, t)}{dt} \epsilon(\boldsymbol{x}, t)^{-2/3} k^{-7/3}, \qquad (1.3)$$

where the superscript T denotes that we consider perturbations due to instationarity. Numerical evidence of this scaling was first obtained by Horiuti & Ozawa (2011) for the case of homogeneous shear flow and by Horiuti & Tamaki (2013) for statistically isotropic turbulence in a periodic box. Further theoretical discussion and a more straightforward derivation of (1.3) can be found in Rubinstein & Clark (2005); Woodruff & Rubinstein (2006); Bos & Rubinstein (2017).

A similar approach is applied in various other configurations where the influence of external effects on isotropic turbulence is modelled as a perturbation to the energy spectrum. The effect of a mean-shear on isotropic turbulence was treated perturbatively by Ishihara et al. (2002). Stratified turbulence was considered in Kaneda & Yoshida (2004) and the limit of weakly compressible turbulence by Bertoglio et al. (2001).

The effect of large-scale temporal fluctuations on the kinetic energy spectrum is thus proportional to $k^{-7/3}$, which decays more rapidly than the equilibrium spectrum (1.1) with the $k^{-5/3}$ scaling. In the remainder of this investigation, we will focus on the influence of inhomogeneity on the scaling of turbulent kinetic energy, which has received little attention. Using Karhunen-Loeve eigenfunctions, it was illustrated that Kolmogorov's equilibrium spectrum can be observed in statistically inhomogeneous flows (Knight & Sirovich 1990; Moser 1994; Liao & Su 2015). By using the SO(3) symmetry group decomposition, Kurien et al. (2000) showed that structure functions contain a subdominant scaling component associated with inhomogeneity.

We further assess at which rate statistical homogeneity is recovered at small scales. To that, in § 2 we derive an analytical prediction of the scaling of $E_1^X(k, \boldsymbol{x})$, where the superscript X denotes the perturbation due to inhomogeneity, in stationary inhomogeneous turbulence. This expression will be the inhomogeneous equivalent of equation (1.3). In § 3, we report the results of Direct Numerical Simulations (DNS) of the three-dimensional Kolmogorov flow to assess the predictions. Section 4 concludes this investigation.

## 2. Derivation of the spectral correction due to inhomogeneity

The main difficulty in the present investigation comes from the fact that we investigate a multi-scale description (the energy spectrum) in an inhomogeneous setting. To simplify as far as possible, we restrict ourselves to a fairly simple setting, where the (statistical) inhomogeneity is periodic in space, and the flow is stationary and far away from boundaries. Before addressing the inhomogeneous multi-scale description, we will first consider the pointwise energy balance of the flow.

### 2.1. Kinetic-energy budget in inhomogeneous turbulence

We consider a statistically inhomogeneous flow kept in a statistically stationary state by a steady forcing $f(z)$. The forcing in the present manuscript consists of a unidirectional steady body force in the x-direction with a sinusoidal dependence in the z-direction. The
of the small scales. This subject was addressed by Yoshizawa (1994), who proposed that the influence of instationarity at large wavenumbers can be described as a perturbation on the energy spectrum as

$$E(k, \boldsymbol{x}, t) = E_0(k, \boldsymbol{x}, t) + E_1(k, \boldsymbol{x}, t), \qquad (1.2)$$

where the equilibrium part $E_0$ is given by (1.1) and the perturbation scales as

$$E_1^T(k, \boldsymbol{x}, t) = C_Y \frac{d\epsilon(\boldsymbol{x}, t)}{dt} \epsilon(\boldsymbol{x}, t)^{-2/3} k^{-7/3}, \qquad (1.3)$$

where the superscript T denotes that we consider perturbations due to instationarity. Numerical evidence of this scaling was first obtained by Horiuti & Ozawa (2011) for the case of homogeneous shear flow and by Horiuti & Tamaki (2013) for statistically isotropic turbulence in a periodic box. Further theoretical discussion and a more straightforward derivation of (1.3) can be found in Rubinstein & Clark (2005); Woodruff & Rubinstein (2006); Bos & Rubinstein (2017).

A similar approach is applied in various other configurations where the influence of external effects on isotropic turbulence is modelled as a perturbation to the energy spectrum. The effect of a mean-shear on isotropic turbulence was treated perturbatively by Ishihara et al. (2002). Stratified turbulence was considered in Kaneda & Yoshida (2004) and the limit of weakly compressible turbulence by Bertoglio et al. (2001).

The effect of large-scale temporal fluctuations on the kinetic energy spectrum is thus proportional to $k^{-7/3}$, which decays more rapidly than the equilibrium spectrum (1.1) with the $k^{-5/3}$ scaling. In the remainder of this investigation, we will focus on the influence of inhomogeneity on the scaling of turbulent kinetic energy, which has received little attention. Using Karhunen-Loeve eigenfunctions, it was illustrated that Kolmogorov's equilibrium spectrum can be observed in statistically inhomogeneous flows (Knight & Sirovich 1990; Moser 1994; Liao & Su 2015). By using the SO(3) symmetry group decomposition, Kurien et al. (2000) showed that structure functions contain a subdominant scaling component associated with inhomogeneity.

We further assess at which rate statistical homogeneity is recovered at small scales. To that, in § 2 we derive an analytical prediction of the scaling of $E_1^X(k, \boldsymbol{x})$, where the superscript X denotes the perturbation due to inhomogeneity, in stationary inhomogeneous turbulence. This expression will be the inhomogeneous equivalent of equation (1.3). In § 3, we report the results of Direct Numerical Simulations (DNS) of the three-dimensional Kolmogorov flow to assess the predictions. Section 4 concludes this investigation.

## 2. Derivation of the spectral correction due to inhomogeneity

The main difficulty in the present investigation comes from the fact that we investigate a multi-scale description (the energy spectrum) in an inhomogeneous setting. To simplify as far as possible, we restrict ourselves to a fairly simple setting, where the (statistical) inhomogeneity is periodic in space, and the flow is stationary and far away from boundaries. Before addressing the inhomogeneous multi-scale description, we will first consider the pointwise energy balance of the flow.

### 2.1. Kinetic-energy budget in inhomogeneous turbulence

We consider a statistically inhomogeneous flow kept in a statistically stationary state by a steady forcing $f(z)$. The forcing in the present manuscript consists of a unidirectional steady body force in the x-direction with a sinusoidal dependence in the z-direction. The
79 Navier-Stokes equations for this specific system write

$$ \frac{D\mathcal{U}(\mathbf{x}, t)}{Dt} = -\nabla\mathcal{P}(\mathbf{x}, t) + \nu\Delta\mathcal{U}(\mathbf{x}, t) + f(z)e_x, \quad (2.1) $$

81 where $D/Dt$ is the material derivative, $\mathcal{P}$ is the pressure (divided by density) ensuring incompressibility $\nabla \cdot \mathcal{U} = 0$, and $e_x$ denotes the unit vector in the x-direction.

83 The equations for the mean flow and the kinetic energy of the fluctuations can be derived by introducing the Reynolds decomposition $\mathcal{U} = \langle \mathcal{U} \rangle + \mathbf{u}$, where $\langle \mathcal{U} \rangle$ is the ensemble-averaged velocity and $\mathbf{u} = (u, v, w)$ the fluctuation. The specific forcing considered in the present investigation leads to a mean flow $\langle \mathcal{U}(\mathbf{x}, t) \rangle = U(z)e_x$. Then, the kinetic energy corresponding to the mean flow can be written as

$$ K_U(z) = U(z)^2/2, \quad (2.2) $$

89 and the kinetic energy of the fluctuations is

$$ K(z) = \frac{1}{2} [\langle u^2 \rangle(z) + \langle v^2 \rangle(z) + \langle w^2 \rangle(z)]. \quad (2.3) $$

91 The equation for the mean-velocity $U(z)$ reduces to,

$$ \frac{D U(z)}{Dt} = -\frac{\partial}{\partial z} \langle uw \rangle(z) + f(z) + \nu \frac{\partial^2 U(z)}{\partial z^2} = 0. \quad (2.4) $$

93 The details are, for instance, provided in Bos (2020). The equation for the turbulent kinetic energy writes, in a steady state,

$$ \frac{DK(z)}{Dt} = p(z) - \epsilon(z) + d(z) = 0, \quad (2.5) $$

95 where the production $p(z)$, dissipation $\epsilon(z)$, and diffusion $d(z)$ terms are given by

$$ p(z) = -\langle uw \rangle(z) \frac{\partial U(z)}{\partial z}, \quad (2.6) $$

$$ \epsilon(z) = \nu \left\langle \frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j} \right\rangle(z), \quad (2.7) $$

$$ d(z) = -\frac{\partial}{\partial z} \left( \langle \mathcal{P}w \rangle(z) + \langle u_i u_i w \rangle(z) - \nu \frac{\partial K(z)}{\partial z} \right), \quad (2.8) $$

96 respectively. The first term $p(z)$ represents the production of turbulent kinetic energy through the interaction of the turbulent fluctuations with the mean-velocity gradient $\partial U(z)/\partial z$. The viscous dissipation term $\epsilon(z)$ involves the gradients of the fluctuating velocity.

99 In statistically homogeneous flows, production and dissipation are the only terms appearing in the turbulent kinetic energy balance. In statistically inhomogeneous flows, we also have spatial diffusion of turbulent kinetic energy $d(z)$. The diffusion contains contributions associated with the turbulent fluctuations of the velocity and pressure (first two terms) and a contribution through viscous diffusion (the last term). This viscous part of the diffusion is generally negligible compared to the contribution of the other two terms and will be dropped in the following.

106 The main question in the present investigation is how such inhomogeneous redistribution processes $d(z)$ affect the scaling of the kinetic energy spectrum $E(k, \mathbf{x})$ in the inertial range of high Reynolds number turbulence.
79 Navier-Stokes equations for this specific system write

$$ \frac{D\mathcal{U}(x,t)}{Dt} = -\nabla\mathcal{P}(x,t) + v\Delta\mathcal{U}(x,t) + f(z)e_x, \quad (2.1) $$

81 where $D/Dt$ is the material derivative, $\mathcal{P}$ is the pressure (divided by density) ensuring incompressibility $\nabla \cdot \mathcal{U} = 0$, and $e_x$ denotes the unit vector in the x-direction.

83 The equations for the mean flow and the kinetic energy of the fluctuations can be derived by introducing the Reynolds decomposition $\mathcal{U} = \langle \mathcal{U} \rangle + \mathbf{u}$, where $\langle \mathcal{U} \rangle$ is the ensemble-averaged velocity and $\mathbf{u} = (u, v, w)$ the fluctuation. The specific forcing considered in the present investigation leads to a mean flow $\langle \mathcal{U}(x,t) \rangle = U(z)e_x$. Then, the kinetic energy corresponding to the mean flow can be written as

$$ K_U(z) = U(z)^2/2, \quad (2.2) $$

89 and the kinetic energy of the fluctuations is

$$ K(z) = \frac{1}{2} [ \langle u^2 \rangle(z) + \langle v^2 \rangle(z) + \langle w^2 \rangle(z) ]. \quad (2.3) $$

91 The equation for the mean-velocity $U(z)$ reduces to,

$$ \frac{D U(z)}{Dt} = - \frac{\partial}{\partial z} \langle uw \rangle(z) + f(z) + v \frac{\partial^2 U(z)}{\partial z^2} = 0. \quad (2.4) $$

93 The details are, for instance, provided in Bos (2020). The equation for the turbulent kinetic energy writes, in a steady state,

$$ \frac{DK(z)}{Dt} = p(z) - \epsilon(z) + d(z) = 0, \quad (2.5) $$

95 where the production $p(z)$, dissipation $\epsilon(z)$, and diffusion $d(z)$ terms are given by

$$ p(z) = -\langle uw \rangle(z) \frac{\partial U(z)}{\partial z}, \quad (2.6) $$

$$ \epsilon(z) = v \left\langle \frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j} \right\rangle(z), \quad (2.7) $$

$$ d(z) = -\frac{\partial}{\partial z} \left( \langle \mathcal{P}w \rangle(z) + \langle u_i u_i w \rangle(z) - v \frac{\partial K(z)}{\partial z} \right), \quad (2.8) $$

96 respectively. The first term $p(z)$ represents the production of turbulent kinetic energy through the interaction of the turbulent fluctuations with the mean-velocity gradient $\partial U(z)/\partial z$. The viscous dissipation term $\epsilon(z)$ involves the gradients of the fluctuating velocity.

99 In statistically homogeneous flows, production and dissipation are the only terms appearing in the turbulent kinetic energy balance. In statistically inhomogeneous flows, we also have spatial diffusion of turbulent kinetic energy $d(z)$. The diffusion contains contributions associated with the turbulent fluctuations of the velocity and pressure (first two terms) and a contribution through viscous diffusion (the last term). This viscous part of the diffusion is generally negligible compared to the contribution of the other two terms and will be dropped in the following.

106 The main question in the present investigation is how such inhomogeneous redistribution processes $d(z)$ affect the scaling of the kinetic energy spectrum $E(k, \mathbf{x})$ in the inertial range of high Reynolds number turbulence.
## 2.2. Fourier-analysis of inhomogeneous turbulence

109

The use of energy spectra in general turbulent flows needs some justification. In principle, Fourier modes are associated with infinite or periodic domains. This property would exclude the use of spatial Fourier analysis of any realistic, non-periodic flow. However, a closer look at the lengthscales involved in turbulent flows permits invoking an assumption of scale separation, allowing us to get around this problem. Indeed, the theoretical basis for practical Fourier modelling of non-periodic turbulent flows can be found in various works (see Jeandel et al. 1978; Yoshizawa 1984; Bertoglio & Jeandel 1987; Laporta & Bertoglio 1995; Besnard et al. 1996). In practice, to develop a spectral description of inhomogeneous flows, one needs to introduce a lengthscale L characterising the inhomogeneity of the flow geometry. Then, one can consider Fourier spectra associated with scales $r \sim k^{-1}$ small compared to L.

110

111

112

113

114

115

116

117

118

119

In the present investigation, we consider a spatially periodic flow without solid boundaries or obstacles to avoid most of these complications. Furthermore, to derive corrections due to statistical inhomogeneity, we consider statistically stationary turbulence with a single inhomogeneous direction $z$. An advantage of the present configuration, where only one inhomogeneous direction is present, is that we can compute energy spectra in planes perpendicular to the z-axis. We thus define

120

121

122

123

124

125

$$E(k_{\perp}, z) = \frac{1}{2} \int u_i(k_{\perp}, z) u_i^*(k_{\perp}, z) \, dA(k_{\perp}), \qquad (2.9)$$

126

where $A(k_{\perp})$ denotes a wavenumber-shell of radius $k_{\perp}$ in the $k_x, k_y$ plane. The velocity field in (2.9) is defined by the two-dimensional Fourier transform,

127

128

$$u_i(k_{\perp}, z) = \int e^{-i(k_x x + k_y y)} u_i(x, y, z) \, dx \, dy. \qquad (2.10)$$

129

The resulting energy spectrum $E(k_{\perp}, z)$ is a function of a perpendicular wavenumber $k_{\perp} = \sqrt{k_x^2 + k_y^2}$ and a vertical coordinate $z$. We note that if isotropy is restored in small scales, $E(k_{\perp}, z)$ is expected to scale like the three-dimensional spectrum $E(k, z)$ (see Appendix A for the definition). In the following subsections in § 2, we will keep the notation $E(k, z)$ for the sake of generality, but it should be kept in mind that the scaling of $E(k, z)$ and $E(k_{\perp}, z)$ should be equivalent in statistically isotropic flow at large $k$.

130

131

132

133

134

135

136

## 2.3. Governing equation and modelling

136

The derivation in this subsection closely follows the rationale used to derive the instationary correction presented in Bos & Rubinstein (2017). This same methodology is here applied to the evolution-equation of the energy spectrum in inhomogeneous turbulence.

137

138

139

The kinetic energy spectrum is associated with the turbulent kinetic energy by the relation

140

$$\int E(k, z) \, dk = K(z). \qquad (2.11)$$

141

The evolution equation for $E(k, z)$ is the multi-scale extension of equation (2.5). This equation reads, for the case of a unidirectional mean flow $U(z)e_x$ as in (2.1),

142

143

$$\frac{DE(k, z)}{Dt} = \underbrace{P(k, z)}_{\text{production}} - \underbrace{2vk^2E(k, z)}_{\text{dissipation}} + \underbrace{T(k, z)}_{\text{transfer}} + \underbrace{D(k, z)}_{\text{diffusion}}. \qquad (2.12)$$

144

For self-consistency, we discuss the derivation of this equation in Appendix A. Except for the viscous dissipation, all the terms in (2.12) are unclosed. In the following, we discuss the different physics and contributions to propose simple models for them.

145

146

147
## 2.2. Fourier-analysis of inhomogeneous turbulence

The use of energy spectra in general turbulent flows needs some justification. In principle, Fourier modes are associated with infinite or periodic domains. This property would exclude the use of spatial Fourier analysis of any realistic, non-periodic flow. However, a closer look at the lengthscales involved in turbulent flows permits invoking an assumption of scale separation, allowing us to get around this problem. Indeed, the theoretical basis for practical Fourier modelling of non-periodic turbulent flows can be found in various works (see Jeandel et al. 1978; Yoshizawa 1984; Bertoglio & Jeandel 1987; Laporta & Bertoglio 1995; Besnard et al. 1996). In practice, to develop a spectral description of inhomogeneous flows, one needs to introduce a lengthscale L characterising the inhomogeneity of the flow geometry. Then, one can consider Fourier spectra associated with scales $r \sim k^{-1}$ small compared to L.

In the present investigation, we consider a spatially periodic flow without solid boundaries or obstacles to avoid most of these complications. Furthermore, to derive corrections due to statistical inhomogeneity, we consider statistically stationary turbulence with a single inhomogeneous direction z. An advantage of the present configuration, where only one inhomogeneous direction is present, is that we can compute energy spectra in planes perpendicular to the z-axis. We thus define

$$E(k_{\perp}, z) = \frac{1}{2} \int u_i(k_{\perp}, z) u_i^*(k_{\perp}, z) dA(k_{\perp}), \quad (2.9)$$

where $A(k_{\perp})$ denotes a wavenumber-shell of radius $k_{\perp}$ in the $k_x, k_y$ plane. The velocity field in (2.9) is defined by the two-dimensional Fourier transform,

$$u_i(k_{\perp}, z) = \int e^{-i(k_x x + k_y y)} u_i(x, y, z) \, dx \, dy. \quad (2.10)$$

The resulting energy spectrum $E(k_{\perp}, z)$ is a function of a perpendicular wavenumber $k_{\perp} = \sqrt{k_x^2 + k_y^2}$ and a vertical coordinate $z$. We note that if isotropy is restored in small scales, $E(k_{\perp}, z)$ is expected to scale like the three-dimensional spectrum $E(k, z)$ (see Appendix A for the definition). In the following subsections in § 2, we will keep the notation $E(k, z)$ for the sake of generality, but it should be kept in mind that the scaling of $E(k, z)$ and $E(k_{\perp}, z)$ should be equivalent in statistically isotropic flow at large $k$.

## 2.3. Governing equation and modelling

The derivation in this subsection closely follows the rationale used to derive the instationary correction presented in Bos & Rubinstein (2017). This same methodology is here applied to the evolution-equation of the energy spectrum in inhomogeneous turbulence.

The kinetic energy spectrum is associated with the turbulent kinetic energy by the relation

$$\int E(k, z) dk = K(z). \quad (2.11)$$

The evolution equation for $E(k, z)$ is the multi-scale extension of equation (2.5). This equation reads, for the case of a unidirectional mean flow $U(z)e_x$ as in (2.1),

$$\frac{DE(k, z)}{Dt} = \underbrace{P(k, z)}_{\text{production}} - \underbrace{2vk^2E(k, z)}_{\text{dissipation}} + \underbrace{T(k, z)}_{\text{transfer}} + \underbrace{D(k, z)}_{\text{diffusion}}. \quad (2.12)$$

For self-consistency, we discuss the derivation of this equation in Appendix A. Except for the viscous dissipation, all the terms in (2.12) are unclosed. In the following, we discuss the different physics and contributions to propose simple models for them.
<figure><img src="image_1.png" /><figcaption>Figure 1: A schematic of the energy spectrum in (k, z) coordinates. Two arrows denote the direction of energy fluxes in wavenumber and physical space, respectively.</figcaption></figure>

Since the flow is statistically stationary and the mean flow is unidirectional, the material derivative on the left-hand side of (2.12) is zero. The first term on the RHS, $P(k, z)$, represents the terms directly proportional to the mean-velocity gradient. It contains two contributions: the production of turbulent kinetic energy and a linear transfer term (Cambon et al. 1981; Briard et al. 2018). These terms are mainly important at large scales and become zero at points in space where the velocity gradient vanishes. The order of magnitude of the production term can be estimated by (Tennekes & Lumley 1972),

$$P(k, z) \sim \left( \frac{\partial U(z)}{\partial z} \right)^2 \tau(k, z) E(k, z), \quad (2.13)$$

with the time scale $\tau(k, z) \sim \epsilon(z)^{-1/3}k^{-2/3}$ in the inertial range. The integral of $P(k, z)$ over wavenumbers yields $p(z)$ in (2.5). Here, $\epsilon(z)$ denotes the profile of the dissipation of kinetic energy through viscous stresses (see (2.7)) and is obtained by the integral of the second term on the RHS of (2.12). At large Reynolds numbers, this term is significant only at large wavenumbers. It is the term $P(k, z)$ which is responsible for energy transfer between the mean velocity field $U(z)$ and the turbulent kinetic energy.

The nonlinear transfer $T(k, z)$ represents the energy flux and is a redistributive term in scale space; thus, its integral over all wavenumbers yields zero. The last term $D(k, z)$ represents the diffusion, or transport, through turbulent fluctuations and viscous diffusion. Note that this term is zero in statistically homogeneous turbulence. The term $D(k, z)$ is also a redistribution term like $T(k, z)$, but in physical space. Its integral over wavenumbers corresponds to $d(z)$ in (2.5).

Both $T(k, z)$ and $D(k, z)$ are a function of triple correlations between Fourier modes at different wavelengths. There is no exact expression of these quantities as a closed function of the kinetic energy spectrum $E(k, z)$. At this moment, we will therefore introduce modelling assumptions. Sophisticated models exist for inhomogeneous spectral dynamics, based on the Test Field Model (Kraichnan 1971) or the Eddy-Damped Quasi-Normal Approximation (Laporta & Bertoglio 1995; Parpais & Bertoglio 1996). However, the resulting closures are quite complicated and do not allow a straightforward analytical perturbation treatment. Therefore, our approach uses simple models that reproduce their main physical features: the redistribution of energy in scale space for $T(k, z)$ and in physical
<figure><img src="image_2.png" /><figcaption>Figure 1: A schematic of the energy spectrum in (k, z) coordinates. Two arrows denote the direction of energy fluxes in wavenumber and physical space, respectively.</figcaption></figure>

Since the flow is statistically stationary and the mean flow is unidirectional, the material derivative on the left-hand side of (2.12) is zero. The first term on the RHS, $P(k, z)$, represents the terms directly proportional to the mean-velocity gradient. It contains two contributions: the production of turbulent kinetic energy and a linear transfer term (Cambon et al. 1981; Briard et al. 2018). These terms are mainly important at large scales and become zero at points in space where the velocity gradient vanishes. The order of magnitude of the production term can be estimated by (Tennekes & Lumley 1972),

$$P(k, z) \sim \left( \frac{\partial U(z)}{\partial z} \right)^2 \tau(k, z) E(k, z), \qquad (2.13)$$

with the time scale $\tau(k, z) \sim \epsilon(z)^{-1/3}k^{-2/3}$ in the inertial range. The integral of $P(k, z)$ over wavenumbers yields $p(z)$ in (2.5). Here, $\epsilon(z)$ denotes the profile of the dissipation of kinetic energy through viscous stresses (see (2.7)) and is obtained by the integral of the second term on the RHS of (2.12). At large Reynolds numbers, this term is significant only at large wavenumbers. It is the term $P(k, z)$ which is responsible for energy transfer between the mean velocity field $U(z)$ and the turbulent kinetic energy.

The nonlinear transfer $T(k, z)$ represents the energy flux and is a redistributive term in scale space; thus, its integral over all wavenumbers yields zero. The last term $D(k, z)$ represents the diffusion, or transport, through turbulent fluctuations and viscous diffusion. Note that this term is zero in statistically homogeneous turbulence. The term $D(k, z)$ is also a redistribution term like $T(k, z)$, but in physical space. Its integral over wavenumbers corresponds to $d(z)$ in (2.5).

Both $T(k, z)$ and $D(k, z)$ are a function of triple correlations between Fourier modes at different wavelengths. There is no exact expression of these quantities as a closed function of the kinetic energy spectrum $E(k, z)$. At this moment, we will therefore introduce modelling assumptions. Sophisticated models exist for inhomogeneous spectral dynamics, based on the Test Field Model (Kraichnan 1971) or the Eddy-Damped Quasi-Normal Approximation (Laporta & Bertoglio 1995; Parpais & Bertoglio 1996). However, the resulting closures are quite complicated and do not allow a straightforward analytical perturbation treatment. Therefore, our approach uses simple models that reproduce their main physical features: the redistribution of energy in scale space for $T(k, z)$ and in physical
space for $D(k, z)$, respectively. We use diffusion approximations for both terms,

$$T(k, z) = -\frac{\partial}{\partial k}\Pi(k, z), \qquad (2.14)$$

$$D(k, z) = -\frac{\partial}{\partial z}\Phi(k, z), \qquad (2.15)$$

where $\Pi(k, z)$ and $\Phi(k, z)$ are turbulent fluxes in wavenumber and physical space, respectively. Figure 1 schematically depicts these two fluxes in $(k, z)$ space. In the absence of inhomogeneity, the flux $\Phi(k, z)$ is zero. In the inhomogeneous case, the presence of this flux will affect the kinetic energy spectrum $E(k, z)$.

We model both fluxes using a gradient-diffusion approximation,

$$\Pi(k, z) = -\rho(k, z) \frac{\partial(k^{-2}E(k, z))}{\partial k} \qquad (2.16)$$

with $\rho(k, z) \sim k^{11/2}E(k, z)^{1/2}$ being a turbulent energy diffusion in Fourier space, and

$$\Phi(k, z) = -\mu(z) \frac{\partial E(k, z)}{\partial z}, \qquad (2.17)$$

where $\mu(z)$ is a turbulent diffusivity in real space (see (2.29)). We have effectively decoupled (and simplified) the transfer terms in scale and physical space. Indeed, both $\Phi(k, z)$ and $\Pi(k, z)$ are determined by the same triple velocity and velocity-pressure correlations (see (A 4)). Decomposing the physical space-scale space flux is a major assumption which seems necessary to obtain an analytically tractable model of energy transfer in inhomogeneous turbulence. The model for $\Pi(k, z)$ (2.16) is known as the Leith model (Leith 1967; Rubinstein & Clark 2022). This model tends to homogenise the kinetic energy in spectral space towards equipartition among wave vectors, corresponding to an energy spectrum proportional to $k^2$. The gradient-diffusion model for the diffusion (2.17) tends to homogenise the energy distribution in physical space and is used in Besnard et al. (1996); Touil et al. (2002); Cadiou et al. (2004), for instance.

Eddy viscosity models are obviously simplified representations of the real transfer terms. For instance, see Pope (2000, § 10) for extensive discussion. However, we think that this kind of modelling is a useful first step before turning to more sophisticated modelling approaches.

## 2.4. Linear perturbation analysis and scaling predictions

Our goal is to derive a prediction for inertial range scaling at large Reynolds numbers in the limit of weak inhomogeneity, where the influence of inhomogeneity can be treated as a perturbation. In the following, the leading order contributions and perturbations are indicated by a subscript 0 and 1, respectively. We define an inertial range $L^{-1} \ll k \ll \eta^{-1}$ with the length $L$ representing the typical length of the largest and energy-containing scales of the flow. Furthermore, in our description, it is associated with the longest wavelength in our flow domain and is chosen constant. We will define this length scale more precisely later, in § 2.5.

We now define the equilibrium about which we expand the equations. To do so, we consider the decomposition

$$E(k, z) = E_0(k, z) + E_1(k, z) \qquad (2.18)$$

with $|E_1| \ll |E_0|$. The other quantities, such as $\Pi(k, z)$ and $\Phi(k, z)$, are decomposed in the same manner. We recall here that in addition to these two contributions to the energy spectrum, the flow also contains the time-averaged velocity profile, which consists of a single wave vector in the $z$-direction in the present case (2.2). This mean flow is not present in the
space for $D(k, z)$, respectively. We use diffusion approximations for both terms,

$$T(k, z) = -\frac{\partial}{\partial k}\Pi(k, z), \qquad (2.14)$$

$$D(k, z) = -\frac{\partial}{\partial z}\Phi(k, z), \qquad (2.15)$$

where $\Pi(k, z)$ and $\Phi(k, z)$ are turbulent fluxes in wavenumber and physical space, respectively. Figure 1 schematically depicts these two fluxes in $(k, z)$ space. In the absence of inhomogeneity, the flux $\Phi(k, z)$ is zero. In the inhomogeneous case, the presence of this flux will affect the kinetic energy spectrum $E(k, z)$.

We model both fluxes using a gradient-diffusion approximation,

$$\Pi(k, z) = -\rho(k, z) \frac{\partial(k^{-2}E(k, z))}{\partial k} \qquad (2.16)$$

with $\rho(k, z) \sim k^{11/2}E(k, z)^{1/2}$ being a turbulent energy diffusion in Fourier space, and

$$\Phi(k, z) = -\mu(z) \frac{\partial E(k, z)}{\partial z}, \qquad (2.17)$$

where $\mu(z)$ is a turbulent diffusivity in real space (see (2.29)). We have effectively decoupled (and simplified) the transfer terms in scale and physical space. Indeed, both $\Phi(k, z)$ and $\Pi(k, z)$ are determined by the same triple velocity and velocity-pressure correlations (see (A 4)). Decomposing the physical space-scale space flux is a major assumption which seems necessary to obtain an analytically tractable model of energy transfer in inhomogeneous turbulence. The model for $\Pi(k, z)$ (2.16) is known as the Leith model (Leith 1967; Rubinstein & Clark 2022). This model tends to homogenise the kinetic energy in spectral space towards equipartition among wave vectors, corresponding to an energy spectrum proportional to $k^2$. The gradient-diffusion model for the diffusion (2.17) tends to homogenise the energy distribution in physical space and is used in Besnard et al. (1996); Touil et al. (2002); Cadiou et al. (2004), for instance.

Eddy viscosity models are obviously simplified representations of the real transfer terms. For instance, see Pope (2000, § 10) for extensive discussion. However, we think that this kind of modelling is a useful first step before turning to more sophisticated modelling approaches.

## 2.4. Linear perturbation analysis and scaling predictions

Our goal is to derive a prediction for inertial range scaling at large Reynolds numbers in the limit of weak inhomogeneity, where the influence of inhomogeneity can be treated as a perturbation. In the following, the leading order contributions and perturbations are indicated by a subscript 0 and 1, respectively. We define an inertial range $L^{-1} \ll k \ll \eta^{-1}$ with the length $L$ representing the typical length of the largest and energy-containing scales of the flow. Furthermore, in our description, it is associated with the longest wavelength in our flow domain and is chosen constant. We will define this length scale more precisely later, in § 2.5.

We now define the equilibrium about which we expand the equations. To do so, we consider the decomposition

$$E(k, z) = E_0(k, z) + E_1(k, z) \qquad (2.18)$$

with $|E_1| \ll |E_0|$. The other quantities, such as $\Pi(k, z)$ and $\Phi(k, z)$, are decomposed in the same manner. We recall here that in addition to these two contributions to the energy spectrum, the flow also contains the time-averaged velocity profile, which consists of a single wave vector in the $z$-direction in the present case (2.2). This mean flow is not present in the
inertial range, on which we will focus in the following. Therefore, in the remainder of this section, we can focus on the contributions $E_0(k, z)$ and $E_1(k, z)$.

For very high Reynolds numbers in the limit of vanishing inhomogeneity, we assume that the equilibrium contributions to the kinetic energy balance (2.12) do not depend on the inhomogeneous turbulent diffusion $D(k, z)$. By integrating the balance between the transfer and dissipation terms in (2.12) from $k$ to $\infty$, we find

$$ \int_k^\infty T(p, z) dp = \int_k^\infty 2\nu p^2 E(p, z) dp \quad (2.19) $$

or, using expression (2.14) and the equilibrium/nonequilibrium decomposition,

$$ \Pi_0(k, z) = \epsilon(z). \quad (2.20) $$

Indeed, this corresponds to the equilibrium between the energy flux and the energy dissipation rate, essential to the inertial range description of Kolmogorov (1941). The constant flux solution of the Leith model is consistent with this framework and is given by

$$ E_0(k, z) \sim \epsilon(z)^{2/3} k^{-5/3}. \quad (2.21) $$

This expression defines our equilibrium solution. We now assess the influence of the inhomogeneity of $\epsilon(z)$ on this scaling as a perturbation.

In the following, we consider the terms in the balance equation (2.12) for the nonequilibrium contributions. The order of magnitude of the production term (2.13) and the diffusion-gradient modelling with the flux (2.17) leads us to deduce that $D(k, z) \gg P(k, z)$ at $k \gg L^{-1}$. Therefore, the first-order perturbation to the equilibrium scaling in the inertial range is due to the inhomogeneous diffusion $D(k, z)$. Then, in the inertial range, we have

$$ T(k, z) = -D(k, z), \quad (2.22) $$

and

$$ -\frac{\partial}{\partial k} \Pi_1(k, z) = \frac{\partial}{\partial z} \Phi_0(k, z), \quad (2.23) $$

since $\partial\Pi_0(k, z)/\partial k = 0$. Thus, the first-order correction of the nonlinear transfer balances the zeroth-order contribution of the inhomogeneous diffusion. The first-order perturbation to the nonlinear flux $\Pi_1(k, z)$ is evaluated as (Rubinstein & Clark 2005)

$$ \Pi_1(k, z) = E_1(k, z) \left. \frac{\delta \Pi}{\delta E} \right|_{E_0}, \quad (2.24) $$

where $\delta\Pi/\delta E|_{E_0}$ is the Fréchet derivative of the total flux $\Pi$ evaluated at $E(k, z) = E_0(k, z)$.

In the inertial range, assuming $E_1$ to scale as a power law, this yields the scaling,

$$ \Pi_1(k, z) \sim \epsilon(z) \frac{E_1(k, z)}{E_0(k, z)}. \quad (2.25) $$

Note that we obtain (2.25) not only for the Leith model, but also for most of the other classical closures such as the Kovaznay and Heisenberg model (Rubinstein & Clark 2022).

Integrating (2.23) from $k$ to $\infty$, we have

$$ \Pi_1(k, z) = \frac{\partial}{\partial z} \int_k^\infty \Phi_0(k, z) dk. \quad (2.26) $$
inertial range, on which we will focus in the following. Therefore, in the remainder of this section, we can focus on the contributions $E_0(k, z)$ and $E_1(k, z)$.

For very high Reynolds numbers in the limit of vanishing inhomogeneity, we assume that the equilibrium contributions to the kinetic energy balance (2.12) do not depend on the inhomogeneous turbulent diffusion $D(k, z)$. By integrating the balance between the transfer and dissipation terms in (2.12) from $k$ to $\infty$, we find

$$ \int_k^\infty T(p, z) dp = \int_k^\infty 2\nu p^2 E(p, z) dp \quad (2.19) $$

or, using expression (2.14) and the equilibrium/nonequilibrium decomposition,

$$ \Pi_0(k, z) = \epsilon(z). \quad (2.20) $$

Indeed, this corresponds to the equilibrium between the energy flux and the energy dissipation rate, essential to the inertial range description of Kolmogorov (1941). The constant flux solution of the Leith model is consistent with this framework and is given by

$$ E_0(k, z) \sim \epsilon(z)^{2/3} k^{-5/3}. \quad (2.21) $$

This expression defines our equilibrium solution. We now assess the influence of the inhomogeneity of $\epsilon(z)$ on this scaling as a perturbation.

In the following, we consider the terms in the balance equation (2.12) for the nonequilibrium contributions. The order of magnitude of the production term (2.13) and the diffusion-gradient modelling with the flux (2.17) leads us to deduce that $D(k, z) \gg P(k, z)$ at $k \gg L^{-1}$. Therefore, the first-order perturbation to the equilibrium scaling in the inertial range is due to the inhomogeneous diffusion $D(k, z)$. Then, in the inertial range, we have

$$ T(k, z) = -D(k, z), \quad (2.22) $$

and

$$ -\frac{\partial}{\partial k} \Pi_1(k, z) = \frac{\partial}{\partial z} \Phi_0(k, z), \quad (2.23) $$

since $\partial\Pi_0(k, z)/\partial k = 0$. Thus, the first-order correction of the nonlinear transfer balances the zeroth-order contribution of the inhomogeneous diffusion. The first-order perturbation to the nonlinear flux $\Pi_1(k, z)$ is evaluated as (Rubinstein & Clark 2005)

$$ \Pi_1(k, z) = E_1(k, z) \left. \frac{\delta\Pi}{\delta E} \right|_{E_0}, \quad (2.24) $$

where $\delta\Pi/\delta E|_{E_0}$ is the Fréchet derivative of the total flux $\Pi$ evaluated at $E(k, z) = E_0(k, z)$. In the inertial range, assuming $E_1$ to scale as a power law, this yields the scaling,

$$ \Pi_1(k, z) \sim \epsilon(z) \frac{E_1(k, z)}{E_0(k, z)}. \quad (2.25) $$

Note that we obtain (2.25) not only for the Leith model, but also for most of the other classical closures such as the Kovaznay and Heisenberg model (Rubinstein & Clark 2022). Integrating (2.23) from $k$ to $\infty$, we have

$$ \Pi_1(k, z) = \frac{\partial}{\partial z} \int_k^\infty \Phi_0(k, z) dk. \quad (2.26) $$
239 By combining this with (2.17) and (2.25), we obtain

$$240 \qquad E_1(k, z) \sim -\frac{E_0(k, z)}{\epsilon(z)} \frac{\partial}{\partial z} \left( \mu(z) \frac{\partial \int_k^\infty E_0(p, z) dp}{\partial z} \right). \qquad (2.27)$$

241 Substituting (2.21), the above expression gives

$$242 \qquad E_1(k, z) \sim -\mu(z)\epsilon(z)^{1/3}k^{-7/3} \left[ \frac{2}{3}\frac{\epsilon_{zz}(z)}{\epsilon(z)} + \frac{2}{3}\frac{\mu_z(z)}{\mu(z)}\frac{\epsilon_z(z)}{\epsilon(z)} - \frac{2}{9}\left(\frac{\epsilon_z(z)}{\epsilon(z)}\right)^3 \right], \qquad (2.28)$$

243 where the subscripts denote derivatives with respect to z, for example, $\epsilon_{zz} = \partial^2\epsilon(z)/\partial z^2$.

244 We will model the unknown eddy diffusivity in its simplest way,

$$245 \qquad \mu(z) \sim L^{4/3} \epsilon(z)^{1/3}. \qquad (2.29)$$

246 Doing so, we obtain

$$247 \qquad E_1(k, z) \sim -\frac{\epsilon_{zz}(z)L^{4/3}}{\epsilon(z)^{1/3}} k^{-7/3}. \qquad (2.30)$$

248 Note that although all the terms involving $\epsilon_z$ and $\mu_z$ vanish exactly for the current definition of $\mu(z)$ in (2.29), this might not be the case for arbitrary choices of $\mu(z)$. It is therefore possible that, in more sophisticated descriptions of turbulent diffusion, terms involving first-order derivatives might appear. In our description this is however not the case at leading order.

252 Considering expressions (2.25) and (2.30), it is observed that the non-equilibrium contribution $\Pi_1(k, z)$ to the energy flux can change sign, according to the sign of $\epsilon_{zz}(z)$. This implies that this part of the flux can be from large to small $k$ values in the inertial range. However, since $\Pi_1$ is small compared to $\Pi_0$, the total flux in the inertial range, which is the sum of the two contributions, will be towards large $k$ as in the classical three-dimensional picture of an inertial range energy cascade.

## 2.5. Case of a sinusoidal dissipation profile

259 The comparison of expressions (2.21) and (2.30) indicates that the inhomogeneous contribution ($\propto k^{-7/3}$) is subdominant compared to the equilibrium energy spectrum ($\propto k^{-5/3}$) at large wavenumbers. Furthermore, the expression is proportional to the second spatial derivative of the dissipation rate $\epsilon_{zz}(z)$ and can thus be both positive and negative. Let us illustrate the implication of this expression by considering a large-scale inhomogeneity characterised by a cosine function with a characteristic wavelength of order $L$,

$$265 \qquad \epsilon(z) = \langle\epsilon\rangle + \tilde{\epsilon}\cos(z/L), \qquad (2.31)$$

266 with $\langle\epsilon\rangle \gg \tilde{\epsilon}$. We consider $L$, first introduced in § 2.4, to be of the order of and proportional to the characteristic large-scale length of the flow. Substituting this expression for $\epsilon(z)$ in (2.30), we find

$$269 \qquad E_1(k, z) = E_1^X(k) \cos(z/L) \qquad (2.32)$$

270 with

$$271 \qquad E_1^X(k) = C_A \tilde{\epsilon} \langle\epsilon\rangle^{-1/3} L^{-2/3} k^{-7/3}, \qquad (2.33)$$

272 where the superscript $X$ indicates the perturbations due to inhomogeneity.

273 Let us now assume that both the equilibrium spectrum $E_0(k, z)$ and $E_1(k, z)$ extend from $k = L^{-1}$ to $\infty$. Integrating the spectra in this range, we find that

$$275 \qquad K_0(z) \sim L^{2/3} \epsilon(z)^{2/3} \qquad (2.34)$$
239 By combining this with (2.17) and (2.25), we obtain

$$
240
E_1(k, z) \sim -\frac{E_0(k, z)}{\epsilon(z)} \frac{\partial}{\partial z} \left( \mu(z) \frac{\partial \int_k^\infty E_0(p, z) dp}{\partial z} \right). \quad (2.27)
$$

241 Substituting (2.21), the above expression gives

$$
242
E_1(k, z) \sim -\mu(z)\epsilon(z)^{1/3}k^{-7/3} \left[ \frac{2}{3} \frac{\epsilon_{zz}(z)}{\epsilon(z)} + \frac{2}{3} \frac{\mu_z(z)}{\mu(z)} \frac{\epsilon_z(z)}{\epsilon(z)} - \frac{2}{9} \left( \frac{\epsilon_z(z)}{\epsilon(z)} \right) \right], \quad (2.28)
$$

243 where the subscripts denote derivatives with respect to z, for example, $\epsilon_{zz} = \partial^2\epsilon(z)/\partial z^2$.

244 We will model the unknown eddy diffusivity in its simplest way,

$$
245
\mu(z) \sim L^{4/3} \epsilon(z)^{1/3}. \tag{2.29}
$$

246 Doing so, we obtain

$$
247
E_1(k, z) \sim - \frac{\epsilon_{zz}(z) L^{4/3}}{\epsilon(z)^{1/3}} k^{-7/3}. \quad (2.30)
$$

248 Note that although all the terms involving $\epsilon_z$ and $\mu_z$ vanish exactly for the current definition of
249 $\mu(z)$ in (2.29), this might not be the case for arbitrary choices of $\mu(z)$. It is therefore possible
250 that, in more sophisticated descriptions of turbulent diffusion, terms involving first-order
derivatives might appear. In our description this is however not the case at leading order.
251

252 Considering expressions (2.25) and (2.30), it is observed that the non-equilibrium contri-
bution Π₁(k, z) to the energy flux can change sign, according to the sign of ϵ<sub>zz</sub>(z). This
253 implies that this part of the flux can be from large to small k values in the inertial range.
254 However, since Π₁ is small compared to Π₀, the total flux in the inertial range, which is the
255 sum of the two contributions, will be towards large k as in the classical three-dimensional
256 picture of an inertial range energy cascade.
257

258 2.5. Case of a sinusoidal dissipation profile

259 The comparison of expressions (2.21) and (2.30) indicates that the inhomogeneous contri-
bution (∝ k⁻⁷/³) is subdominant compared to the equilibrium energy spectrum (∝ k⁻⁵/³)
260 at large wavenumbers. Furthermore, the expression is proportional to the second spatial
261 derivative of the dissipation rate ϵ<sub>zz</sub>(z) and can thus be both positive and negative. Let
262 us illustrate the implication of this expression by considering a large-scale inhomogeneity
263 characterised by a cosine function with a characteristic wavelength of order L,
264

$$
265
\epsilon(z) = \langle \epsilon \rangle + \tilde{\epsilon} \cos(z/L), \tag{2.31}
$$

266 with $\langle\epsilon\rangle \gg \tilde{\epsilon}$. We consider $L$, first introduced in § 2.4, to be of the order of and proportional to
267 the characteristic large-scale length of the flow. Substituting this expression for $\epsilon(z)$ in (2.30),
268 we find

$$
269
E_1(k, z) = E_1^X(k) \cos(z/L) \tag{2.32}
$$

270 with

$$
271
E_1^X(k) = C_A \tilde{\epsilon} \langle\epsilon\rangle^{-1/3} L^{-2/3} k^{-7/3}, \quad (2.33)
$$

272 where the superscript X indicates the perturbations due to inhomogeneity.

273 Let us now assume that both the equilibrium spectrum $E_0(k, z)$ and $E_1(k, z)$ extend from
274 $k = L^{-1}$ to $\infty$. Integrating the spectra in this range, we find that

$$
275
K_0(z) \sim L^{2/3} \epsilon(z)^{2/3} \qquad (2.34)
$$
<table><caption>Table 1: DNS parameters and statistical quantities. The resolution N and kinematic viscosity ν are the control parameters. The remaining statistical quantities are: the fluctuating isotropic RMS velocity $u' = \sqrt{2K'/3}$ where energy of the temporal fluctuating velocity $K' = \langle u'_i u'_i \rangle_{x,t} / 2$ and $u'_i(x,t) = u_i(x,t) - \langle u_i \rangle_t(x)$; the Taylor microscale $\lambda = 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 $Re_\lambda = u'\lambda/\nu$; the integral time scale $T = L/u'$ with $L = k_f^{-1} = 1$; the simulation time in the statistically steady state $T_{\text{total}}$ as a function of T.</caption><thead><tr><th>N</th><th>ν</th><th>u'</th><th>λ</th><th>Re<sub>λ</sub></th><th>T<sub>total</sub>/T</th></tr></thead><tbody><tr><td>128</td><td>0.07</td><td>1.31</td><td>0.371</td><td>69.6</td><td>959</td></tr><tr><td>256</td><td>0.028</td><td>1.35</td><td>0.233</td><td>113</td><td>645</td></tr><tr><td>512</td><td>0.01</td><td>1.33</td><td>0.138</td><td>184</td><td>170</td></tr></tbody></table>

and

$$K_1(z) \sim -\frac{\epsilon_{zz}(z)L^2}{\langle \epsilon \rangle} K_0(z). \qquad (2.35)$$

Comparing these last two expressions illustrates that the formal expansion parameter in our system is

$$\gamma = \frac{\epsilon_{zz}(z)L^2}{\langle \epsilon \rangle}. \qquad (2.36)$$

The main analytical results of the present investigation [(2.33)–(2.35)] are obtained by phenomenological modelling based on gradient-diffusion assumptions of nonlinear transfer in both physical and scale space. The models and their consequences are, at best, crude approximations of the intricate nonlinear interactions in the actual flow. Therefore, the resulting expressions need verification by experiments or direct numerical simulations.

## 3. Assessment of the inhomogeneous scaling

### 3.1. Numerical set-up

In order to verify the theoretical predictions, in particular expression (2.33), we carry out DNS of three-dimensional Kolmogorov flow in a triple-periodic box. Such flow has the convenient properties of being statistically inhomogeneous in one direction and free of solid boundaries. Furthermore, its properties have been widely investigated numerically (Borue & Orszag 1996; Musacchio & Boffetta 2014; Wu et al. 2021).

The dynamics of the Kolmogorov flow in the present investigation are governed by (2.1) with $f(z) = \sin(k_f z)$. The numerical domain is a cube of size $2\pi$. These choices imply that the forcing wavelength is equal to the width of the cubic domain, and we set $k_f = L^{-1} = 1$. Simulations are carried out using a standard pseudo-spectral solver (Delache et al. 2014) with a third-order Adams-Bashfort time-integration scheme. The details of the simulations are reported in table 1. Since we focus on the effect of inhomogeneity, we attempt to obtain statistics in a steady state over a long-enough time interval to allow the effects of the temporal variations to become as small as possible (see the last column in table 1).

### 3.2. Visualisation and dissipation profile

In the following, we will discuss the simulation at the highest considered Reynolds number $Re_\lambda = 184$. A flow visualisation is shown in figure 2 (a) with the x-component of the velocity field $U_x(x,t)$. The influence of the large-scale mean flow, proportional to the sinusoidal
<table><caption>Table 1: DNS parameters and statistical quantities. The resolution *N* and kinematic viscosity *ν* are the control parameters. The remaining statistical quantities are: the fluctuating isotropic RMS velocity $u' = \sqrt{2K'/3}$ where energy of the temporal fluctuating velocity $K' = \langle u'_i u'_i \rangle_{x,t} / 2$ and $u'_i(x,t) = u_i(x,t) - \langle u_i \rangle_t(x)$; the Taylor microscale $\lambda = 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 $Re_\lambda = u'\lambda/\nu$; the integral time scale $T = L/u'$ with $L = k_f^{-1} = 1$; the simulation time in the statistically steady state $T_{\text{total}}$ as a function of $T$.</caption><thead><tr><th>N</th><th>ν</th><th>u'</th><th>λ</th><th>Re<sub>λ</sub></th><th>T<sub>total</sub>/T</th></tr></thead><tbody><tr><td>128</td><td>0.07</td><td>1.31</td><td>0.371</td><td>69.6</td><td>959</td></tr><tr><td>256</td><td>0.028</td><td>1.35</td><td>0.233</td><td>113</td><td>645</td></tr><tr><td>512</td><td>0.01</td><td>1.33</td><td>0.138</td><td>184</td><td>170</td></tr></tbody></table>

276 and

277

$$K_1(z) \sim -\frac{\epsilon_{zz}(z)L^2}{\langle \epsilon \rangle} K_0(z). \qquad (2.35)$$

278 Comparing these last two expressions illustrates that the formal expansion parameter in our system is

279

280

$$\gamma = \frac{\epsilon_{zz}(z)L^2}{\langle \epsilon \rangle}. \qquad (2.36)$$

281 The main analytical results of the present investigation [(2.33)–(2.35)] are obtained by phenomenological modelling based on gradient-diffusion assumptions of nonlinear transfer in both physical and scale space. The models and their consequences are, at best, crude approximations of the intricate nonlinear interactions in the actual flow. Therefore, the resulting expressions need verification by experiments or direct numerical simulations.

285

## 3. Assessment of the inhomogeneous scaling

287

### 3.1. Numerical set-up

288 In order to verify the theoretical predictions, in particular expression (2.33), we carry out DNS of three-dimensional Kolmogorov flow in a triple-periodic box. Such flow has the convenient properties of being statistically inhomogeneous in one direction and free of solid boundaries. Furthermore, its properties have been widely investigated numerically (Borue & Orszag 1996; Musacchio & Boffetta 2014; Wu et al. 2021).

289 The dynamics of the Kolmogorov flow in the present investigation are governed by (2.1) with $f(z) = \sin(k_f z)$. The numerical domain is a cube of size $2\pi$. These choices imply that the forcing wavelength is equal to the width of the cubic domain, and we set $k_f = L^{-1} = 1$. 295 Simulations are carried out using a standard pseudo-spectral solver (Delache et al. 2014) with a third-order Adams-Bashfort time-integration scheme. The details of the simulations are reported in table 1. Since we focus on the effect of inhomogeneity, we attempt to obtain statistics in a steady state over a long-enough time interval to allow the effects of the temporal variations to become as small as possible (see the last column in table 1).

300

301

### 3.2. Visualisation and dissipation profile

302 In the following, we will discuss the simulation at the highest considered Reynolds number $Re_\lambda = 184$. A flow visualisation is shown in figure 2 (a) with the *x*-component of the velocity field $U_x(x,t)$. The influence of the large-scale mean flow, proportional to the sinusoidal

304
<figure><img src="image_2.png" /><figcaption>Figure 2: (a) Instantaneous distributions of $U_x(x, t)$ at $Re_λ = 184$. Blue (red) corresponds to the negative (positive) value of $U_x$. (b) Instantaneous profiles of $U_x(z, t) = \langle U_x(x, t) \rangle_⊥$ in grey. Time-averaged profile $U(z) = \langle U_x(x, t) \rangle_{⊥,t}$ is indicated by a thick line.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: (a) Instantaneous profile of $\epsilon(z, t) = \langle\epsilon(x, t)\rangle_⊥$. Time-averaged profile $\epsilon(z) = \langle\epsilon(x, t)\rangle_{⊥,t}$ is also shown. The red dashed line denotes $\bar{\epsilon}(z)$, a sinusoidal fitting of $\epsilon(z)$ by (2.31). (b) Time-averaged profile of kinetic energy with fluctuating velocity $K(z)$ and its equilibrium $K_0(z)$ and nonequilibrium $K_1(z)$ contributions. See the main text and Appendix B for the definition.</figcaption></figure>

forcing along the $z$ axis, is distinguishable. Figure 2 (b) shows the instantaneous profile of $U_x(z, t) = \langle U_x(x, t) \rangle_⊥$. The single curve corresponds to the horizontal average of a snapshot, as shown in figure 2 (a). Its time average, $U(z) = \langle U_x(x, t) \rangle_{⊥,t}$, is also shown in figure 2 (b) with a smooth sinusoidal profile.

In figure 3 (a), the instantaneous profile of the energy dissipation rate $\epsilon(z, t) = \langle\epsilon(x, t)\rangle_⊥$ is shown along with its time average $\epsilon(z) = \langle\epsilon(x, t)\rangle_{⊥,t}$. The instantaneous profile shows large fluctuations in comparison to the velocity profile (figure 2 (b)). Its time average, in contrast, shows a smooth sinusoidal profile. This property allows us to use the approximations in § 2.5. As expected, the dissipation peaks at values where the mean velocity gradient is strongest (at $z = 0$ and $\pm\pi$). For numerical convenience, we perform a sinusoidal fitting $\bar{\epsilon}(z)$ introduced in (2.31). This profile is also shown in figure 3 (a).

Figure 3 (b) shows the kinetic energy profile of the fluctuating velocity field. The fluctuating energy profile is defined by $K(z, t) = \mathcal{K}(z, t) - K_U(z)$, where the total energy is $\mathcal{K}(z, t) = U(z)U(z)/2$ and the mean flow energy is $K_U(z) = U(z)^2/2$. We consider
<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Instantaneous distributions of $U_x(x, t)$ at $Re_λ = 184$. Blue (red) corresponds to the negative (positive) value of $U_x$. (b) Instantaneous profiles of $U_x(z, t) = \langle U_x(x, t) \rangle_⊥$ in grey. Time-averaged profile $U(z) = \langle U_x(x, t) \rangle_{⊥,t}$ is indicated by a thick line.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Instantaneous profile of $\epsilon(z, t) = \langle \epsilon(x, t) \rangle_⊥$. Time-averaged profile $\epsilon(z) = \langle \epsilon(x, t) \rangle_{⊥,t}$ is also shown. The red dashed line denotes $\bar{\epsilon}(z)$, a sinusoidal fitting of $\epsilon(z)$ by (2.31). (b) Time-averaged profile of kinetic energy with fluctuating velocity $K(z)$ and its equilibrium $K_0(z)$ and nonequilibrium $K_1(z)$ contributions. See the main text and Appendix B for the definition.</figcaption></figure>

forcing along the $z$ axis, is distinguishable. Figure 2 (b) shows the instantaneous profile of $U_x(z, t) = \langle U_x(x, t) \rangle_⊥$. The single curve corresponds to the horizontal average of a snapshot, as shown in figure 2 (a). Its time average, $U(z) = \langle U_x(x, t) \rangle_{⊥,t}$, is also shown in figure 2 (b) with a smooth sinusoidal profile.

In figure 3 (a), the instantaneous profile of the energy dissipation rate $\epsilon(z, t) = \langle \epsilon(x, t) \rangle_⊥$ is shown along with its time average $\epsilon(z) = \langle \epsilon(x, t) \rangle_{⊥,t}$. The instantaneous profile shows large fluctuations in comparison to the velocity profile (figure 2 (b)). Its time average, in contrast, shows a smooth sinusoidal profile. This property allows us to use the approximations in § 2.5. As expected, the dissipation peaks at values where the mean velocity gradient is strongest (at $z = 0$ and $\pm\pi$). For numerical convenience, we perform a sinusoidal fitting $\bar{\epsilon}(z)$ introduced in (2.31). This profile is also shown in figure 3 (a).

Figure 3 (b) shows the kinetic energy profile of the fluctuating velocity field. The fluctuating energy profile is defined by $K(z, t) = \mathcal{K}(z, t) - K_U(z)$, where the total energy is $\mathcal{K}(z, t) = U(z)U(z)/2$ and the mean flow energy is $K_U(z) = U(z)^2/2$. We consider
<figure><img src="image_4.png" /><figcaption>Figure 4: Time-averaged three-dimensional isotropic energy spectrum $E(k) = \langle E(k, t) \rangle_t$, normalised by Kolmogorov variables. Results are shown at $Re_\lambda = 69.6, 113$ and $184$ (see table 1). The red dashed line denotes the $k^{-5/3}$ scaling for reference.</figcaption></figure>

319 the decomposition, see (2.34)–(2.35) and Appendix B,

$$K(z, t) = K_0(z, t) + K_1(z, t). \qquad (3.1)$$

320 In figure 3 (b), we observe that the equilibrium $K_0(z) = \langle K_0(z, t) \rangle_t$ and the nonequilibrium $K_1(z) = \langle K_1(z, t) \rangle_t$ profiles share the same phase, consistent with the prediction that the spectrum $E_1(k, z)$ is proportional to $-\partial^2\epsilon(z)/\partial z^2$.

324 3.3. Equilibrium and nonequilibrium spectra

325 Figure 4 shows the isotropic energy spectrum $E(k, t)$ (see Appendix A for the definition) at 326 three different Taylor-length Reynolds numbers. For simplicity, we denote its time-average by 327 $E(k) = \langle E(k, t) \rangle_t$. Normalisation using $\nu$ and $\epsilon = \langle \epsilon(x, t) \rangle_{x,t}$ allows an excellent collapse 328 for large values of $k$.

329 Next, we assess energy spectra in statistically homogeneous planes perpendicular to the $z$ axis, as defined in (2.9). In the following, we analyse the time-averaged inhomogeneous energy spectrum $E(k_\perp, z) = \langle E(k_\perp, z, t) \rangle_t$ in a statistically steady state (see table 1). 330 Figure 5 (a) shows $E(k_\perp, z)$ nondimensionalised by $\bar{\epsilon}(z)^{2/3}L^{2/3}$. The fluctuations at small 331 scales are small, and variations are barely visible.

332 Then, we assume that the energy spectrum in the inertial range can be written, in a general 333 form,

$$E_0(k_{\perp}, z) \sim \epsilon(z)^{2/3} k_{\perp}^{-5/3} f_L[k_{\perp}L] f_{\eta}[k_{\perp}\eta(z)], \qquad (3.2)$$

334 where $f_L$ and $f_\eta$ determine the shape at small and large $k$, respectively. Therefore, there are 335 two distinct choices to collapse the spectra. See Appendix B for the detail of two normalisation 336 procedures. In this study, we employ the large-scale normalisation and evaluate

$$f_L[k_{\perp}L] = \langle E(k_{\perp}, z) \bar{\epsilon}(z)^{-2/3} k_{\perp}^{5/3} \rangle_z, \qquad (3.3)$$

340 as shown in figure 5 (a). Note that this expression is valid for $k_\perp\eta \ll 1$ where $f_\eta(k_\perp\eta)$ tends 341 to a constant value. Then, the equilibrium spectrum can be defined as

$$E_0(k_{\perp}, z) \equiv \bar{\epsilon}(z)^{2/3} f_L[k_{\perp}L] k_{\perp}^{-5/3}. \qquad (3.4)$$

344 Now, we can evaluate the nonequilibrium spectrum by $E_1(k_\perp, z) = E(k_\perp, z) - E_0(k_\perp, z)$.
<figure><img src="image_5.png" /><figcaption>Figure 4: Time-averaged three-dimensional isotropic energy spectrum $E(k) = \langle E(k, t) \rangle_t$, normalised by Kolmogorov variables. Results are shown at $Re_\lambda = 69.6, 113$ and $184$ (see table 1). The red dashed line denotes the $k^{-5/3}$ scaling for reference.</figcaption></figure>

319 the decomposition, see (2.34)–(2.35) and Appendix B,

$$K(z, t) = K_0(z, t) + K_1(z, t). \qquad (3.1)$$

320 In figure 3 (b), we observe that the equilibrium $K_0(z) = \langle K_0(z, t) \rangle_t$ and the nonequilibrium $K_1(z) = \langle K_1(z, t) \rangle_t$ profiles share the same phase, consistent with the prediction that the spectrum $E_1(k, z)$ is proportional to $-\partial^2\epsilon(z)/\partial z^2$.

324 3.3. Equilibrium and nonequilibrium spectra

325 Figure 4 shows the isotropic energy spectrum $E(k, t)$ (see Appendix A for the definition) at 326 three different Taylor-length Reynolds numbers. For simplicity, we denote its time-average by 327 $E(k) = \langle E(k, t) \rangle_t$. Normalisation using $\nu$ and $\epsilon = \langle \epsilon(x, t) \rangle_{x,t}$ allows an excellent collapse 328 for large values of $k$.

329 Next, we assess energy spectra in statistically homogeneous planes perpendicular to the $z$ axis, as defined in (2.9). In the following, we analyse the time-averaged inhomogeneous energy spectrum $E(k_\perp, z) = \langle E(k_\perp, z, t) \rangle_t$ in a statistically steady state (see table 1). 330 Figure 5 (a) shows $E(k_\perp, z)$ nondimensionalised by $\bar{\epsilon}(z)^{2/3}L^{2/3}$. The fluctuations at small 331 scales are small, and variations are barely visible.

332 Then, we assume that the energy spectrum in the inertial range can be written, in a general 333 form,

$$E_0(k_{\perp}, z) \sim \epsilon(z)^{2/3} k_{\perp}^{-5/3} f_L[k_{\perp}L] f_{\eta}[k_{\perp}\eta(z)], \qquad (3.2)$$

334 where $f_L$ and $f_\eta$ determine the shape at small and large $k$, respectively. Therefore, there are 335 two distinct choices to collapse the spectra. See Appendix B for the detail of two normalisation 336 procedures. In this study, we employ the large-scale normalisation and evaluate

$$f_L[k_{\perp}L] \equiv \langle E(k_{\perp}, z) \bar{\epsilon}(z)^{-2/3} k_{\perp}^{5/3} \rangle_z, \qquad (3.3)$$

340 as shown in figure 5 (a). Note that this expression is valid for $k_\perp\eta \ll 1$ where $f_\eta(k_\perp\eta)$ tends 341 to a constant value. Then, the equilibrium spectrum can be defined as

$$E_0(k_{\perp}, z) \equiv \bar{\epsilon}(z)^{2/3} f_L[k_{\perp}L] k_{\perp}^{-5/3}. \qquad (3.4)$$

344 Now, we can evaluate the nonequilibrium spectrum by $E_1(k_\perp, z) \equiv E(k_\perp, z) - E_0(k_\perp, z)$.
<figure><img src="image_5.png" /><figcaption>Figure 5: (a) Nondimensionalised two-dimensional energy spectrum. Note that $E(k_{\perp}, z) = \langle E(k_{\perp}, z, t) \rangle_t$. Dark (light) colour represents the small (large) value of z coordinate. The thick black line denotes (3.3), the average over z coordinate.</figcaption></figure>

(b) Time-averaged nonequilibrium energy spectrum with specific signs: $E_1^+(k_{\perp}) = \langle E_1(k_{\perp}, z) > 0 \rangle_z$ and $E_1^-(k_{\perp}) = \langle E_1(k_{\perp}, z) < 0 \rangle_z$. Red dashed lines denote the $k_{\perp}^{-7/3}$ slope.

<figure><img src="image_6.png" /><figcaption>Figure 6: (a) Absolute value of the time-averaged nonequilibrium energy spectrum $|E_1(k_{\perp}, z)| = |\langle E_1(k_{\perp}, z, t) \rangle_t|$ for three values of the Taylor-length Reynolds numbers. The red dashed line represents the $k_{\perp}^{-7/3}$ scaling. (b) Compensated spectrum of the panel (a). The red dashed line denotes the compensated $k_{\perp}^{-7/3}$ scaling.</figcaption></figure>

Note that i) this quantity is defined by the time-averaged spectra and ii) since this quantity can be regarded as a perturbation of $E(k_{\perp}, z)$ around $E_0(k_{\perp}, z)$, it can be both positive and negative. Figure 5 (b) shows the z-average of $E(k_{\perp}, z)$ for specific signs.

Similar plots are shown in figure 10 of Horiuti & Ozawa (2011) and figure 2 of Horiuti & Tamaki (2013). The scaling is consistent with the one derived in § 2.4,

$$ \langle |E_1(k_{\perp}, z)| \rangle_z \propto k_{\perp}^{-7/3}. \qquad (3.5) $$

Figure 6 (a) compares $\langle |E_1(k_{\perp}, z)| \rangle_z$ for three different Reynolds numbers as in figure 4. For smaller values of $Re_{\lambda}$, the spectrum exhibits steeper scaling than $k_{\perp}^{-7/3}$. At larger $Re_{\lambda}$, the
<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Nondimensionalised two-dimensional energy spectrum. Note that $E(k_{\perp}, z) = \langle E(k_{\perp}, z, t) \rangle_t$. Dark (light) colour represents the small (large) value of z coordinate. The thick black line denotes (3.3), the average over z coordinate.</figcaption></figure>

(b) Time-averaged nonequilibrium energy spectrum with specific signs: $E_1^+(k_{\perp}) = \langle E_1(k_{\perp}, z) > 0 \rangle_z$ and $E_1^-(k_{\perp}) = \langle E_1(k_{\perp}, z) < 0 \rangle_z$. Red dashed lines denote the $k_{\perp}^{-7/3}$ slope.

<figure><img src="image_7.png" /><figcaption>Figure 6: (a) Absolute value of the time-averaged nonequilibrium energy spectrum $|E_1(k_{\perp}, z)| = |\langle E_1(k_{\perp}, z, t) \rangle_t|$ for three values of the Taylor-length Reynolds numbers. The red dashed line represents the $k_{\perp}^{-7/3}$ scaling. (b) Compensated spectrum of the panel (a). The red dashed line denotes the compensated $k_{\perp}^{-7/3}$ scaling.</figcaption></figure>

Note that i) this quantity is defined by the time-averaged spectra and ii) since this quantity can be regarded as a perturbation of $E(k_{\perp}, z)$ around $E_0(k_{\perp}, z)$, it can be both positive and negative. Figure 5 (b) shows the z-average of $E(k_{\perp}, z)$ for specific signs.

Similar plots are shown in figure 10 of Horiuti & Ozawa (2011) and figure 2 of Horiuti & Tamaki (2013). The scaling is consistent with the one derived in § 2.4,

$$ \langle |E_1(k_{\perp}, z)| \rangle_z \propto k_{\perp}^{-7/3}. \qquad (3.5) $$

Figure 6 (a) compares $\langle |E_1(k_{\perp}, z)| \rangle_z$ for three different Reynolds numbers as in figure 4. For smaller values of $Re_{\lambda}$, the spectrum exhibits steeper scaling than $k_{\perp}^{-7/3}$. At larger $Re_{\lambda}$, the
353 slope approaches the $k_{\perp}^{-7/3}$ scaling. At the same time, the spectrum in the higher $k_{\perp}$ range
354 exhibits a bump associated with shallower scaling than $k_{\perp}^{-7/3}$.

355 We plot the compensated spectra in figure 6 (b). Although the scaling range extends for
356 less than a decade, the emergence of the $k_{\perp}^{-7/3}$ scaling range is well captured using this
357 normalisation. In Appendix B, we confirm that the bump in the compensated spectra is due
358 to our choice of the nondimensional function (3.2).

## 4. Conclusion and discussion

360 The numerical simulations in the previous section support our prediction,

$$E(k, z) = C_K \epsilon(z)^{2/3} k^{-5/3} - C_A \frac{\epsilon_{zz}(z) L^{4/3}}{\epsilon(z)^{1/3}} k^{-7/3}, \quad (4.1)$$

362 of the energy spectrum for turbulence with inhomogeneity in the z direction. This scaling
363 quantifies the influence of spatial inhomogeneity in wavenumber space. In particular, the
364 special case where the dissipation fluctuates as a sinusoidal function around a mean value,
365 discussed in §2.5, gives us a useful estimate of the influence of inhomogeneity (2.32)–(2.33).
366 Indeed, introducing an average dissipation $\langle\epsilon\rangle$ and smooth spatial fluctuations $\tilde{\epsilon}$ around $\langle\epsilon\rangle$,
367 so that $\epsilon_{zz}/\tilde{\epsilon} \sim L^{-2}$, we obtain that

$$\frac{E_1^X(k, z)}{E_0(k, z)} \sim \frac{\tilde{\epsilon}}{\langle\epsilon\rangle} (kL)^{-2/3}. \quad (4.2)$$

369 This expression shows that the influence of large-scale inhomogeneity is negligible for

$$k \gg L^{-1} \left( \frac{\tilde{\epsilon}}{\langle \epsilon \rangle} \right)^{3/2}. \quad (4.3)$$

371 Therefore, if this requirement is fulfilled in a statistically stationary flow, far enough away
372 from walls, Kolmogorov's equilibrium spectrum is expected to be dominant compared to the
373 contributions associated with spatial inhomogeneity.

## Appendix A. Governing equation of the inhomogeneous energy spectrum

374 In this Appendix, we define the spectrum and the governing equations for the energy spectrum
375 tensor in inhomogeneous flow (see (2.12)). The generalised spectrum $E(k, x) = E_{ii}(k, x)$
376 is defined by

$$E(k, x) = \frac{1}{2} \iint e^{-ik \cdot r} \left\langle u_i \left( x + \frac{r}{2} \right) u_i \left( x - \frac{r}{2} \right) \right\rangle dr \, d\Omega_k, \quad (A 1)$$

378 where $\int d\Omega_k$ denotes the integral over spherical shells of radius $k$. The brackets in this
379 Appendix section denote the ensemble average. For statistically homogeneous turbulence,
380 this definition is equivalent to the expression

$$E(k) = \int \frac{1}{2} \langle u_i(k) u_i^*(k) \rangle d\Omega_k. \quad (A 2)$$

383 The evolution equation for $E_{ij}(k, x, t)$ formally reads

$$\frac{\partial E_{ij}(k, x, t)}{\partial t} = \frac{1}{2} \int \left[ \int e^{-ik \cdot r} \Psi_{ij}(x+r/2, x-r/2, t) dr \right] d\Omega_k. \quad (A 3)$$
353 slope approaches the $k_{\perp}^{-7/3}$ scaling. At the same time, the spectrum in the higher $k_{\perp}$ range
354 exhibits a bump associated with shallower scaling than $k_{\perp}^{-7/3}$.

355 We plot the compensated spectra in figure 6 (b). Although the scaling range extends for
356 less than a decade, the emergence of the $k_{\perp}^{-7/3}$ scaling range is well captured using this
357 normalisation. In Appendix B, we confirm that the bump in the compensated spectra is due
358 to our choice of the nondimensional function (3.2).

## 4. Conclusion and discussion

360 The numerical simulations in the previous section support our prediction,

$$E(k, z) = C_K \epsilon(z)^{2/3} k^{-5/3} - C_A \frac{\epsilon_{zz}(z) L^{4/3}}{\epsilon(z)^{1/3}} k^{-7/3}, \quad (4.1)$$

362 of the energy spectrum for turbulence with inhomogeneity in the z direction. This scaling
363 quantifies the influence of spatial inhomogeneity in wavenumber space. In particular, the
364 special case where the dissipation fluctuates as a sinusoidal function around a mean value,
365 discussed in §2.5, gives us a useful estimate of the influence of inhomogeneity (2.32)–(2.33).
366 Indeed, introducing an average dissipation $\langle\epsilon\rangle$ and smooth spatial fluctuations $\tilde{\epsilon}$ around $\langle\epsilon\rangle$,
367 so that $\epsilon_{zz}/\tilde{\epsilon} \sim L^{-2}$, we obtain that

$$\frac{E_1^X(k, z)}{E_0(k, z)} \sim \frac{\tilde{\epsilon}}{\langle\epsilon\rangle} (kL)^{-2/3}. \quad (4.2)$$

369 This expression shows that the influence of large-scale inhomogeneity is negligible for

$$k \gg L^{-1} \left( \frac{\tilde{\epsilon}}{\langle \epsilon \rangle} \right)^{3/2}. \quad (4.3)$$

371 Therefore, if this requirement is fulfilled in a statistically stationary flow, far enough away
372 from walls, Kolmogorov's equilibrium spectrum is expected to be dominant compared to the
373 contributions associated with spatial inhomogeneity.

## Appendix A. Governing equation of the inhomogeneous energy spectrum

374 In this Appendix, we define the spectrum and the governing equations for the energy spectrum
375 tensor in inhomogeneous flow (see (2.12)). The generalised spectrum $E(k, x) = E_{ii}(k, x)$
376 is defined by

$$E(k, x) = \frac{1}{2} \iint e^{-ik \cdot r} \left\langle u_i \left( x + \frac{r}{2} \right) u_i \left( x - \frac{r}{2} \right) \right\rangle dr \, d\Omega_k, \quad (A 1)$$

378 where $\int d\Omega_k$ denotes the integral over spherical shells of radius $k$. The brackets in this
379 Appendix section denote the ensemble average. For statistically homogeneous turbulence,
380 this definition is equivalent to the expression

$$E(k) = \int \frac{1}{2} \langle u_i(k) u_i^*(k) \rangle d\Omega_k. \quad (A 2)$$

383 The evolution equation for $E_{ij}(k, x, t)$ formally reads

$$\frac{\partial E_{ij}(k, x, t)}{\partial t} = \frac{1}{2} \int \left[ \int e^{-ik \cdot r} \Psi_{ij}(x+r/2, x-r/2, t) dr \right] d\Omega_k. \quad (A 3)$$
For the tensor on the RHS, we have

$$ \Psi_{ij}(x_1, x_2) = v(\nabla_1^2 + \nabla_2^2)R_{ij}(x_1, x_2) \\ - \left[ \frac{\partial}{\partial x_{1n}} U_n(x_1)R_{ij}(x_1, x_2) + \frac{\partial}{\partial x_{2n}} U_n(x_2)R_{ij}(x_1, x_2) \right. \\ + \frac{\partial}{\partial x_{1n}} U_i(x_1)R_{nj}(x_1, x_2) + \frac{\partial}{\partial x_{2n}} U_j(x_2)R_{ni}(x_1, x_2) \\ + \frac{\partial}{\partial x_{1i}} \langle \mathcal{P}(x_1)u_j(x_2) \rangle + \frac{\partial}{\partial x_{2j}} \langle \mathcal{P}(x_2)u_i(x_1) \rangle \\ \left. + \frac{\partial}{\partial x_{1n}} \langle u_i(x_1)u_n(x_1)u_j(x_2) \rangle + \frac{\partial}{\partial x_{2n}} \langle u_j(x_2)u_n(x_2)u_i(x_1) \rangle \right]. \quad (A4) $$

In this expression and the following, the argument *t* for time is omitted for visibility. The two-point velocity tensor is defined by

$$ R_{ij}(x_1, x_2) \equiv \langle u_i(x_1)u_j(x_2) \rangle, \quad (A5) $$

and the equation (A4) is completed by incompressibility conditions for the mean field and the fluctuations.

Both inhomogeneous turbulence diffusion and spectral transfer are associated with the last two lines of expression (A4). An assumption of weak inhomogeneity must be invoked to dissociate them to obtain a closed expression. Even retaining only the leading order terms in an expansion about inhomogeneity, the resulting equations become quite cumbersome (see Laporta (1995); Besnard *et al.* (1996)).

Subsequently, the different terms in (A4) need to be modelled to close the triple correlations. We will not proceed in this direction and will directly model them by their physical effects. See (2.13)–(2.16).

## Appendix B. Normalisation to extract the non-equilibrium spectrum and kinetic energy profile

In this Appendix, we investigate the different normalisation procedures mentioned in § 3.3.

We state that the energy spectrum is decomposed into equilibrium (labelled by the subscript 0) and nonequilibrium (labelled by the subscript 1) contributions as

$$ E(k, z) = E_0(k, z) + E_1(k, z). \quad (B1) $$

We further assume that the nonequilibrium contributions are zero-mean,

$$ \langle E_1(k, z) \rangle_z = 0. \quad (B2) $$

We have therefore

$$ \langle E(k, z) \rangle_z = \langle E_0(k, z) \rangle_z. \quad (B3) $$

In order to compute the nonequilibrium contributions $E_1(k, z) = E(k, z) - E_0(k, z)$, we need to know the z-dependence of $E_0(k, z)$. For this purpose, we use self-similarity assumptions and Kolmogorov's equilibrium hypothesis.

Scaling ranges in turbulence spectra appear when scale separation is attained, i.e., at sufficiently high Reynolds numbers. In general, one can write the energy spectrum to scale as (see (3.2)),

$$ E_0(k, z) \sim \epsilon(z)^{2/3} k^{-5/3} f_L[kL] f_\eta[k\eta(z)], \quad (B4) $$

at high-Reynolds numbers. We have two nondimensional functions in (B4); the $f_L$ determines
For the tensor on the RHS, we have

$$ \Psi_{ij}(x_1, x_2) = v(\nabla_1^2 + \nabla_2^2) R_{ij}(x_1, x_2) \\ - \left[ \frac{\partial}{\partial x_{1n}} U_n(x_1) R_{ij}(x_1, x_2) + \frac{\partial}{\partial x_{2n}} U_n(x_2) R_{ij}(x_1, x_2) \right. \\ + \frac{\partial}{\partial x_{1n}} U_i(x_1) R_{nj}(x_1, x_2) + \frac{\partial}{\partial x_{2n}} U_j(x_2) R_{ni}(x_1, x_2) \\ + \frac{\partial}{\partial x_{1i}} \langle \mathcal{P}(x_1) u_j(x_2) \rangle + \frac{\partial}{\partial x_{2j}} \langle \mathcal{P}(x_2) u_i(x_1) \rangle \\ \left. + \frac{\partial}{\partial x_{1n}} \langle u_i(x_1) u_n(x_1) u_j(x_2) \rangle + \frac{\partial}{\partial x_{2n}} \langle u_j(x_2) u_n(x_2) u_i(x_1) \rangle \right]. \quad (A4) $$

In this expression and the following, the argument *t* for time is omitted for visibility. The two-point velocity tensor is defined by

$$ R_{ij}(x_1, x_2) \equiv \langle u_i(x_1) u_j(x_2) \rangle, \quad (A5) $$

and the equation (A4) is completed by incompressibility conditions for the mean field and the fluctuations.

Both inhomogeneous turbulence diffusion and spectral transfer are associated with the last two lines of expression (A4). An assumption of weak inhomogeneity must be invoked to dissociate them to obtain a closed expression. Even retaining only the leading order terms in an expansion about inhomogeneity, the resulting equations become quite cumbersome (see Laporta (1995); Besnard et al. (1996)).

Subsequently, the different terms in (A4) need to be modelled to close the triple correlations. We will not proceed in this direction and will directly model them by their physical effects. See (2.13)–(2.16).

## Appendix B. Normalisation to extract the non-equilibrium spectrum and kinetic energy profile

In this Appendix, we investigate the different normalisation procedures mentioned in § 3.3.

We state that the energy spectrum is decomposed into equilibrium (labelled by the subscript 0) and nonequilibrium (labelled by the subscript 1) contributions as

$$ E(k, z) = E_0(k, z) + E_1(k, z). \quad (B1) $$

We further assume that the nonequilibrium contributions are zero-mean,

$$ \langle E_1(k, z) \rangle_z = 0. \quad (B2) $$

We have therefore

$$ \langle E(k, z) \rangle_z = \langle E_0(k, z) \rangle_z. \quad (B3) $$

In order to compute the nonequilibrium contributions $E_1(k, z) = E(k, z) - E_0(k, z)$, we need to know the z-dependence of $E_0(k, z)$. For this purpose, we use self-similarity assumptions and Kolmogorov's equilibrium hypothesis.

Scaling ranges in turbulence spectra appear when scale separation is attained, i.e., at sufficiently high Reynolds numbers. In general, one can write the energy spectrum to scale as (see (3.2)),

$$ E_0(k, z) \sim \epsilon(z)^{2/3} k^{-5/3} f_L[kL] f_\eta[k\eta(z)], \quad (B4) $$

at high-Reynolds numbers. We have two nondimensional functions in (B4); the $f_L$ determines
<figure><img src="image_7.png" /><figcaption>Figure 7: Absolute value of the time-averaged nonequilibrium energy spectrum $|E_1(k_\perp, z)| = |\langle E_1(k_\perp, z, t) \rangle_t|$ for the highest Reynolds number dataset. Different nondimensionalised functions are employed to compute the nonequilibrium spectrum; panel (a) with (B 7) and (b) with (B 9), respectively. The black solid and red dashed lines denote the z-average and $k_\perp^{-7/3}$ scaling, respectively.</figcaption></figure>

the shape of the spectrum for small $k$ (large-scale) and $f_\eta$ for large $k$ (small-scale), respectively. These functions satisfy the framework of Kolmogorov (1941),

$$ \lim_{x \to 0} f_{\eta}[x] = \lim_{x \to \infty} f_{L}[x] = 1. \qquad (B 5) $$

Therefore, we retrieve (1.1) for scales $L^{-1} \ll k \ll \eta^{-1}$ in the limit of $L/\eta \to \infty$. Multiplying both sides of (B 4) by $\eta^{-5/3}$ and dividing by $\epsilon^{2/3}$, we obtain

$$ \frac{E_0(k, z)}{\epsilon(z)^{1/4} \nu^{5/4}} = F_{\eta}[k\eta(z)] f_L[kL], \qquad (B 6) $$

with

$$ F_{\eta}[k\eta(z)] = (k\eta)^{-5/3} f_{\eta}[k\eta(z)]. \qquad (B 7) $$

Since $F(kL)$ tends to unity for $k \gg 1/L$, the equilibrium spectra $E_0(k, z)$ should collapse when normalised by (B 6) for any $z$, for large $kL$.

Similarly, if the large scales are characterised by a length scale $L$, we can propose an alternative normalisation for (B 6),

$$ \frac{E_0(k, z)}{\epsilon(z)^{2/3} L^{5/3}} \sim F_L[kL] f[k\eta], \qquad (B 8) $$

with

$$ F_L[kL] = (kL)^{-5/3} f_L[kL]. \qquad (B 9) $$

It should scale purely as a function of $kL$ for $k\eta \ll 1$.

There are, therefore, two normalisation possibilities. One focuses on the high wavenumber limit of the inertial range close to the dissipation range (B 6)–(B 7), the other one on the low $k$ range close to the energy-range (B 8)–(B 9). In the limit of infinite Reynolds number, we should find them to be equivalent in the inertial range, since

$$ \lim_{x \to \infty} F_L[x] = \lim_{x \to 0} F_\eta[x] = x^{-5/3}. \qquad (B 10) $$

Figure 7 plots the absolute value of the nonequilibrium energy spectrum $E_1(k_\perp, z)$ for these two normalisations. It is observed that using $F[kL]$ we reveal a larger inertial range. We will therefore use this normalisation in the present investigation.
<figure><img src="image_8.png" /><figcaption>Figure 7: Absolute value of the time-averaged nonequilibrium energy spectrum $|E_1(k_\perp, z)| = |\langle E_1(k_\perp, z, t) \rangle_t|$ for the highest Reynolds number dataset. Different nondimensionalised functions are employed to compute the nonequilibrium spectrum; panel (a) with (B 7) and (b) with (B 9), respectively. The black solid and red dashed lines denote the z-average and $k_\perp^{-7/3}$ scaling, respectively.</figcaption></figure>

416 the shape of the spectrum for small $k$ (large-scale) and $f_\eta$ for large $k$ (small-scale),
417 respectively. These functions satisfy the framework of Kolmogorov (1941),

$$ \lim_{x \to 0} f_{\eta}[x] = \lim_{x \to \infty} f_{L}[x] = 1. \qquad (B 5) $$

419 Therefore, we retrieve (1.1) for scales $L^{-1} \ll k \ll \eta^{-1}$ in the limit of $L/\eta \to \infty$. Multiplying both sides of (B 4) by $\eta^{-5/3}$ and dividing by $\epsilon^{2/3}$, we obtain

$$ \frac{E_0(k, z)}{\epsilon(z)^{1/4} \nu^{5/4}} = F_{\eta}[k\eta(z)] f_L[kL], \qquad (B 6) $$

422 with

$$ F_{\eta}[k\eta(z)] = (k\eta)^{-5/3} f_{\eta}[k\eta(z)]. \qquad (B 7) $$

424 Since $F(kL)$ tends to unity for $k \gg 1/L$, the equilibrium spectra $E_0(k, z)$ should collapse when normalised by (B 6) for any $z$, for large $kL$.

426 Similarly, if the large scales are characterised by a length scale $L$, we can propose an
427 alternative normalisation for (B 6),

$$ \frac{E_0(k, z)}{\epsilon(z)^{2/3} L^{5/3}} \sim F_L[kL] f[k\eta], \qquad (B 8) $$

429 with

$$ F_L[kL] = (kL)^{-5/3} f_L[kL]. \qquad (B 9) $$

431 It should scale purely as a function of $kL$ for $k\eta \ll 1$.

432 There are, therefore, two normalisation possibilities. One focuses on the high wavenumber
433 limit of the inertial range close to the dissipation range (B 6)–(B 7), the other one on the low
434 $k$ range close to the energy-range (B 8)–(B 9). In the limit of infinite Reynolds number, we
435 should find them to be equivalent in the inertial range, since

$$ \lim_{x \to \infty} F_L[x] = \lim_{x \to 0} F_\eta[x] = x^{-5/3}. \qquad (B 10) $$

437 Figure 7 plots the absolute value of the nonequilibrium energy spectrum $E_1(k_\perp, z)$ for these
438 two normalisations. It is observed that using $F[kL]$ we reveal a larger inertial range. We will
439 therefore use this normalisation in the present investigation.
440 We have checked the data and we verified that the spectra $E_1(k_\perp, z)$ are for a given $z$ of
441 a single sign, with a very occasional sign-change when the statistics are determined over a
442 too short time-interval. Perturbation spectra are, as predicted, proportional to $-\epsilon_{zz}(z)$ as is
443 also the case for $K_1(z)$, as can be seen by comparing Figs. 3(a) and (b)). The procedure to
444 determine $K_1(z)$ is now outlined.

445 Since the kinetic energy is dominantly determined by large scales, (B 8) allows us to
446 determine the equilibrium kinetic energy profile,

$$447 \qquad K_0(z) = \int E_0(k, z) dk = C_L \epsilon(z)^{2/3} L^{2/3} \qquad (B 11)$$

448 with $C_L = \int F_L[x] dx$. Then, we define the decomposition

$$449 \qquad K_0(z) = \langle K(z) \rangle_z + \tilde{K}_0(z), \qquad (B 12)$$

450 where $\langle K_0(z) \rangle_z = \langle K(z) \rangle_z$ follows from the assumption that $\langle K_1(z) \rangle_z = 0$. By employing
451 the decomposition for the energy dissipation rate profile

$$452 \qquad \epsilon(z)^{2/3} = \langle \epsilon(z)^{2/3} \rangle_z + \overline{\epsilon(z)^{2/3}}, \qquad (B 13)$$

453 it follows from (B 11) that

$$454 \qquad \frac{\tilde{K}_0(z)}{\langle K_0(z) \rangle_z} = \frac{\overline{\epsilon(z)^{2/3}}}{\langle \epsilon(z)^{2/3} \rangle_z}, \qquad (B 14)$$

455 and by (B 12),

$$456 \qquad K_0(z) = \left( 1 + \frac{\overline{\epsilon(z)^{2/3}}}{\langle \epsilon(z)^{2/3} \rangle_z} \right) \langle K(z) \rangle_z . \qquad (B 15)$$

457 Since all the terms on the RHS are known, one can evaluate the nonequilibrium kinetic
458 energy profile $K_1(z) = K(z) - K_0(z)$ (see Equation (3.1)).

459 Acknowledgements. We dedicate this work to Jean-Pierre Bertoglio, on the occasion of his retirement. All
460 DNS simulations were carried out using the facilities of the PMCS2I (École Centrale de Lyon). R.A. is
461 supported by the Takenaka Scholarship Foundation.

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

465 Declaration of interests. The authors report no conflict of interest.

REFERENCES

466 BERTOGLIO, JEAN-PIERRE, BATAILLE, FRANÇOISE & MARION, JEAN-DENIS 2001 Two-point closures for weakly compressible turbulence. *Physics of Fluids* **13** (1), 290–310.

468 BERTOGLIO, JEAN-PIERRE & JEANDEL, DENIS 1987 A simplified spectral closure for inhomogeneous turbulence: application to the boundary layer. In *Turbulent Shear Flows* 5, pp. 19–30. Springer.

470 BESNARD, DIDIER C., HARLOW, FRANCIS H., RAUENZAHN, RICK M. & ZEMACH, CHARLES 1996 Spectral transport model for turbulence. *Theoretical and Computational Fluid Dynamics* **8** (1), 1–35.

472 BORUE, VADIM & ORSZAG, STEVEN A. 1996 Numerical study of three-dimensional Kolmogorov flow at high Reynolds numbers. *Journal of Fluid Mechanics* **306**, 293–323.

474 BOS, WOUTER J. T. 2020 Production and dissipation of kinetic energy in grid turbulence. *Physical Review Fluids* **5** (10), 104607.

476 BOS, WOUTER J. T. & RUBINSTEIN, ROBERT 2017 Dissipation in unsteady turbulence. *Physical Review Fluids* **2** (2), 022601.

478 BRIARD, ANTOINE, GRÉA, BENOÎT-JOSEPH, MONS, VINCENT, CAMBON, CLAUDE, GOMEZ, THOMAS & SAGAUT,
440 We have checked the data and we verified that the spectra $E_1(k_\perp, z)$ are for a given $z$ of
441 a single sign, with a very occasional sign-change when the statistics are determined over a
442 too short time-interval. Perturbation spectra are, as predicted, proportional to $-\epsilon_{zz}(z)$ as is
443 also the case for $K_1(z)$, as can be seen by comparing Figs. 3(a) and (b)). The procedure to
444 determine $K_1(z)$ is now outlined.

445 Since the kinetic energy is dominantly determined by large scales, (B 8) allows us to
446 determine the equilibrium kinetic energy profile,

$$
447 \qquad K_0(z) = \int E_0(k, z) dk = C_L \epsilon(z)^{2/3} L^{2/3} \qquad (B 11)
$$

448 with $C_L = \int F_L[x] \, dx$. Then, we define the decomposition

$$
449 \qquad K_0(z) = \langle K(z) \rangle_z + \tilde{K}_0(z), \qquad (B 12)
$$

450 where $\langle K_0(z) \rangle_z = \langle K(z) \rangle_z$ follows from the assumption that $\langle K_1(z) \rangle_z = 0$. By employing
451 the decomposition for the energy dissipation rate profile

$$
452 \qquad \epsilon(z)^{2/3} = \langle \epsilon(z)^{2/3} \rangle_z + \overline{\epsilon(z)^{2/3}}, \qquad (B 13)
$$

453 it follows from (B 11) that

$$
454 \qquad \frac{\tilde{K}_0(z)}{\langle K_0(z) \rangle_z} = \frac{\overline{\epsilon(z)^{2/3}}}{\langle \epsilon(z)^{2/3} \rangle_z}, \qquad (B 14)
$$

455 and by (B 12),

$$
456 \qquad K_0(z) = \left( 1 + \frac{\overline{\epsilon(z)^{2/3}}}{\langle \overline{\epsilon(z)^{2/3}} \rangle_z} \right) \langle K(z) \rangle_z . \qquad (B 15)
$$

457 Since all the terms on the RHS are known, one can evaluate the nonequilibrium kinetic
458 energy profile $K_1(z) = K(z) - K_0(z)$ (see Equation (3.1)).

459 Acknowledgements. We dedicate this work to Jean-Pierre Bertoglio, on the occasion of his retirement. All
460 DNS simulations were carried out using the facilities of the PMCS2I (École Centrale de Lyon). R.A. is
461 supported by the Takenaka Scholarship Foundation.

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

465 Declaration of interests. The authors report no conflict of interest.

REFERENCES

466 BERTOGLIO, JEAN-PIERRE, BATAILLE, FRANÇOISE & MARION, JEAN-DENIS 2001 Two-point closures for
467 weakly compressible turbulence. *Physics of Fluids* **13** (1), 290–310.

468 BERTOGLIO, JEAN-PIERRE & JEANDEL, DENIS 1987 A simplified spectral closure for inhomogeneous turbulence: application to the boundary layer. In *Turbulent Shear Flows* 5, pp. 19–30. Springer.

469 BESNARD, DIDIER C., HARLOW, FRANCIS H., RAUENZAHN, RICK M. & ZEMACH, CHARLES 1996 Spectral transport model for turbulence. *Theoretical and Computational Fluid Dynamics* **8** (1), 1–35.

470 BORUE, VADIM & ORSZAG, STEVEN A. 1996 Numerical study of three-dimensional Kolmogorov flow at high Reynolds numbers. *Journal of Fluid Mechanics* **306**, 293–323.

471 BOS, WOUTER J. T. 2020 Production and dissipation of kinetic energy in grid turbulence. *Physical Review Fluids* **5** (10), 104607.

472 BOS, WOUTER J. T. & RUBINSTEIN, ROBERT 2017 Dissipation in unsteady turbulence. *Physical Review Fluids* **2** (2), 022601.

473 BRIARD, ANTOINE, GRÉA, BENOÎT-JOSEPH, MONS, VINCENT, CAMBON, CLAUDE, GOMEZ, THOMAS & SAGAUT,
479 PIERRE 2018 Advanced spectral anisotropic modelling for shear flows. *Journal of Turbulence* **19** (7), 570–599.

481 CADIOU, ANNE, HANJALIĆ, KEMAL & STAWIARSKI, K 2004 A two-scale second-moment turbulence closure based on weighted spectrum integration. *Theoretical and Computational Fluid Dynamics* **18** (1), 1–26.

484 CAMBON, CLAUDE, JEANDEL, DENIS & MATHIEU, JEAN 1981 Spectral modelling of homogeneous non-isotropic turbulence. *Journal of Fluid Mechanics* **104**, 247–262.

486 DELACHE, ALEXANDRE, CAMBON, CLAUDE & GODEFERD, FABIEN 2014 Scale by scale anisotropy in freely decaying rotating turbulence. *Physics of Fluids* **26** (2), 025104.

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

490 Horiuti, Kiyosi & TAMAKI, TAKAHiro 2013 Nonequilibrium energy spectrum in the subgrid-scale one-equation model in large-eddy simulation. *Physics of Fluids* **25** (12), 125104.

492 Ishihara, Takashi, Yoshida, Kyo & Kaneda, Yukio 2002 Anisotropic velocity correlation spectrum at small scales in a homogeneous turbulent shear flow. *Physical Review Letters* **88** (15), 154501.

494 JEANDEL, DENIS, BRISON, J. F. & MATHIEU, JEAN 1978 Modeling methods in physical and spectral space. *The Physics of Fluids* **21** (2), 169–182.

496 KANEDA, YUKIO & YOSHIDA, Kyo 2004 Small-scale anisotropy in stably stratified turbulence. *New Journal of Physics* **6** (1), 34.

498 KNIGHT, BRUCE & SIROVICH, LAWRENCE 1990 Kolmogorov inertial range for inhomogeneous turbulent flows. *Physical Review Letters* **65** (11), 1356.

500 KOLMOGOROV, ANDREY NIKOLAEVICH 1941 The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. *Doklady Akademii Nauk SSSR* **30**, 299–303.

502 KRAICHNAN, ROBERT H. 1971 An almost-Markovian Galilean-invariant turbulence model. *Journal of Fluid Mechanics* **47** (3), 513–524.

504 KURIEN, SUSAN, L'VOV, VICTOR S., PROCACCIA, ITAMAR & SREENIVASAN, KATEPALLI R. 2000 Scaling structure of the velocity statistics in atmospheric boundary layers. *Physical Review E* **61** (1), 407.

506 LAPORTA, ANDRÉ 1995 Etude spectrale et modélisation de la turbulence inhomogène. PhD thesis, Ecole Centrale de Lyon.

508 LAPORTA, ANDRÉ & BERTOGLIO, JEAN-PIERRE 1995 A model for inhomogeneous turbulence based on two-point correlations. In *Advances in Turbulence V*, pp. 286–297. Springer.

510 LEITH, CECIL E. 1967 Diffusion approximation to inertial energy transfer in isotropic turbulence. *The Physics of Fluids* **10** (7), 1409–1416.

512 LIAO, Zi-Ju & Su, Wei-Dong 2015 Kolmogorov's hypotheses and global energy spectrum of turbulence. *Physics of Fluids* **27** (4), 041701.

514 MOSER, ROBERT D. 1994 Kolmogorov inertial range spectra for inhomogeneous turbulence. *Physics of Fluids* **6** (2), 794–801.

516 MUSACCHIO, STEFANO & BOFFETTA, GUIDO 2014 Turbulent channel without boundaries: The periodic Kolmogorov flow. *Physical Review E* **89** (2), 023004.

518 PARPAIS, SYLVAIN & BERTOGLIO, JEAN-PIERRE 1996 A spectral closure for inhomogeneous turbulence applied to turbulent confined flow. In *Advances in Turbulence VI*, pp. 75–76. Springer.

520 POPE, STEPHEN B. 2000 *Turbulent Flows*. Cambridge University Press.

521 RUBINSTEIN, ROBERT & CLARK, TIMOTHY 2022 *Reassessment of the classical turbulence closures*. Cambridge Scholars Publishing.

523 RUBINSTEIN, ROBERT & CLARK, TIMOTHY T. 2005 Self-similar turbulence evolution and the dissipation rate transport equation. *Physics of Fluids* **17** (9), 095104.

525 TENNEKES, HENDRIK & LUMLEY, JOHN L. 1972 *A First Course in Turbulence*. MIT Press.

526 TOUIL, HATEM, BERTOGLIO, JEAN-PIERRE & SHAO, LIANG 2002 The decay of turbulence in a bounded domain. *Journal of Turbulence* **3** (1), 049.

528 WOODRUFF, STEPHEN L. & RUBINSTEIN, ROBERT 2006 Multiple-scale perturbation analysis of slowly evolving turbulence. *Journal of Fluid Mechanics* **565**, 95–103.

530 WU, WENWEI, SCHMITT, FRANÇOIS G, CALZAVARINI, ENRICO & WANG, LIPO 2021 A quadratic Reynolds stress development for the turbulent Kolmogorov flow. *Physics of Fluids* **33** (12), 125129.

532 YOSHIZAWA, AKIRA 1984 Statistical analysis of the deviation of the Reynolds stress from its eddy-viscosity representation. *The Physics of Fluids* **27** (6), 1377–1387.

535 YOSHIZAWA, AKIRA 1994 Nonequilibrium effect of the turbulent-energy-production process on the inertial-range energy spectrum. *Physical Review E* **49** (5), 4065.
479 PIERRE 2018 Advanced spectral anisotropic modelling for shear flows. *Journal of Turbulence* **19** (7), 570–599.

481 CADIOU, ANNE, HANJALIĆ, KEMAL & STAWIARSKI, K 2004 A two-scale second-moment turbulence closure based on weighted spectrum integration. *Theoretical and Computational Fluid Dynamics* **18** (1), 1–26.

484 CAMBON, CLAUDE, JEANDEL, DENIS & MATHIEU, JEAN 1981 Spectral modelling of homogeneous non-isotropic turbulence. *Journal of Fluid Mechanics* **104**, 247–262.

486 DELACHE, ALEXANDRE, CAMBON, CLAUDE & GODEFERD, FABIEN 2014 Scale by scale anisotropy in freely decaying rotating turbulence. *Physics of Fluids* **26** (2), 025104.

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

490 Horiuti, Kiyosi & Tamaki, TAKAHiro 2013 Nonequilibrium energy spectrum in the subgrid-scale one-equation model in large-eddy simulation. *Physics of Fluids* **25** (12), 125104.

492 Ishihara, Takashi, Yoshida, Kyo & Kaneda, Yukio 2002 Anisotropic velocity correlation spectrum at small scales in a homogeneous turbulent shear flow. *Physical Review Letters* **88** (15), 154501.

494 JEANDEL, DENIS, BRISON, J. F. & MATHIEU, JEAN 1978 Modeling methods in physical and spectral space. *The Physics of Fluids* **21** (2), 169–182.

496 KANEDA, YUKIO & YOSHIDA, Kyo 2004 Small-scale anisotropy in stably stratified turbulence. *New Journal of Physics* **6** (1), 34.

498 KNIGHT, BRUCE & SIROVICH, LAWRENCE 1990 Kolmogorov inertial range for inhomogeneous turbulent flows. *Physical Review Letters* **65** (11), 1356.

500 KOLMOGOROV, ANDREY NIKOLAEVICH 1941 The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. *Doklady Akademii Nauk SSSR* **30**, 299–303.

502 KRAICHNAN, ROBERT H. 1971 An almost-Markovian Galilean-invariant turbulence model. *Journal of Fluid Mechanics* **47** (3), 513–524.

504 KURIEN, SUSAN, L'VOV, VICTOR S., PROCACCIA, ITAMAR & SREENIVASAN, KATEPALLI R. 2000 Scaling structure of the velocity statistics in atmospheric boundary layers. *Physical Review E* **61** (1), 407.

506 LAPORTA, ANDRÉ 1995 Etude spectrale et modélisation de la turbulence inhomogène. PhD thesis, Ecole Centrale de Lyon.

508 LAPORTA, ANDRÉ & BERTOGLIO, JEAN-PIERRE 1995 A model for inhomogeneous turbulence based on two-point correlations. In *Advances in Turbulence V*, pp. 286–297. Springer.

510 LEITH, CECIL E. 1967 Diffusion approximation to inertial energy transfer in isotropic turbulence. *The Physics of Fluids* **10** (7), 1409–1416.

512 LIAO, Zi-Ju & Su, Wei-Dong 2015 Kolmogorov's hypotheses and global energy spectrum of turbulence. *Physics of Fluids* **27** (4), 041701.

514 MOSER, ROBERT D. 1994 Kolmogorov inertial range spectra for inhomogeneous turbulence. *Physics of Fluids* **6** (2), 794–801.

516 MUSACCHIO, STEFANO & BOFFETTA, GUIDO 2014 Turbulent channel without boundaries: The periodic Kolmogorov flow. *Physical Review E* **89** (2), 023004.

518 PARPAIS, SYLVAIN & BERTOGLIO, JEAN-PIERRE 1996 A spectral closure for inhomogeneous turbulence applied to turbulent confined flow. In *Advances in Turbulence VI*, pp. 75–76. Springer.

520 POPE, STEPHEN B. 2000 *Turbulent Flows*. Cambridge University Press.

521 RUBINSTEIN, ROBERT & CLARK, TIMOTHY 2022 *Reassessment of the classical turbulence closures*. Cambridge Scholars Publishing.

523 RUBINSTEIN, ROBERT & CLARK, TIMOTHY T. 2005 Self-similar turbulence evolution and the dissipation rate transport equation. *Physics of Fluids* **17** (9), 095104.

525 TENNEKES, HENDRIK & LUMLEY, JOHN L. 1972 *A First Course in Turbulence*. MIT Press.

526 TOUIL, HATEM, BERTOGLIO, JEAN-PIERRE & SHAO, LIANG 2002 The decay of turbulence in a bounded domain. *Journal of Turbulence* **3** (1), 049.

528 WOODRUFF, STEPHEN L. & RUBINSTEIN, ROBERT 2006 Multiple-scale perturbation analysis of slowly evolving turbulence. *Journal of Fluid Mechanics* **565**, 95–103.

530 WU, WENWEI, SCHMITT, FRANÇOIS G, CALZAVARINI, ENRICO & WANG, LIPO 2021 A quadratic Reynolds stress development for the turbulent Kolmogorov flow. *Physics of Fluids* **33** (12), 125129.

532 YOSHIZAWA, AKIRA 1984 Statistical analysis of the deviation of the Reynolds stress from its eddy-viscosity representation. *The Physics of Fluids* **27** (6), 1377–1387.

535 YOSHIZAWA, AKIRA 1994 Nonequilibrium effect of the turbulent-energy-production process on the inertial-range energy spectrum. *Physical Review E* **49** (5), 4065.