2.1. Basic ingredients of the phenomenology of fluid turbulence and the implied statistical spatial structure

Let us consider the velocity vector field ${\boldsymbol{u}}^\nu (t,{\boldsymbol{x}})\in {\mathbb{R}}^3$ for ${\boldsymbol{x}}\in {\mathbb{R}}^3$ and $t\geqslant 0$ , which reads componentwise ${\boldsymbol{u}}^\nu =(u_i^\nu )_{1\leqslant i \leqslant 3}$ , solution of the Navier–Stokes equations, for incompressible fluids of unit density:

(2.1) \begin{align} \frac {\partial {\boldsymbol{u}}^\nu }{\partial t}+({\boldsymbol{u}}^\nu \boldsymbol{\cdot }\boldsymbol{\nabla }) {\boldsymbol{u}}^\nu = -\boldsymbol{\nabla }\!p^\nu +\nu \Delta {\boldsymbol{u}}^\nu +\boldsymbol{f}. \end{align}

Here $\nu$ is the kinematic viscosity, $p^\nu (t,{\boldsymbol{x}})$ the pressure field obtained by the Poisson equation using the incompressibility constraint $\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{u}}^\nu =0$ and $\boldsymbol{f}(t,{\boldsymbol{x}})$ a forcing vector assumed to be divergence free, smooth in space, correlated on a characteristic length scale $L$ of the order of the integral length scale and independent of viscosity. We will denote the velocity gradient $\partial _{\!j} u_i^\nu$ of the component $u_i^\nu$ along $x_{\!j}$ .

For quite a general form of the forcing field $\boldsymbol{f}$ , observations suggest that the velocity field ${\boldsymbol{u}}^\nu$ reaches a statistically homogeneous and stationary regime, in which the velocity variance and average dissipation depend very weakly on viscosity, such that

(2.2) \begin{align} \sigma ^2= \lim _{\nu \to 0} {\mathbb{E}}\big [|{\boldsymbol{u}}^\nu |^2\big ] \end{align}

and

(2.3) \begin{align} \varepsilon =\lim _{\nu \to 0} \lim _{t\to \infty }\nu \sum _{\textit{ij}}{\mathbb{E}} \big [\big ( \partial _{\!j} u^\nu _i \big )^2\big ] \end{align}

are positive numbers. Starting from the randomly forced Navier–Stokes equations (2.1), the expectations in (2.2) and (2.3) are taken over the forcing realisations. From a practical point of view, we assume that expectations can be accurately replaced by empirical averages over space and/or time, which is in particular guaranteed by the statistical homogeneity and stationarity observed in DNS and experiments. In the asymptotic limit of infinite Reynolds numbers, i.e. in the vanishing viscosity limit $\nu \to 0$ , the finite quantities $\sigma$ (2.2) and $\varepsilon$ (2.3) are expected to be related according to the dimensional prediction (see the classical textbooks by Frisch (Reference Frisch1995) and Pope (Reference Pope2000) for instance)

(2.4) \begin{align} \varepsilon \propto \frac {\sigma ^3}{L}. \end{align}

The independence of the velocity variance on viscosity (2.2) is far from being obvious and is at the heart of the phenomenology of fluid turbulence. Recall that energy is injected into the fluid in a statistically stationary way through the forcing term $\boldsymbol{f}$ (2.1). In order to maintain a velocity variance that becomes independent of viscosity (2.2) as $\nu \to 0$ , the fluid develops small scales by populating Fourier modes located at higher wavenumbers than those excited by the forcing.

To characterise this behaviour, usually referred to as the cascade phenomenon, the PSD $E_\nu ^{\textit{E}}({\boldsymbol{k}})$ has been measured and estimated in many studies. This quantity is defined to be the Fourier transform of the trace of the covariance tensor of the velocity vector field, that is (using Einstein’s summation convention over repeated indices),

(2.5) \begin{align} E_\nu ^{\textit{E}}({\boldsymbol{k}}) = \int _{{\boldsymbol{r}}\in {\mathbb{R}}^3} \textit{e}^{-2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{r}}}{\mathbb{E}} \left [ u^\nu _i(t,{\boldsymbol{x}})u^\nu _i(t,{\boldsymbol{x}}+{\boldsymbol{r}}) \right ]{\mathrm{d}} {\boldsymbol{r}}. \end{align}

Note that the correlation $ {\mathbb{E}} [ u_i(t,{\boldsymbol{x}})u_i(t,{\boldsymbol{x}}+{\boldsymbol{r}}) ]$ is independent of both position $\boldsymbol{x}$ and of time $t$ due to statistical homogeneity and stationarity, respectively. Furthermore, by virtue of statistical isotropy, $E_\nu ^{\textit{E}}({\boldsymbol{k}})=E_\nu ^{\textit{E}}(|{\boldsymbol{k}}|)$ , i.e. the PSD is a function of the norm $|{\boldsymbol{k}}|$ only, and is independent of the direction of $\boldsymbol{k}$ . Note that the upper letter ` $^{\textit{E}}$ ’ indicates that $E_\nu ^{\textit{E}}({\boldsymbol{k}})$ is a genuinely spatial, i.e. Eulerian, quantity. As shown by Pope (Reference Pope2000), the PSD $E_\nu ^{\textit{E}}({\boldsymbol{k}})$ defined in (2.5) fully determines the covariance function of the velocity field owing to the identity

(2.6) \begin{align} {\mathbb{E}} \big [ u^\nu _i(t,{\boldsymbol{x}})u^\nu _{\!j}(t,{\boldsymbol{x}}+{\boldsymbol{r}}) \big ] = \frac {1}{2}\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3} \textit{e}^{2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{r}}} E_\nu ^{\textit{E}}(|{\boldsymbol{k}}|) \widehat {P}_{\textit{ij}}({\boldsymbol{k}} ){\mathrm{d}} {\boldsymbol{k}} , \end{align}

where Leray’s projector on divergence-free vector fields $ \widehat {P}_{\textit{ij}}({\boldsymbol{k}} )$ has been used and is defined by

(2.7) \begin{align} \widehat {P}_{\textit{ij}}({\boldsymbol{k}} ) = \delta _{\textit{ij}} - \frac {k_ik_{\!j}}{|{\boldsymbol{k}}|^2}, \end{align}

with $\delta _{\textit{ij}}$ the Kronecker delta symbol.

Although being accessible in DNS of the Navier–Stokes equations, the 3 $d$ PSD $E_\nu ^{\textit{E}}({\boldsymbol{k}})$ (2.5) is difficult to obtain from experimental measurements. Instead, one-dimensional longitudinal PSD has traditionally been estimated. By virtue of statistical isotropy, this quantity is defined as the one-dimensional Fourier transform of the correlation function of one of the components of velocity, say $u_1$ , along the direction $\hat {e}_1$ , that is,

(2.8) \begin{align} E_\nu ^{\textit{E,long}}(k) = \int _{\ell \in {\mathbb{R}}} \textit{e}^{-2i\pi k\ell }{\mathbb{E}} \left [ u^\nu _1(t,{\boldsymbol{x}})u^\nu _1(t,{\boldsymbol{x}}+\ell \hat {e}_1) \right ]{\mathrm{d}}\ell . \end{align}

It can be checked that $E_\nu ^{\textit{E,long}}(k)$ is a positive and even function of $k$ . The three-dimensional PSD $E_\nu ^{\textit{E}}({\boldsymbol{k}})$ (2.5), which is a function of the norm $k = |{\boldsymbol{k}}|$ of the wave vector $\boldsymbol{k}$ , and the one-dimensional longitudinal PSD $E_\nu ^{\textit{E,long}}(k)$ (2.8) are related as

(2.9) \begin{align} E_\nu ^{\textit{E}}(k) = \frac {k}{2\pi }\frac {{\mathrm{d}}}{{\mathrm{d}} k}\left (\frac {1}{k}\frac {{\mathrm{d}} E_\nu ^{\textit{E,long}}(k)}{{\mathrm{d}} k }\right )\!. \end{align}

Equation (2.9) is a consequence of incompressibility and statistical isotropy. A complete proof of relation (2.9) is provided by Pope (Reference Pope2000).

Going back to the aforementioned phenomenological aspects of fluid turbulence, a manifestation of the underlying cascading process of energy can be read on the PSD, which will follow, in the limit of infinite Reynolds numbers, a power-law behaviour at large wavenumbers. This is known as the Kolmogorov ‘ $k^{-5/3}$ ’ spectrum (Frisch Reference Frisch1995). At a given finite Reynolds number, a realistic parameterisation of the longitudinal PSD $E_\nu ^{\textit{E,long}}(k)$ (2.8) would require a regularisation at small wavenumbers of order $1/L$ , i.e. around the inverse integral length scale stemming from the forcing term, and a dissipative cutoff at large wavenumbers $1/\eta _K$ , where $\eta _K \propto (\nu ^3/\varepsilon )^{1/4}$ is known as Kolmogorov’s dissipative length scale. Sophisticated models that do so have been proposed in the literature (Pope Reference Pope2000; Meyers & Meneveau Reference Meyers and Meneveau2008). For the sake of simplicity, we will be working with the following simple but realistic parameterisation:

(2.10) \begin{align} E_\nu ^{\textit{E,long}}(k) =D_2 |k|_{1/L}^{-5/3} \textit{e}^{-\eta _d k}.\end{align}

Here we have introduced a regularised norm $|k|_{1/L} = \sqrt {k^2 + 1/L^2}$ , a multiplicative constant $D_2\gt 0$ that is related to velocity variance (2.2) according to the identity

(2.11) \begin{align} D_2 = \frac {\sigma ^2}{3L^{2/3}}\frac {\varGamma (5/6)}{\sqrt {\pi }\varGamma (1/3)} \end{align}

and

(2.12) \begin{align} \eta _d = 2\pi \left ( \frac {10C_2\varGamma (4/3)}{\varGamma (1/3)}\right )^{3/4}\left (\frac {\nu ^3}{\varepsilon }\right )^{1/4}\!, \end{align}

which is of order of the Kolmogorov dissipative length scale $\eta _K$ . $\varGamma ( z )$ denotes here the gamma function. The expression of the prefactor in the definition of $\eta _d$ (2.12) is derived below, and is chosen such that to obtain ultimately (2.22). The so-called Kolmogorov constant $C_2$ has been estimated on experimental data and is close to $C_2 \approx 2$ ; see Pope (Reference Pope2000). From the proposed expression for $E_\nu ^{\textit{E,long}}(k)$ (2.10), we deduce from (2.9) the expression of the three-dimensional PSD $E_\nu ^{\textit{E}}(k)$ (2.5), which has a particularly simple expression in the vanishing viscosity limit, i.e.

(2.13) \begin{align} E^{\textit{E}}(k) =\lim _{\nu \to 0}E_\nu ^{\textit{E}}(k) = \frac {1}{2\pi }\frac {55}{9} D_2 k^2 |{\boldsymbol{k}}|_{1/L}^{-17/3}. \end{align}

Note that the proportionality constant appearing in the expression of $D_2$ (2.11) ensures that the velocity variance $\sigma ^2$ (2.2) is given by the integral of the PSD (2.10). We therefore have the identity

(2.14) \begin{align} \sigma _\nu ^2=\lim _{t\to \infty }{\mathbb{E}}\big [|{\boldsymbol{u}}^\nu (t,{\boldsymbol{x}})|^2\big ]=\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3} E_\nu ^{\textit{E}}({\boldsymbol{k}}) {\mathrm{d}} {\boldsymbol{k}}=3\int _{k\in {\mathbb{R}}} E_\nu ^{\textit{E,long}}(k) {\mathrm{d}} k . \end{align}

Consequently, using the limiting value of $\sigma _\nu ^2$ as $\nu \to 0$ as expressed in (2.14) and the parameterisation of the longitudinal PSD provided in (2.10) and recalling that $\eta _d$ vanishes as $\nu \to 0$ (2.12), we obtain

(2.15) \begin{align} \sigma ^2=\lim _{\nu \to 0}\sigma _\nu ^2=6D_2\int _{k=0}^{\infty } |k|_{1/L}^{-5/3} {\mathrm{d}} k , \end{align}

where $D_2$ is given by (2.11). Also, as claimed above, the regularisation parameter $L$ appearing in the expression of the longitudinal PSD (2.10) can be interpreted as the correlation length scale of the longitudinal component. More precisely, $L$ can be related to the so-called integral length scale $L_{\textit{int}}$ , i.e. the integral of the correlation function, as

(2.16) \begin{align} L_{\textit{int}}\equiv \int _{\ell \gt 0}\frac {3}{\sigma ^2}{\mathbb{E}} \left [ u^\nu _1(t,{\boldsymbol{x}})u^\nu _1(t,{\boldsymbol{x}}+\ell \hat {e}_1) \right ]{\rm d}\ell = \frac {3}{2\sigma ^2}E_\nu ^{\textit{E,long}}(0) = \frac {\varGamma (5/6)}{2\sqrt {\pi }\varGamma (1/3)}L. \end{align}

Concerning the average dissipation $\varepsilon _\nu$ , by making use of the isotropy condition, we have

(2.17) \begin{align} \varepsilon _\nu =15\nu {\mathbb{E}} \big [\left ( \partial _1 u_1 \right )^2\big ]&=15\nu \int _{k\in {\mathbb{R}}} (2\pi )^2k^2E_\nu ^{\textit{E,long}}(k) {\mathrm{d}} k\notag \\ &=30(2\pi )^2D_2\nu \eta _d^{-4/3}\int _{k=0}^{\infty } k^2|k|_{\eta _d/L}^{-5/3} \textit{e}^{-k}{\mathrm{d}} k\notag \\ &\mathrel {\mathop {\sim }\limits _{\nu \to 0}}30(2\pi )^2D_2\nu \eta _d^{-4/3}\varGamma (4/3). \end{align}

Taking into account the expressions of $D_2$ (2.11) and $\eta _d$ (2.12), we then obtain

(2.18) \begin{align} \varepsilon =\lim _{\nu \to 0}\varepsilon _\nu = 2\pi \left ( \frac {\varGamma (5/6)}{\sqrt {\pi }C_2}\right )^{3/2}\frac {\sigma ^3}{L}. \end{align}

In the vanishing viscosity limit, a Gaussian random field that possesses a PSD that decays as a power law of the type proposed in (2.10) is known in the literature as a fractional Gaussian field of Hölder (or Hurst) parameter $H=1/3$ . In particular, such a field is continuous and its variance is finite (2.15), but it is nowhere differentiable, and must be regarded as being Hölder continuous. To characterise the Hölder continuity, instead of using gradients that are expected to be infinite, one needs to consider velocity increments that are well posed. Precisely, the longitudinal velocity increment $\delta _{\boldsymbol{\ell }}^{\textit{long}} {\boldsymbol{u}}^\nu$ is defined by

(2.19) \begin{align} \delta _\ell ^{\textit{long}} {\boldsymbol{u}}^\nu (t,{\boldsymbol{x}}) = \left ({\boldsymbol{u}}^\nu (t,{\boldsymbol{x}}+{\boldsymbol{\ell }}) - {\boldsymbol{u}}^\nu (t,{\boldsymbol{x}})\right )\boldsymbol{\cdot }\frac {{\boldsymbol{\ell }}}{\ell } \end{align}

for all ${\boldsymbol{\ell }} \neq 0$ . The variance of the longitudinal increment (2.19) is known as the longitudinal second-order structure function and is given by

(2.20) \begin{align} S_{2,\nu }^{\textit{long}}(\ell )={\mathbb{E}} \big [ \big (\delta _\ell ^{\textit{long}} {\boldsymbol{u}}^\nu \big )^2\big ]. \end{align}

In the vanishing viscosity limit, weobtain

(2.21) \begin{align} S_2^{\textit{long}}(\ell ) &= \lim _{\nu \to 0} {\mathbb{E}} \big [ \big (\delta _\ell ^{\textit{long}} {\boldsymbol{u}}^\nu \big )^2\big ] \notag \\ &= 2\int _{k\in {\mathbb{R}}} \big ( 1-\textit{e}^{2i\pi k\ell }\big ) E^{\textit{E,long}}(k) {\mathrm{d}} k = 4D_2\int _{k=0}^{\infty } \left ( 1-\cos (2\pi k\ell )\right ) |k|_{1/L}^{-5/3} {\mathrm{d}} k\notag \\ &\mathrel {\mathop {\sim }\limits _{\ell \to 0}}4D_2 \ell ^{2/3}\int _{k=0}^{\infty } \left ( 1-\cos (2\pi k)\right ) k^{-5/3} {\mathrm{d}} k=3D_2(2\pi )^{2/3}\varGamma (1/3)\ell ^{2/3}. \end{align}

As can be seen from the behaviour of $S_2^{\textit{long}}(\ell )$ (2.21) on small scales $\ell \to 0$ , the second-order structure function exhibits a power-law behaviour with respect to the scale $\ell$ with the exponent $2/3$ . Then using the expression of $D_2$ (2.11) and expressing the remaining multiplicative constant in units of the average dissipation $\varepsilon$ (2.18), we finally get

(2.22) \begin{align} S_2^{\textit{long}}(\ell ) \mathrel {\mathop {\sim }\limits _{\ell \to 0}} C_2|\varepsilon \ell |^{2/3}, \end{align}

which justifies the introduction of the Kolmogorov constant $C_2$ in the definition of the dissipative cutoff $\eta _d$ (2.12).

More generally, let us go back to the full vector field. Its covariance is provided in (2.6). Let us introduce the increment of a given velocity component $u_i^\nu$ ,

(2.23) \begin{align} \delta _{\ell }u_i^\nu (t,{\boldsymbol{x}}) = u_i^\nu (t,{\boldsymbol{x}}+{\boldsymbol{\ell }})-u_i^\nu (t,{\boldsymbol{x}}), \end{align}

and the respective second-order structure function

(2.24) \begin{align} S_{2,ij}(\ell ) = \lim _{\nu \to 0}{\mathbb{E}}\big [\delta _{\ell }u_i^\nu (t,{\boldsymbol{x}}) \delta _{\ell }u_{\!j}^\nu (t,{\boldsymbol{x}}) \big ]. \end{align}

Then, under the assumption of statistical isotropy and incompressibility, from (Pope Reference Pope2000), we obtain

(2.25) \begin{align} S_{2,ij}(\ell ) \mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}} C_2|\varepsilon \ell |^{2/3}\left [-\frac {1}{3}\frac {\ell _i\ell _{\!j}}{\ell ^2}+\frac {4}{3}\delta _{\textit{ij}} \right ]\!. \end{align}

As a consequence, we have

(2.26) \begin{align} S_2(\ell )\equiv S_{2,ii}(\ell ) = \lim _{\nu \to 0} {\mathbb{E}} \big [ \big |\delta _\ell {\boldsymbol{u}}^\nu \big |^2\big ] &= 2\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3} \big ( 1-\textit{e}^{2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}}\big ) E^{\textit{E}}(|{\boldsymbol{k}}|) {\mathrm{d}} {\boldsymbol{k}}\notag \\ &\mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}}\frac {11}{3}C_2|\varepsilon \ell |^{2/3}. \end{align}

Firstly taking the limit $\ell \to 0$ and secondly the limit $\eta _d\to 0$ (or equivalently $\nu \to 0$ ), we then obtain the power-law behaviour

(2.27) \begin{align} {\mathbb{E}} \big [ \big (\delta _\ell ^{\textit{long}} {\boldsymbol{u}}^\nu \big )^2\big ] \mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}} 2\ell ^2\int _{k\in {\mathbb{R}}} \left (2\pi k \right )^2 E^{\textit{E,long}}_\nu (k) {\mathrm{d}} k\mathrel {\mathop {\sim }\limits _{\nu \to 0}} \frac {\varepsilon _\nu }{15\nu } \ell ^2, \end{align}

where the expression of $\varepsilon _\nu$ is given in (2.17). The power law (2.27) is associated with differentiable velocity fields instead of being Hölder continuous with fractional exponent $H=1/3$ .

2.2. A first comparison concerning the spatial structure of DNS and the Gaussian model

We show in figure 1 a comparison between the velocity field ${\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})$ (figure 1 a) extracted from a DNS of the forced Navier–Stokes equations (2.1) and a Gaussian model (figure 1b). Concerning the DNS, it is provided by the Johns Hopkins turbulence database (Li et al. Reference Li, Perlman, Wan, Yang, Burns, Meneveau, Burns, Chen, Szalay and Eyink2008). The data set concerns a simulation of the Navier–Stokes equations in a homogeneous and isotropic situation, using a pseudo-spectral method fully dealiased with $N=1024$ nodes in each direction, over 5028 time steps, corresponding to roughly speaking ten turnover time scales. The Taylor-based Reynolds number is estimated to be 433, corresponding to a fully developed turbulent state. Detailed properties of the simulation are provided on their website. For our purpose, we have downloaded the three velocity components over three spatial slices , $(0,y,z)$ , $(x,0,z)$ and $(x,y,0)$ , at all available times. The Gaussian model ${\boldsymbol{u}}^{\textit{fGv},\nu }({\boldsymbol{x}})$ corresponds to a fractional Gaussian vector (fGv) field, in a periodic setting, presented in § 3.2.1. At this stage, it is enough to state that this is the unique statistically homogeneous, isotropic and incompressible Gaussian vector field consistent with the longitudinal PSD $E_\nu ^{\textit{E,long}}(k)$ provided in (2.10). Following a suitable truncation of the Fourier modes, as described in § 4.1, we have similarly performed a simulation of this Gaussian vector field on a periodic box of size $L_{\textit{tot}}=2\pi$ using $N=1024$ nodes in each direction, with the remaining free parameters entering in (2.10) chosen to match the DNS dataset. This corresponds to $D_2 = 0.021$ , $L= L_{\textit{tot}}$ and $\eta _d = 0.085$ . The results of our comparison are collected in figure 1.

Figure 1. Comparison of the instantaneous and statistical spatial structure of a DNS velocity field ${\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})$ and of a realisation of the model that coincides with a fractional Gaussian vector field at a given instant of time in the statistically stationary regime. (a) Instantaneous spatial profile of the norm of the DNS velocity field $|{\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})|$ in the plane $y=0$ , at the initial time of the DNS dataset. (b) Same as (a) for the model, we used the same colourbar in both representations. (c) Estimation of the longitudinal PSDs $E_\nu ^{\textit{E,long}}(k)$ (2.8) based on the variance of the Fourier modes of the one-dimensional discrete Fourier transform. Statistics (DNS in yellow and model in blue) are estimated by averaging both in space and time. We superimpose with a black dashed line the inviscid limit of the functional form provided in (2.10), corresponding to the power law $D_2|{\boldsymbol{k}}|^{-5/3}$ . (d) Statistical estimation of the second-order longitudinal structure function $S_{2,\nu }^{\textit{long}}(\ell )$ (2.20); same colours as (c). We superimpose with dashed lines the inviscid predictions in the three ranges of scales of interest: (i) at large scales of order $L_{\textit{tot}}$ , $S_2^{\textit{long}}(\ell )$ reaches the plateau $({2}/{3})\sigma ^2$ , where $\sigma ^2$ is related to the free parameter $D_2$ according to (2.11); (ii) in the intermediate inertial range of scales, we represent the prediction made in (2.22); and (iii) by the smooth behaviour $\propto \ell ^2$ in the dissipative range, as predicted in (2.27).

We begin by displaying in figure 1(a) a snapshot of the norm of the DNS velocity field $|{\boldsymbol{u}}^{\textit{DNS},\nu }(t = 0,x,y=0,z)|$ in the plane defined by the Cartesian coordinates $y=0$ at the initial time of the dataset, which has been chosen in the statistically stationary regime. Similarly, in figure 1(b) we display for comparison the same snapshot but for the model $|{\boldsymbol{u}}^\nu (t=0,x,y=0,z)|$ . The colourbars are the same in both panels, where dark corresponding to the lowest values and light corresponding to the highest values of the norm. We observe similar trends concerning the amplitude of the fluctuations and also the statistical symmetries, i.e. homogeneity and isotropy. However, it appears that, while ${\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})$ exhibits filamentary structures, the Gaussian field looks more patchy. Despite their different appearances, we will show that these two vector fields have consistent second-order statistics.

The statistical analysis of these two vector fields is shown in figures 1(c) and 1(d). The respective longitudinal PSDs $E_\nu ^{\textit{E,long}}(k)$ are shown in panel (c). Concerning the DNS, it is obtained as the variance of the respective Fourier modes obtained by performing one-dimensional discrete Fourier transforms (DFTs) of the longitudinal components, say, for instance, the $x$ component of velocity along the $x$ direction. The statistical sample is averaged over the three longitudinal components, spatially over the slices, and over time. Results of our estimation are displayed with an yellow solid line. The estimation of the Gaussian model statistics is also performed by averaging over the three components, over space and time, as has been done for the DNS. We display the result of our estimation using a solid blue line. As we can observe, both longitudinal PSDs superimpose in a very satisfactory manner, showing that, based on the variance of Fourier modes, these two fields are statistically almost indiscernible. We mention that we could have, instead of the estimate on the Gaussian model, displayed the theoretical functional form provided in (2.10) without a noticeable difference, as expected. The result of our statistical estimation of the longitudinal structure function $S_{2,\nu }^{\textit{long}}(\ell )$ (2.20) is displayed in figure 1(d), using an yellow line for DNS and a blue line (dark grey) for the Gaussian model, averages have been taken on the same statistical sample as the one considered for figure 1(c). The three ranges of scales, including decorrelation length scales of the order of the size of the box $L_{\textit{tot}}$ , inertial and dissipative length scales, are reproduced in a very satisfactory manner by the model, as expected from the behaviours of the PSDs in figure 1(c).

We conclude this section noting that, whereas it is clear that instantaneous snapshots of DNS and spatial representations of the model, as they are displayed in figures 1(a) and 1(b), exhibit differences regarding in particular filamentary structures, they are nonetheless barely indiscernible from a second-order statistical point of view. For this reason, we claim that a model based on a Gaussian approximation is a fairly good starting point. We are planning to improve this in future investigations.

2.3. The statistical temporal structure of the turbulent velocity field and the sweeping effect

Fluctuations in time of turbulent velocity can be considered counter-intuitive from a dimensional point of view (Tennekes & Lumley Reference Tennekes and Lumley1972; Tennekes Reference Tennekes1975). In particular, the sweeping effect that can be seen, loosely speaking, as the advection of small scales by large scales, has been investigated and defined in very different ways by several groups (Kraichnan Reference Kraichnan1964; Belinicher & L’vov Reference Belinicher and L’vov1987; Gotoh et al. Reference Gotoh, Rogallo, Herring and Kraichnan1993; Kaneda et al. Reference Kaneda, Ishihara and Gotoh1999; Chaves et al. Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003; Chevillard et al. Reference Chevillard, Roux, Lévêque, Mordant, Pinton and Arnéodo2005; Favier et al. Reference Favier, Godeferd and Cambon2010; Biferale, Calzavarini & Toschi Reference Biferale, Calzavarini and Toschi2011; Drivas et al. Reference Drivas, Johnson, Lalescu and Wilczek2017; Canet et al. Reference Canet, Rossetto, Wschebor and Balarac2017; Gorbunova et al. Reference Gorbunova, Balarac, Canet, Eyink and Rossetto2021; Matsumoto et al. Reference Matsumoto, Otsuki, Ooshida and Goto2021). We should bear in mind that such a mechanism should not be confused with Taylor’s hypothesis related to the presence of a mean flow (Wilczek & Narita Reference Wilczek and Narita2012; He, Jin & Yang Reference He, Jin and Yang2017). As far as we are concerned, let us discuss this phenomenon in the spirit of Tennekes (Reference Tennekes1975) who focuses on the dimensional aspects of this phenomenon and its statistical signature.

From the statistical spatial structure of the velocity field formerly presented, in a large part given by the scaling relation provided in (2.22) and (2.26), let us infer the statistical temporal structure of the velocity field ${\boldsymbol{u}}^\nu$ in the asymptotic limit $\nu \to 0$ . The Hölder regularity pinpointed in (2.22) and (2.26) suggests that eddies of characteristic length scale $\ell$ have a characteristic turnover time scale $\tau _{\text{e}}(\ell )$ given by

(2.28) \begin{align} \tau _{\text{e}}(\ell ) \equiv \frac {\ell }{\sqrt {S_2(\ell )}}\mathrel {\mathop {\propto }\limits _{\ell \to 0}^{}} \frac {\ell ^{2/3}}{\sqrt {C_2}\varepsilon ^{1/3}}. \end{align}

The square root of the structure function $\sqrt {S_2(\ell )}$ can be interpreted as the typical velocity of eddies of size $\ell$ . The characteristic turnover time $\tau _{\text{e}}(\ell )$ (2.28) of eddies of typical size $\ell$ would then follow a power law with exponent $2/3$ as a function of the length scale $\ell$ . One may wonder whether this prediction is correct. The proposed definition of the time scale $\tau _{\text{e}}(\ell )$ (2.28) clearly depends solely on the spatial structure of the velocity field and is independent of the temporal structure of the velocity field. In order to take into account the temporal structure, let us thus define a characteristic time scale $\tau _c(\ell )$ , which is based on the correlation in time of the spatial increment given by

(2.29) \begin{align} \tau _{\textit{c}}(\ell ) \equiv \frac {1}{S_2(\ell )} \lim _{\nu \to 0} \int _{0}^{\infty }{\mathbb{E}} \left [\delta _\ell {\boldsymbol{u}}^\nu (t,{\boldsymbol{x}})\boldsymbol{\cdot }\delta _\ell {\boldsymbol{u}}^\nu (t+\tau ,{\boldsymbol{x}}) \right ]{\mathrm{d}} \tau . \end{align}

As it will be argued below, we could also define a similar time scale as $\tau _{\textit{c}}(\ell )$ (2.29) to characterise the turnover time scale of an eddy of size $\ell$ using instead the temporal structure of the velocity Fourier modes at a given wave vector $k$ . This indicates that one associates in a loose sense the length scale $\ell$ to $|{\boldsymbol{k}}|^{-1}$ . This must be treated with great care since we have assumed in this section that the spatial domain is infinite, i.e. ${\boldsymbol{x}}\in {\mathbb{R}}^3$ , and thus, Fourier modes must be considered as distributions. Nevertheless, we would define in a formal way the time scale

(2.30) \begin{align} \tau _{\textit{f}}(k) \equiv \text{`}\lim _{\nu \to 0} \frac {1}{{\mathbb{E}} \big [\left |\widehat {{\boldsymbol{u}}}^\nu (t,{\boldsymbol{k}})\right |^2\big ]}\int _{0}^{\infty }{\mathbb{E}} \left [\widehat {{\boldsymbol{u}}}^\nu (t,{\boldsymbol{k}})\boldsymbol{\cdot }\overline {\widehat {{\boldsymbol{u}}}^\nu (t+\tau ,{\boldsymbol{k}})} \right ]{\mathrm{d}} \tau \text{'}, \end{align}

where $\widehat {{\boldsymbol{u}}}^\nu (t,{\boldsymbol{k}})$ is the Fourier transform of the velocity field ${\boldsymbol{u}}^\nu (t,{\boldsymbol{x}})$ and $\overline {\boldsymbol{\cdot }}$ denotes the complex conjugation. The definition of the time scale $\tau _{\textit{f}}(k)$ (2.30) is for the time being put in quotes because, as we will see, the Fourier transform over all space has a distributional nature and necessitates the formulation of the model in a periodic setting.

Note that the time scale $\tau _{\textit{f}}(k)$ is expected to depend solely on the norm $k$ of the wave vector $\boldsymbol{k}$ by statistical isotropy and can be interpreted as the typical correlation duration at a given wavenumber. The two time scales $\tau _{\textit{c}}(\ell )$ (2.29) and $\tau _{\textit{f}}(k)$ (2.30) are related according to

(2.31) \begin{align} \tau _{\textit{c}}(\ell ) =\text{`}\frac {\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3}\left [ 1-\cos \left ( 2\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}\right )\right ]\tau _{\textit{f}}(k){\mathbb{E}} \big [\left |\widehat {{\boldsymbol{u}}}(t,{\boldsymbol{k}})\right |^2\big ]{\mathrm{d}}{\boldsymbol{k}}}{\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3}\left [ 1-\cos \left ( 2\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}\right )\right ]{\mathbb{E}}\big [ \left |\widehat {{\boldsymbol{u}}}(t,{\boldsymbol{k}})\right |^2\big ]{\mathrm{d}}{\boldsymbol{k}}}\text{'}, \end{align}

where we have denoted ${\mathbb{E}} [ \left |\widehat {{\boldsymbol{u}}}(t,{\boldsymbol{k}})\right |^2 ]$ the limit as $\nu \to 0$ of the PSD ${\mathbb{E}} [ \left |\widehat {{\boldsymbol{u}}}^\nu (t,{\boldsymbol{k}})\right |^2 ]$ . Once again, we recall that expressions such as (2.31) should be taken in a formal way. Only the consideration of periodic boundary conditions gives a proper meaning to such quantities.

With its definition (2.30), the characteristic time scale $\tau _{\textit{f}}(k)$ has been estimated in DNS by several groups (Kaneda et al. Reference Kaneda, Ishihara and Gotoh1999; Favier et al. Reference Favier, Godeferd and Cambon2010; Gorbunova et al. Reference Gorbunova, Balarac, Canet, Eyink and Rossetto2021) while assuming periodic boundary conditions instead of the full space ${\mathbb{R}}^3$ , leading to well-defined Fourier modes. The numerical investigations from the references quoted above show, in a consistent manner with the dimensional analysis of Tennekes (Reference Tennekes1975), that

(2.32) \begin{align} \tau _{\textit{f}}(k) \mathrel {\mathop {\propto }\limits _{k\to \infty }^{}} \frac {1}{\sigma k}. \end{align}

The behaviour of the time scale $\tau _{\textit{f}}(k)$ (2.32) at large wavenumbers as $k^{-1}$ suggests, based on the dimensional ground, that the typical velocity entering in the motion of eddies of size $k^{-1}$ is the standard deviation of the velocity $\sigma$ , i.e. the characteristic velocity of large scales. To see this more clearly, let us see the implications of this behaviour on the time scale $\tau _{\textit{c}}(\ell )$ (2.29). To do so, consider the identity displayed in (2.31) and rescale the dummy integration variable $\boldsymbol{k}$ by $\ell$ , then we obtain

(2.33) \begin{align} \tau _{\textit{c}}(\ell ) &=\text{``}\frac {\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3}\left [ 1-\cos \left ( 2\pi {\boldsymbol{k}} \boldsymbol{\cdot }\frac {{\boldsymbol{\ell }}}{\ell }\right )\right ]\tau _{\textit{f}}(|{\boldsymbol{k}}|/\ell ){\mathbb{E}} \big [\left |\widehat {{\boldsymbol{u}}}(t,{\boldsymbol{k}}/\ell )\right |^2\big ]{\mathrm{d}}{\boldsymbol{k}}}{\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3}\left [ 1-\cos \left ( 2\pi {\boldsymbol{k}} \boldsymbol{\cdot }\frac {{\boldsymbol{\ell }}}{\ell }\right )\right ]{\mathbb{E}}\big [ \left |\widehat {{\boldsymbol{u}}}(t,{\boldsymbol{k}}/\ell )\right |^2\big ]{\mathrm{d}}{\boldsymbol{k}}}\text{''}\notag \\[3pt] &\mathrel {\mathop {\propto }\limits _{\ell \to 0}^{}}\frac {\ell }{\sigma }\frac {\int _{\rho \geqslant 0,\theta \in [0,\pi ]}\left [ 1-\cos \left ( 2\pi \rho \cos \theta \right )\right ]\rho ^{-8/3}{\mathrm{d}}\rho {\mathrm{d}}\theta }{\int _{\rho \geqslant 0,\theta \in [0,\pi ]}\left [ 1-\cos \left ( 2\pi \rho \cos \theta \right )\right ]\rho ^{-5/3}{\mathrm{d}}\rho {\mathrm{d}}\theta }, \end{align}

where we have used the asymptotic behaviours of $\tau _{\textit{f}}(k)$ (2.32) and of the PSD ${\mathbb{E}} [\left |\widehat {{\boldsymbol{u}}}(t,{\boldsymbol{k}})\right |^2 ]\sim |{\boldsymbol{k}}|^{-11/3}$ (2.13) at large argument. Note that the remaining ratio of double integrals entering in (2.33) can be shown to be finite (an exact expression can even be obtained using a symbolic calculation software). The relevant result given in (2.33) is the fact that, as a consequence of the observation that $\tau _{\textit{f}}(k)$ becomes inversely proportional to $k$ at large $k$ (2.32), the genuine characteristic time scale $\tau _{\textit{c}}(\ell )$ of eddies of typical size $\ell$ (2.29) becomes proportional to $\ell$ (2.33). Consequently, the relevant characteristic velocity appearing in (2.33) is $\sigma$ , i.e. the characteristic velocity of the largest eddies, while it may have been more natural to obtain their typical velocity $\sqrt {S_2(\ell )}$ as in the expression (2.28) of $\tau _{\text{e}}(\ell )$ : this is the phenomenon of advection of the small scales by the large eddies, referred to here and in the turbulence literature as the sweeping effect.

2.4. Implications for the spatio-temporal covariance matrix of the velocity field

Now let us present the general form of the covariance structure in time and space of the velocity field that encompasses both the spatial structure detailed in § 2.1 and the peculiar temporal structure that is governed by the sweeping effect, as is the case in § 2.3. To do so, let us begin by defining the relevant spatio-temporal covariance of the velocity field $\mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }})$ as

(2.34) \begin{align} \mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }}) = {\mathbb{E}} \big [ u_i^{\nu }(t,{\boldsymbol{x}})u_{\!j}^{\nu }(t+\tau ,{\boldsymbol{x}}+{\boldsymbol{\ell }})\big ]. \end{align}

Again, $\mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }})$ is independent of position $\boldsymbol{x}$ and time $t$ by statistical homogeneity and stationarity, respectively. We mention that, moreover, the spatiotemporal covariance function (2.34) is an even function of the temporal argument $\tau \in {\mathbb{R}}$ , $\mathcal C_{\textit{ij}}^{\nu }(-\tau ,{\boldsymbol{\ell }}) = \mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }})$ . In the following, we also assume statistical isotropy and explore the consequences of this assumption on the structure of the covariance function. The natural extension of the spatial covariance matrix (2.6) to a spatio-temporal framework $\mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }})$ (2.34) can be conveniently written in Fourier space as

(2.35) \begin{align} \mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }}) = \frac {1}{2} \int _{{\boldsymbol{k}}\in {\mathbb{R}}^3} \textit{e}^{2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}}F(\tau /T_{\boldsymbol{k}}) E_\nu ^{\textit{E}}({\boldsymbol{k}}) \widehat {P}_{\textit{ij}}({\boldsymbol{k}} ){\mathrm{d}} {\boldsymbol{k}} , \end{align}

where $F$ encodes the temporal dependence, assumed to be an integrable function normalised such that $F(0)=1$ , $T_{\boldsymbol{k}}$ is a time scale that depends on the wave vector $\boldsymbol{k}$ , $ E_\nu ^{\textit{E}}({\boldsymbol{k}})$ is the PSD of the velocity field (2.5) studied above and $\widehat {P}_{\textit{ij}}({\boldsymbol{k}} )$ the Leray projector on divergence-free vector fields (2.7). The covariance structure (2.35) was initially proposed by Chaves et al. (Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003) when $F$ is chosen such that $F(\tau )=\exp (-|\tau |)$ . As explained below, this choice leads to a Markovian evolution for the Fourier modes, as noted earlier by Komorowski & Peszat (Reference Komorowski and Peszat2004).

As we have already mentioned, it is clear from (2.35) that the PSD $E_\nu ^{\textit{E}}(k)$ corresponds to the Fourier transform of $\mathcal C_{\textit{ii}}^{\nu }(0,\ell )$ . As a consequence, we recover from (2.35) the scaling behavious of the second-order structure functions (2.22) and (2.26). To make a connection with the discussion developed in § 2.3, let us establish the link between the characteristic time scale $\tau _{\textit{c}}(\ell )$ (2.29) and the $k$ -dependent time scale $T_k$ entering in $\mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }})$ (2.35). We have

(2.36) \begin{align} \tau _{\textit{c}}(\ell ) &= \frac {\int _{\tau =0}^{\infty }\left [\mathcal C_{\textit{ii}}^{\nu }(\tau ,0)-\mathcal C_{\textit{ii}}^{\nu }(\tau ,{\boldsymbol{\ell }})\right ]{\mathrm{d}} \tau }{\mathcal C_{\textit{ii}}^{\nu }(0,0)-\mathcal C_{\textit{ii}}^{\nu }(0,{\boldsymbol{\ell }})}\notag \\[3pt] &=\frac {\int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \left (1-\textit{e}^{2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}} \right )T_{\boldsymbol{k}} E^{\textit{E}}(|{\boldsymbol{k}}|){\mathrm{d}} {\boldsymbol{k}}}{\int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \left (1-\textit{e}^{2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}} \right ) E^{\textit{E}}(|{\boldsymbol{k}}|){\mathrm{d}} {\boldsymbol{k}}}\int _0^{\infty } F(s){\mathrm{d}} s. \end{align}

As suggested, in a formal manner, by the dependence on the wavenumber $k$ of the characteristic time scale $\tau _{\textit{f}}(k)$ (2.30), let us similarly assume, as has been done by Chaves et al. (Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003), that

(2.37) \begin{align} T_k\mathrel {\mathop {\sim }\limits _{k\to \infty }^{}}\frac {1}{D_3 k}, \end{align}

where $D_3$ is a free parameter of the description that has the dimension of a velocity. Similarly to the derivation of the relationship between $\tau _{\textit{c}}(\ell )$ and $\tau _{\textit{f}}(k)$ obtained in (2.33), using the scaling behaviour obtained in (2.37) and the behaviour of the PSD in the same limit (2.13), we get from (2.36), while performing the remaining ratio of integrals expressed in spherical variables,

(2.38) \begin{align} \tau _{\textit{c}}(\ell ) &\mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}}\frac {\ell }{D_3}\frac {2^{1/3} \sqrt {\pi } \varGamma ^2(1/6)}{15 \varGamma (11/6)}\int _0^{\infty } F(s){\mathrm{d}} s. \end{align}

The asymptotic behaviour in (2.38), assuming the asymptotic behaviour of $T_k$ (2.37), is consistent with the expected behaviour given in (2.33). This suggests that the free parameter $D_3$ can be expressed in units of the velocity standard deviation $\sigma$ (2.2).

Although physically insightful, the behaviour of the characteristic time scale $\tau _{\textit{c}}(\ell )$ with the eddy’s size $\ell$ (2.38) is not traditionally estimated in experiments and numerical simulations. Instead, the temporal second-order structure function and the time spectrum are usually studied.

Concerning the temporal second-order structure function $S_2^{\text{T}}(\tau )$ , let us define the velocity time increment, i.e. the variation of velocity in time at a fixed position $\boldsymbol{x}$ ,

(2.39) \begin{align} \delta _\tau {\boldsymbol{u}}^\nu (t,{\boldsymbol{x}})={\boldsymbol{u}}^\nu (t+\tau ,{\boldsymbol{x}})-{\boldsymbol{u}}^\nu (t,{\boldsymbol{x}}). \end{align}

Using the covariance structure of the velocity field (2.35), we obtain

(2.40) \begin{align} S_2^{\text{T}}(\tau )\equiv \lim _{\nu \to 0}{\mathbb{E}}\big [ \left |\delta _\tau {\boldsymbol{u}}^\nu \right |^2\big ]&=2\left [ \mathcal C_{\textit{ii}}(0,0) -\mathcal C_{\textit{ii}}(\tau ,0) \right ]\notag \\ &=2\int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \left [1-F(\tau /T_{\boldsymbol{k}}) \right ]E^{\textit{E}}(|{\boldsymbol{k}}|){\mathrm{d}} {\boldsymbol{k}}\notag \\ &= 2\int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \left [1-F(\tau /T_{{\boldsymbol{k}}/\tau }) \right ]E^{\textit{E}}(|{\boldsymbol{k}}|/\tau ){\mathrm{d}} {\boldsymbol{k}}/\tau ^3, \end{align}

such that, in the limit of vanishing time scales $\tau \to 0^+$ , using the asymptotic behaviours of $T_{\boldsymbol{k}}$ (2.37) and $E^{\textit{E}}(|{\boldsymbol{k}}|)$ (2.13) at a large wave-vector amplitude, we obtain

(2.41) \begin{align} S_2^{\text{T}}(\tau )&\mathrel {\mathop {\sim }\limits _{\tau \to 0^+}^{}} 4\frac {55}{9}D_2 D_3^{2/3}\tau ^{2/3} \int _{0}^{\infty } \left [1-F(\rho ) \right ]\rho ^{-5/3} {\mathrm{d}}\rho . \end{align}

The integral on the right-hand side of (2.41) is finite as long as $1-F(\rho )$ goes to 0 as $\rho ^a$ with $a\gt 2/3$ , which will be the case in all future situations. Expressing $D_2$ in units of the average dissipation $\varepsilon$ (compare, for instance, the expressions provided in (2.21) and (2.22)) and $D_3$ in units of velocity standard deviation $\sigma$ , the scaling behaviour obtained for the temporal second-order structure function (2.41) can be alternatively written as

(2.42) \begin{align} S_2^{\text{T}}(\tau )&\mathrel {\mathop {\propto }\limits _{\tau \to 0^+}^{}} \sigma ^{2/3} (\varepsilon \tau )^{2/3}, \end{align}

which coincides with the dimensional prediction of Tennekes (Reference Tennekes1975).

Similarly, let us consider the corresponding temporal (or frequency) spectrum $E_\nu ^{\text{T}}(\omega )$ defined by

(2.43) \begin{align} E_\nu ^{\text{T}}(\omega ) &= \int _{\tau \in {\mathbb{R}}}\textit{e}^{-2i\pi \omega \tau } {\mathbb{E}} \left [ {\boldsymbol{u}}^{\nu }(t,{\boldsymbol{x}})\boldsymbol{\cdot }{\boldsymbol{u}}^{\nu }(t+\tau ,{\boldsymbol{x}})\right ]{\mathrm{d}} \tau = \int _{\tau \in {\mathbb{R}}}\textit{e}^{-2i\pi \omega \tau } \mathcal C_{\textit{ii}}^{\nu }(\tau ,0){\mathrm{d}} \tau \\[-12pt]\nonumber\end{align}

(2.44) \begin{align} &= \int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \int _{\tau \in {\mathbb{R}}} \textit{e}^{-2i\pi \omega \tau } F(\tau /T_{\boldsymbol{k}}) E_\nu ^{\textit{E}}({\boldsymbol{k}}){\mathrm{d}} {\boldsymbol{k}}{\mathrm{d}} \tau \notag \\ &= \int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} T_{\boldsymbol{k}}\widehat {F}(\omega T_{\boldsymbol{k}}) E_\nu ^{\textit{E}}({\boldsymbol{k}}){\mathrm{d}} {\boldsymbol{k}}\notag \\ &= \int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} T_{\omega {\boldsymbol{k}}}\widehat {F}(\omega T_{\omega {\boldsymbol{k}}}) E_\nu ^{\textit{E}}(\omega |{\boldsymbol{k}}|) \omega ^3 {\mathrm{d}} {\boldsymbol{k}} , \end{align}

where $\widehat {F}$ is the Fourier transform of the function $F$ . Once again, using the expressions of $T_{\boldsymbol{k}}$ (2.37) and $E^{\textit{E}}(|{\boldsymbol{k}}|)$ (2.13) at a large wave-vector amplitude, we obtain

(2.45) \begin{align} \lim _{\nu \to 0}E_\nu ^{\text{T}}(\omega ) &= \int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \omega T_{\omega {\boldsymbol{k}}}\widehat {F}(\omega T_{\omega {\boldsymbol{k}}}) E^{\textit{E}}(\omega |{\boldsymbol{k}}|) \omega ^2{\mathrm{d}}{\boldsymbol{k}}\notag \\[3pt] &\mathrel {\mathop {\sim }\limits _{\omega \to \infty }^{}}2\frac {55}{9}D_2D_3^{2/3}\omega ^{-5/3}\int _{0}^{\infty } \widehat {F}\big(\rho ^{-1}\big) \rho ^{-8/3} {\mathrm{d}}\rho . \end{align}

The integral on the right-hand side of (2.45) will be shown to be finite for the classes of functions $F$ considered below. After being reformulated in units of $\varepsilon$ and $\sigma$ , similarly to what has been done in (2.42), the frequency dependence of the time spectrum (2.45) is given by

(2.46) \begin{align} \lim _{\nu \to 0}E_\nu ^{\text{T}}(\omega ) \mathrel {\mathop {\propto }\limits _{\omega \to \infty }^{}} \sigma ^{2/3} \varepsilon ^{2/3}\omega ^{-5/3}, \end{align}

which is consistent with the dimensional predictions of Tennekes (Reference Tennekes1975).

3.1. Formulation of the model with periodic boundary conditions: statistical structure of Fourier modes

For practical and numerical reasons, we formulate a version of the Gaussian model described in the former section on the three-dimensional periodic space $[-{L_{\textit{tot}}}/2\ ;{L_{\textit{tot}}}/2]^3$ of period $L_{\textit{tot}}$ , as has been done in DNS. The model is Gaussian, it is thus fully determined in the statistically homogeneous and stationary regime by its covariance function (2.34). The respective velocity field will thus be $L_{\textit{tot}}$ periodic $u^\nu _i(t,{\boldsymbol{x}})=u^\nu _i(t,{\boldsymbol{x}}+L_{\textit{tot}}\hat {e}_p)$ in any Cartesian direction $p=1$ , $2$ or $3$ . Let us mention that future simulations will be performed using $L_{\textit{tot}}= 2\pi$ in order to be consistent with the DNS. A fundamental consequence of considering a periodic setting is the possibility to expand the velocity field using a Fourier series and the respective Fourier modes $\widehat {u}^\nu _i(t,{\boldsymbol{k}}_{\boldsymbol{n}})$ , depending on ${\boldsymbol{k}}_{\boldsymbol{n}} = {\boldsymbol{n}}/L_{\textit{tot}}$ for ${\boldsymbol{n}}=(n_1,n_2,n_3)\in {\mathbb{Z}}^3$ . It can then be checked that the variance of the Fourier modes is finite. This solves the issue exhibited in the former section, where Fourier modes in the whole space need to be interpreted in a distributional sense. Since the velocity field ${\boldsymbol{u}}^\nu$ is real-valued, Fourier modes are Hermitian symmetric, i.e. $\overline {\widehat {u}^\nu _i(t,{\boldsymbol{k}}_{\boldsymbol{n}})} = \widehat {u}^\nu _i(t,-{\boldsymbol{k}}_{\boldsymbol{n}})= \widehat {u}^\nu _i(t,{\boldsymbol{k}}_{-{\boldsymbol{n}}})$ . This implies that the Fourier mode for ${\boldsymbol{n}}=(0,0,0)$ is real-valued.

The purpose of this section is the formulation in a periodic setting of the statistical quantities defined in the whole space in § 2.4. To do so, let us begin by defining the spatio-temporal covariance of the velocity field $\mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }})$ as

(3.1) \begin{align} \mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }}) = {\mathbb{E}} \big [ u_i^{\nu }(t,{\boldsymbol{x}})u_{\!j}^{\nu }(t+\tau ,{\boldsymbol{x}}+{\boldsymbol{\ell }})\big ]. \end{align}

The definition is identical to (2.34), however $\mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }})$ given in (3.1) is periodic. If it is furthermore required that $\mathcal C_{\textit{ij}}^{\nu }$ (3.1) characterises a statistically homogeneous and stationary incompressible velocity field, the right-hand side of (3.1) is thus independent of the position $\boldsymbol{x}$ and time $t$ . The covariance structure of the Fourier modes is then written as

(3.2) \begin{align} {\mathbb{E}} \big [ \widehat {u}_i^{\nu }(t,{\boldsymbol{k}}_{\boldsymbol{n}})\widehat {u}_{\!j}^{\nu }(t+\tau ,{\boldsymbol{k}}_{\boldsymbol{m}})\big ] = \frac {L^3_{\textit{tot}}}{2}F(\tau /T_{k_{\boldsymbol{n}}}) E_\nu ^{\textit{E}}(k_{\boldsymbol{n}})\widehat {P}_{\textit{ij}}({\boldsymbol{k}}_{\boldsymbol{n}})\delta _{{\boldsymbol{n}},-{\boldsymbol{m}}}^{(3)}, \end{align}

where $\widehat {P}_{\textit{ij}}({\boldsymbol{k}}_{\boldsymbol{n}})$ is the Leray projector (2.7) at the discrete wave vector ${\boldsymbol{k}}_{\boldsymbol{n}}$ (whose norm is $k_{\boldsymbol{n}}$ and also corresponds to $k_n$ in virtue of isotropy), $\delta _{{\boldsymbol{n}},{\boldsymbol{m}}}^{(3)}=\prod _i\delta _{n_i,m_i}$ the three-dimensional Kronecker delta function over the indices $\boldsymbol{n}$ and $\boldsymbol{m}$ , scalar functions $F$ and $T_{k_{\boldsymbol{n}}}$ encoding the temporal structure similarly to what we have seen in § 2.4, and the respective PSD defined as the expectation of the norm square of the Fourier mode, i.e.

(3.3) \begin{align} {\mathbb{E}} \big [\left |\widehat {{\boldsymbol{u}}}^{\nu }(t,{\boldsymbol{k}}_{\boldsymbol{n}})\right |^2\big ]= L^3_{\textit{tot}}E_\nu ^{\textit{E}}(k_{\boldsymbol{n}}). \end{align}

This allows us to write the covariance of the velocity field (3.1) as

(3.4) \begin{align} \mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }}) = \frac {1}{2L^3_{\textit{tot}}}\sum _{ {\boldsymbol{n}}\in {\mathbb{Z}}^3} \textit{e}^{2i\pi {\boldsymbol{k}}_{\boldsymbol{n}}\boldsymbol{\cdot }{\boldsymbol{\ell }}}F(\tau /T_{k_{\boldsymbol{n}}}) E_\nu ^{\textit{E}}(k_{\boldsymbol{n}})\widehat {P}_{\textit{ij}}({\boldsymbol{k}}_{\boldsymbol{n}}). \end{align}

For the sake of generality, let us consider an expression of the uni-dimensional PSD $E_\nu ^{\textit{E,long}}(k)$ (2.10), which we recall to deal asymptotically with Hölder-continuous functions of parameter $H=1/3$ , as observed in fluid turbulence, to any parameter $0\lt H\lt 1$ . For these regularities, the expected PSD would also exhibit a power-law behaviour as in (2.10), but with an exponent that would depend on $H$ , of the form

(3.5) \begin{align} E_\nu ^{\textit{E,long}}(k) =D_2 |k|_{1/L}^{-(2H+1)} \textit{e}^{-\eta _d k}, \end{align}

where the dissipative cutoff $\eta _d$ is a parameter of the formulation. Let us recall that the usual phenomenology of turbulence, as depicted in a self-consistent manner in § 2.1 is recovered only for $H=1/3$ . In the following, we consider the length scale $\eta _d$ as a free parameter. As shown in § 2.1, imposing a given PSD in the longitudinal case $E_\nu ^{\textit{E,long}}(k)$ (2.10) implies for statistical isotropic reasons a form of the PSD $E_\nu ^{\textit{E}}(k)$ in dimension $d=3$ through the relation provided in (2.9), giving a particularly simple form in the asymptotic limit of vanishing viscosity $E^{\textit{E}}(k)$ (2.13). Using any $0\lt H\lt 1$ , instead of the particular value $H=1/3$ , the one-dimensional spectrum (3.5) leads to, through the relation (2.9),

(3.6) \begin{align} \lim _{\nu \to 0}E_\nu ^{\textit{E}}(k) = \frac {1}{2\pi }(1+2H)(3+2H) D_2 k^2 |k|_{1/L}^{-(2H+5)}, \end{align}

which is consistent with (2.13) for $H=1/3$ . Also, for the sake of generality, and to make a connection with Chaves et al. (Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003), the time scale $T_k$ will be chosen to be

(3.7) \begin{align} T_k = \frac {1}{D_3 |k|_{1/L}^{2\beta }}, \end{align}

where $\beta$ is an additional new free parameter of the model, which takes the particular value $\beta =1/2$ to ensure the asymptotic behaviour observed in turbulence (2.37). The remaining free parameter $D_3$ is expected to be proportional to the velocity standard deviation only in the case $\beta =1/2$ and $H=1/3$ .

In the limit of an infinitely large periodic box $L_{\textit{tot}}\to \infty$ , we obtain the correspondence for any appropriate function $f$ ,

(3.8) \begin{align} \lim _{L_{\textit{tot}}\to \infty } \frac {1}{L^3_{\textit{tot}}}\sum _{ {\boldsymbol{n}}\in {\mathbb{Z}}^3} f({\boldsymbol{k}}_{\boldsymbol{n}})=\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3}f({\boldsymbol{k}}){\mathrm{d}} {\boldsymbol{k}} , \end{align}

such that, in particular, the velocity covariance function (3.4) tends to its continuous counterpart (2.35). As a consequence, all relevant statistical quantities such as variance predictions (2.14) and average dissipation (2.17) remain valid in the periodic setting for $H=1/3$ , and can easily be generalised for any $0\lt H\lt 1$ .

Let us also explore the statistical structure of the spatial structure functions at small scales implied by the power-law behaviour of the PSD given in (3.5). Concerning the variance of the longitudinal velocity increment (2.19), instead of (2.21), we obtain

(3.9) \begin{align} S_2^{\textit{long}}(\ell ) &= \lim _{L_{\textit{tot}}\to \infty } \lim _{\nu \to 0} {\mathbb{E}} \big [ \big (\delta _\ell ^{\textit{long}} {\boldsymbol{u}}^\nu \big )^2\big ]= 2\int _{k\in {\mathbb{R}}} \big ( 1-\textit{e}^{2i\pi k\ell }\big ) E^{\textit{E,long}}(k) {\mathrm{d}} k\notag \\[3pt] &\mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}} 2D_2(2 \pi )^{2 H} \cos (\pi H) \frac {\varGamma (2-2 H)}{H(1-2H)}\ell ^{2H}. \end{align}

Using the well-known relation $\varGamma (z+1)=z\varGamma (z)$ , one checks that when $H=1/3$ , the prediction (3.9) coincides with (2.21). Concerning the increment of a velocity component (2.23) and the second-order structure function (2.24), we derive under the assumption of statistical isotropy and incompressibility a generalisation of (2.25), i.e.

(3.10) \begin{align} S_{2,ij}({\boldsymbol{\ell }}) \mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}} S_2^{\textit{long}}(\ell )\left ( -H\frac {\ell _i\ell _{\!j}}{\ell ^2}+(1+H)\delta _{\textit{ij}}\right )\!, \end{align}

with the consequence that

(3.11) \begin{align} S_2(\ell )\equiv S_{2,ii}(\ell ) = \lim _{L_{\textit{tot}}\to \infty } \lim _{\nu \to 0} {\mathbb{E}} \big [ \left |\delta _\ell {\boldsymbol{u}}^\nu \right |^2\big ] &= 2\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3} \big ( 1-\textit{e}^{2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}}\big ) E^{\textit{E}}(|{\boldsymbol{k}}|){\mathrm{d}} {\boldsymbol{k}}\notag \\ &\mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}}(3+2H)S_2^{\textit{long}}(\ell ). \end{align}

Correspondingly, the characteristic time scale $\tau _{\textit{c}}(\ell )$ (2.29) reads, using the expression provided in (2.36) with the generalised forms of $ E^{\textit{E}}(|{\boldsymbol{k}}|)$ (3.6) and $T_{\boldsymbol{k}}$ (3.7), and for any appropriate function $F$ , in the limit of an infinitely large periodic box allowing us to use the correspondence depicted in (3.8),

(3.12) \begin{align} \lim _{L_{\textit{tot}}\to \infty }\tau _{\textit{c}}(\ell ) &=\frac {\int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \left (1-\textit{e}^{2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}} \right )T_{\boldsymbol{k}} E^{\textit{E}}(|{\boldsymbol{k}}|){\mathrm{d}}{\boldsymbol{k}}}{\int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \left (1-\textit{e}^{2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{\ell }}} \right ) E^{\textit{E}}(|{\boldsymbol{k}}|){\mathrm{d}}{\boldsymbol{k}}}\int _{0}^\infty F(s){\mathrm{d}} s. \end{align}

At small scales, the characteristic time scale $\tau _{\textit{c}}(\ell )$ (3.12) undergoes a transition depending on the values of $H$ and $\beta$ , following the behaviours

(3.13) \begin{align} \lim _{L_{\textit{tot}}\to \infty }\tau _{\textit{c}}(\ell ) \mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}}\frac {\ell ^{2\beta }}{D_3}\frac {\int _{\rho \geqslant 0} \left (1-\dfrac {\sin (2\pi \rho )}{2\pi \rho } \right )\rho ^{-2H-1-2\beta }{\mathrm{d}} \rho }{\int _{\rho \geqslant 0} \left (1-\dfrac {\sin (2\pi \rho )}{2\pi \rho } \right )\rho ^{-2H-1}{\mathrm{d}}\rho }\int _{0}^\infty F(s){\mathrm{d}} s \quad \mbox{ for } \quad H\lt 1-\beta\end{align}

and

(3.14) \begin{align} \lim _{L_{\textit{tot}}\to \infty }&\tau _{\textit{c}}(\ell ) \notag \\ &\mathrel {\mathop {\sim }\limits _{\ell \to 0}^{}}\frac {4\pi ^2}{3}\frac {\ell ^{2(1-H)}}{D_3}\frac {\int _{\rho \geqslant 0} \rho ^6\sqrt {\rho ^2+1/L^2}^{-2H-5-2\beta }{\mathrm{d}} \rho }{\int _{\rho \geqslant 0} \left (1-\dfrac {\sin (2\pi \rho )}{2\pi \rho } \right )\rho ^{-2H-1}{\mathrm{d}} \rho }\int _{0}^\infty F(s){\mathrm{d}} s \quad \mbox{ for } \quad H\gt 1-\beta . \end{align}

Expressions provided in (3.13) and (3.14) could be further written in an explicit way with the help of a symbolic calculation software.

3.2. Markovian evolution of the Fourier modes

Now let us discuss the time evolution of the Fourier modes $\widehat {u}^\nu _i(t,{\boldsymbol{k}}_{\boldsymbol{n}})$ , such that the correlation structure is given by (3.2), leading to the expected covariance structure $\mathcal C_{\textit{ij}}^{\nu }(\tau ,{\boldsymbol{\ell }})$ given in (3.4). This cannot be done in full generality for any function $F$ since we are asking for a Markovian evolution. Actually, very few choices for $F$ can be achieved for a given Markovian evolution, and we will only explore a very limited set of them. As far as we are concerned, we will present choices of $F$ that give a Markovian evolution in which additional underlying layers will be involved.

Before proceeding, we have to define a space-periodic and temporal Gaussian white noise ` ${\mathrm{d}} W_i(t,{\boldsymbol{x}})$ ’, a notation that we put in quotes because it has to be considered as a distribution. More precisely, let us focus on its Fourier modes ` $\widehat {{\mathrm{d}} W}_i(t,{\boldsymbol{k}}_{\boldsymbol{n}})$ ’. Being complex Gaussian and Hermitian, it is eventually fully determined by the covariance function

(3.15) \begin{align} {\mathbb{E}}\left [\widehat {{\mathrm{d}} W}_i(t,{\boldsymbol{k}}_{\boldsymbol{n}})\widehat {{\mathrm{d}} W}_{\!j}(t',{\boldsymbol{k}}_{\boldsymbol{m}}) \right ]=L^3_{\textit{tot}}\delta _{i,j}\delta _{{\boldsymbol{n}},-{\boldsymbol{m}}}^{(3)}\delta (t-t'){\mathrm{d}} t{\mathrm{d}} t', \end{align}

where $\delta (t)$ stands for the Dirac-delta distribution.

3.2.1. The statistical structure of the stationary solution: spatial fractional Gaussian vector field

The stochastic evolutions that are proposed below are designed such that their statistically stationary solutions are well defined, of finite variance and all coincide with a so-called fractional Gaussian vector field, that we denote as ${\boldsymbol{u}}^{\textit{fGv},\nu }({\boldsymbol{x}})$ . The corresponding Fourier modes are denoted by $\widehat {{\boldsymbol{u}}}^{\textit{fGv},\nu }({\boldsymbol{k}}_{\boldsymbol{n}})$ .

Let us first consider the purely spatial white noise ` ${\mathrm{d}} {\boldsymbol{W}}({\boldsymbol{x}})$ ’, with (Hermitian symmetric) Fourier modes $\widehat {{\mathrm{d}} {\boldsymbol{W}}}({\boldsymbol{k}}_{\boldsymbol{n}})$ covariance given by

(3.16) \begin{align} {\mathbb{E}}\left [\widehat {{\mathrm{d}} W}_i({\boldsymbol{k}}_{\boldsymbol{n}})\widehat {{\mathrm{d}} W}_{\!j}({\boldsymbol{k}}_{\boldsymbol{m}}) \right ]=L^3_{\textit{tot}}\delta _{i,j}\delta _{{\boldsymbol{n}},-{\boldsymbol{m}}}^{(3)}. \end{align}

The real random field ${\boldsymbol{u}}^{\textit{fGv},\nu }({\boldsymbol{x}})$ is fully defined by its Fourier mode, and is given, componentwise, by

(3.17) \begin{align} \widehat {u}^{\textit{fGv},\nu }_i({\boldsymbol{k}}_{\boldsymbol{n}}) = \sqrt {\frac {E_\nu ^{\textit{E}}(k_{\boldsymbol{n}})}{2}}\widehat {P}_{ip}({\boldsymbol{k}}_{\boldsymbol{n}})\widehat {{\mathrm{d}} W}_p({\boldsymbol{k}}_{\boldsymbol{n}}). \end{align}

The covariance coincides with (3.2) when $\tau =0$ , i.e.

(3.18) \begin{align} {\mathbb{E}} \big [\widehat {u}^{\textit{fGv},\nu }_i({\boldsymbol{k}}_{\boldsymbol{n}})\widehat {u}^{\textit{fGv},\nu }_{\!j}({\boldsymbol{k}}_{\boldsymbol{m}})\big ] = \frac {L^3_{\textit{tot}}}{2} E_\nu ^{\textit{E}}(k_{\boldsymbol{n}})\widehat {P}_{\textit{ij}}({\boldsymbol{k}}_{\boldsymbol{n}})\delta _{{\boldsymbol{n}},-{\boldsymbol{m}}}^{(3)}. \end{align}

3.2.2. The Markovian evolution of Chaves et al. (Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003)

The first natural evolution for Fourier modes, which coincides with the proposition of Chaves et al. (Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003), as recognised earlier by Komorowski & Peszat (Reference Komorowski and Peszat2004), is of an Ornstein–Uhlenbeck type, genuinely Markovian, and reads as the stochastic differential equation

(3.19) \begin{equation} d\widehat {u}^{(1),\nu }_i(t,{\boldsymbol{k}}) = -\frac {1}{T_{{\boldsymbol{k}}}}\widehat {u}^{(1),\nu }_i(t,{\boldsymbol{k}}){\mathrm{d}} t+\sqrt {\frac {E_\nu ^{\textit{E}}({\boldsymbol{k}})}{T_{{\boldsymbol{k}}}}}\widehat {P}_{ip}({\boldsymbol{k}})\widehat {{\mathrm{d}} W}_p(t,{\boldsymbol{k}}) \end{equation}

for a given discrete wave vector ${\boldsymbol{k}}={\boldsymbol{k}}_{\boldsymbol{n}}$ . Here $T_{\boldsymbol{k}}$ is the characteristic time of the Fourier mode $\boldsymbol{k}$ defined in (3.7), $\widehat {P}_{\textit{ij}}({\boldsymbol{k}} )$ the Leray projector ensuring the incompressibility defined in (2.7) and $E_\nu ^{\textit{E}}(k)$ the previously defined three-dimensional PSD obtained from (3.5) using the isotropic relation (2.9). The superscript ` $^{(1)}$ ’ appearing in the notation of the resulting velocity field $\widehat {u}^{(1),\nu }_i(t,{\boldsymbol{k}})$ (3.19) illustrates that it is obtained using a single evolution with respect to the white noise. The justification for this notation will become clear in the next paragraph. The random contribution to the evolution (3.19) is provided by the stochastic vector term $\widehat {{\mathrm{d}} W}_p(t,{\boldsymbol{k}})$ characterised by the covariance (3.15). Let us mention a very recent work by Letournel et al. (Reference Letournel, Morhain, Massot and Vié2025) devoted to an equivalent formulation of Chaves et al. (Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003) for a decomposition of the velocity field in a wavelet dyadic basis.

The unique statistically stationary solution of the stochastic dynamics proposed in (3.19) reads

(3.20) \begin{equation} \widehat {u}^{(1),\nu }_i(t,{\boldsymbol{k}}_{\boldsymbol{n}}) = \sqrt {\frac {E_\nu ^{\textit{E}}(k_{\boldsymbol{n}})}{T_{{\boldsymbol{k}}_{\boldsymbol{n}}}}}\widehat {P}_{ip}({\boldsymbol{k}}_{\boldsymbol{n}})\int _{s=-\infty }^t \textit{e}^{-\frac {t-s}{T_{{\boldsymbol{k}}_{\boldsymbol{n}}}}}\widehat {{\mathrm{d}} W}_p(s,{\boldsymbol{k}}_{\boldsymbol{n}}). \end{equation}

From the former expression (3.20), it is then straightforward to derive the covariance structure of the velocity Fourier modes that are Hermitian symmetric and characterised by their covariance

(3.21) \begin{align} {\mathbb{E}} \big [ \widehat {u}_i^{(1),\nu }(t,{\boldsymbol{k}}_{\boldsymbol{n}})\widehat {u}_{\!j}^{(1),\nu }(t+\tau ,{\boldsymbol{k}}_{\boldsymbol{m}})\big ] = \frac {L^3_{\textit{tot}}}{2}\textit{e}^{-|\tau |/T_{{\boldsymbol{k}}_{\boldsymbol{n}}}} E_\nu ^{\textit{E}}(k_{\boldsymbol{n}})\widehat {P}_{\textit{ij}}({\boldsymbol{k}}_{\boldsymbol{n}})\delta _{{\boldsymbol{n}},-{\boldsymbol{m}}}^{(3)}, \end{align}

showing that, when compared with the general model provided in (3.2), the temporal structure is given by the function

(3.22) \begin{align} F^{(1)}(\tau )=\textit{e}^{-|\tau |}, \end{align}

where once again the superscript ` $^{(1)}$ ’ entering in the particular notation $F^{(1)}$ will become clear in the next paragraph.

Let us now explore the consequences of the particular temporal kernel $F^{(1)}$ given in (3.22) on the regularity in time of the corresponding velocity field, and more precisely on the scaling behaviour of the respective temporal structure function $S_2^{(1),\text{T}}(\tau )$ , i.e. the variance of the velocity time increment $\delta _\tau {\boldsymbol{u}}^{(1),\nu }$ (2.39).

Recall that we found the expression (2.40), that reads, using the generalised forms of $ E^{\textit{E}}(|{\boldsymbol{k}}|)$ (3.6) and $T_{\boldsymbol{k}}$ (3.7),

(3.23) \begin{align} S_2^{(1),\text{T}}(\tau )&=\lim _{L_{\textit{tot}}\to \infty }\lim _{\nu \to 0}{\mathbb{E}}\Big [ \big |\delta _\tau {\boldsymbol{u}}^{(1),\nu }\big |^2\Big ]\notag \\ &= \frac {1}{\pi }(1+2H)(3+2H) D_2 \int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \Big [1-\textit{e}^{-D_3|{\boldsymbol{k}}|_{1/L}^{2\beta }|\tau |} \Big ]|{\boldsymbol{k}}|^2 |{\boldsymbol{k}}|_{1/L}^{-(2H+5)} {\mathrm{d}} {\boldsymbol{k}} . \end{align}

Similarly to the characteristic time scale $\tau _{\textit{c}}(\ell )$ (3.12), the temporal structure function $S_2^{(1),\text{T}}(\tau )$ (3.23) undergoes a transition at small scales depending on the values of $H$ and $\beta$ . We have, while rescaling the dummy integration variable $\boldsymbol{k}$ by the appropriate power of the time scale $\tau$ for $H\lt \beta$ ,

(3.24) \begin{align} S_2^{(1),\text{T}}(\tau )&\mathrel {\mathop {\sim }\limits _{|\tau |\to 0}^{}}4(1+2H)(3+2H) D_2 |D_3\tau |^{H/\beta }\int _{\rho =0}^{\infty } \Big [1-\textit{e}^{-\rho ^{2\beta }} \Big ]\rho ^{-(2H+1)} {\mathrm{d}} \rho . \end{align}

Note that when $H=1/3$ and $\beta =1/2$ , one retrieves (2.41) from the general formulation (3.24). For $H\gt \beta$ , we have

(3.25) \begin{align} S_2^{(1),\text{T}}(\tau )\mathrel {\mathop {\sim }\limits _{|\tau |\to 0}^{}}4(1+2H)(3+2H) D_2 |D_3\tau |\int _{\rho =0}^{\infty } \rho ^4\sqrt {\rho ^2+1/L^2}^{2\beta -(2H+5)} {\mathrm{d}} \rho . \end{align}

The scaling behaviours derived in (3.24) and (3.25) pertain to the inertial range. We would also like to explore the behaviour of the structure function in the dissipative range, first by taking the limit of vanishing scales $\tau \to 0$ and second by taking the limit of vanishing viscosity $\nu \to 0$ . When the viscosity $\nu$ is not zero, instead of (3.23), we have

(3.26) \begin{align} S_{2,\nu }^{(1),\text{T}}(\tau )&=\lim _{L_{\textit{tot}}\to \infty }{\mathbb{E}}\Big [ \big |\delta _\tau {\boldsymbol{u}}^{(1),\nu }\big |^2\Big ]= 2 \int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} \Big [1-\textit{e}^{-D_3|{\boldsymbol{k}}|_{1/L}^{2\beta }|\tau |} \Big ] E_\nu ^{\textit{E}}({\boldsymbol{k}}) {\mathrm{d}} {\boldsymbol{k}} , \end{align}

which will clearly behave proportionally to $|\tau |$ for any couple of values of $H$ and $\beta$ due to the behaviour of $1-F^{(1)}(\tau /T_{\boldsymbol{k}})=|\tau |/T_{\boldsymbol{k}} + o(\tau )$ in the vicinity of $\tau =0$ . In this regard, the spatio-temporal velocity field ${\boldsymbol{u}}^\nu (t,{\boldsymbol{x}})$ is not smooth in time and, thus, does not provide the expected temporal dissipative range scaling in $\tau ^2$ . This is related to the non-differentiability of the temporal Fourier mode correlation function at zero. Thus, we have to improve the model while proposing another temporal kernel in addition to the simple Markovian model of Chaves et al. (Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003) that has the exponential temporal correlation (3.22).

Let us close this paragraph by exploring the corresponding temporal spectrum $E_\nu ^{(1),\text{T}}(\omega )$ . As already considered in (2.43), for the resulting velocity field ${\boldsymbol{u}}^{(1),\nu }$ , we have

(3.27) \begin{align} E_\nu ^{(1),\text{T}}(\omega ) &= \lim _{L_{\textit{tot}}\to \infty }\int _{\tau \in {\mathbb{R}}}\textit{e}^{-2i\pi \omega \tau } {\mathbb{E}} \big [ {\boldsymbol{u}}^{(1),\nu }(t,{\boldsymbol{x}})\boldsymbol{\cdot }{\boldsymbol{u}}^{(1),\nu }(t+\tau ,{\boldsymbol{x}})\big ]{\mathrm{d}} \tau \notag \\ &= \int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} T_{\boldsymbol{k}}\widehat {F}^{(1)}(\omega T_{\boldsymbol{k}}) E_\nu ^{\textit{E}}({\boldsymbol{k}}) {\mathrm{d}} {\boldsymbol{k}} , \end{align}

where the one-dimensional Fourier transform $\widehat {F}^{(1)}$ of the temporal kernel $F^{(1)}$ (3.22) is given by

(3.28) \begin{align} \widehat {F}^{(1)}(\omega )= \frac {2}{1+4\pi ^2\omega ^2}. \end{align}

Similarly to the structure function $S_2^{(1),\text{T}}(\tau )$ (3.23), the temporal spectrum $E_\nu ^{(1),\text{T}}(\omega )$ (3.27) also undergoes a transition depending on the values of $H$ and $\beta$ . For $H\lt \beta$ , rescaling the dummy wavenumber $k$ by the appropriate power of the frequency $\omega$ and looking for the equivalent at high frequencies $\omega \to \infty$ , we obtain, for $H\lt \beta$ ,

(3.29) \begin{align} &E^{(1),\text{T}}(\omega ) \equiv \lim _{\nu \to 0} E_\nu ^{(1),\text{T}}(\omega )= \int _{{\boldsymbol{k}} \in {\mathbb{R}}^3} T_{\omega ^{\frac {1}{2\beta }} {\boldsymbol{k}}}\widehat {F}^{(1)}\left(\omega T_{\omega ^{\frac {1}{2\beta }} {\boldsymbol{k}}}\right) E^{\textit{E}}\Big(\omega ^{\frac {1}{2\beta }} |{\boldsymbol{k}}|\Big) \omega ^{\frac {3}{2\beta }} {\mathrm{d}} {\boldsymbol{k}} ,\notag \\ &\mathrel {\mathop {\sim }\limits _{\omega \to \infty }^{}}2(1+2H)(3+2H) D_2D_3^{H/\beta }\omega ^{-\left (1+H/\beta \right )}\int _{0}^\infty \rho ^{-\left (2\beta +2H+1 \right )}\widehat {F}^{(1)}(\rho ^{-2\beta }){\mathrm{d}} \rho \quad \text{ for } \quad H\lt \beta , \end{align}

which coincides with (2.45) when $H=1/3$ and $\beta =1/2$ . For $H\gt \beta$ , however, we obtain

(3.30) \begin{align} E^{(1),\text{T}}(\omega ) &\mathrel {\mathop {\sim }\limits _{\omega \to \infty }^{}}\frac {1}{4\pi ^3} (1+2H)(3+2H) D_2D_3\omega ^{-2}\int _{{\boldsymbol{k}}\in {\mathbb{R}}^3} |{\boldsymbol{k}}|^2|{\boldsymbol{k}}|_{1/L}^{2\beta -(2H+5)}{\mathrm{d}} {\boldsymbol{k}} \quad\text{ for } \quad H\gt \beta . \end{align}

Both these asymptotic behaviours (3.29) and (3.30) are consistent with the non-differentiability in time of the velocity field. In the first case (3.29), the regularity in time is determined by the value of $H/\beta$ , whereas in the second scenario (3.30), this temporal regularity coincides with that of a Brownian motion, independently of the parameters $H$ and $\beta$ .

3.2.3. Generalisation to a Markovian and differentiable-in-time framework

In § 3.2.2 we described the statistical properties in time of the Markovian model of Chaves et al. (Reference Chaves, Gawedzki, Horvai, Kupiainen and Vergassola2003) that can be seen as an Ornstein–Uhlenbeck type of evolution for Fourier modes (3.19). Fourier modes of the velocity field $\widehat {u}^{(1),\nu }_i(t,{\boldsymbol{k}})$ at a given wavenumber $\boldsymbol{k}$ are not differentiable in time because they solve a first-order evolution equation driven by a noise that is white-in-time. As a consequence, the temporal paths of $\widehat {u}^{(1),\nu }_i(t,{\boldsymbol{k}})$ share the same regularity as the Brownian motion and are only Hölder continuous. This non-differentiable nature of Fourier modes has a counterpart in the physical space: the corresponding velocity field $u^{(1),\nu }_i(t,{\boldsymbol{x}})$ is not differentiable in time. This is consistent with the behaviour of the temporal second-order structure functions $S_{2,\nu }^{(1),\text{T}}(\tau )$ , which behave proportionally to the time scale $\tau$ (3.26), even for a finite viscosity $\nu \gt 0$ . If the velocity field were differentiable, we would expect instead a quadratic behaviour, i.e. $S_{2,\nu }^{(1),\text{T}}(\tau )=O_\nu (\tau ^2)$ , where $O_\nu (\tau ^2)$ stands for a term of order $\tau ^2$ with a viscosity-dependent multiplicative factor. The necessity of proposing a dynamic capable of achieving this goal in a Markovian way is of tremendous importance, both from a physical and numerical perspective. To do so, we will exploit an idea first proposed by Sawford (Reference Sawford1991) and then generalised to an infinitely differentiable framework by Viggiano et al. (Reference Viggiano, Friedrich, Volk, Bourgoin, Cal and Chevillard2020) that consists of replacing the white noise entering in the dynamics (3.19) by a continuous random force, which is itself of Ornstein–Uhlenbeck type and/or itself governed by a dynamics stirred by a continuous random force, which would ensure the finiteness of the first and/or subsequent time derivatives.

We are now in a position to introduce the evolution for the Fourier mode $\widehat {u}^{(N),\nu }_i(t,{\boldsymbol{k}})$ . Given an integer $N\geqslant 2$ , we consider the following set of $N$ ebedded evolution equations:

(3.31) \begin{align} \frac {{\mathrm{d}} \widehat {u}^{\,(N),\nu }_i(t,{\boldsymbol{k}})}{{\mathrm{d}} t} &= -\frac {\sqrt {4N}}{T_{\boldsymbol{k}}}\widehat {u}^{\,(N),\nu }_i(t,{\boldsymbol{k}})+\widehat {f}^{\,(N-1)}_i(t,{\boldsymbol{k}}), \\[-12pt]\nonumber \end{align}

(3.32) \begin{align} \frac {{\mathrm{d}} \widehat {f}^{\,(N-1)}_i(t,{\boldsymbol{k}})}{{\mathrm{d}} t} &= -\frac {\sqrt {4N}}{T_{\boldsymbol{k}}}\widehat {f}^{\,(N-1)}_i(t,{\boldsymbol{k}})+\widehat {f}^{\,(N-2)}_i(t,{\boldsymbol{k}}), \\[-12pt]\nonumber\end{align}

(3.33) \begin{align} &\vdots \notag \\ \frac {{\mathrm{d}}\widehat {f}^{\,(2)}_i(t,{\boldsymbol{k}})}{{\mathrm{d}} t} &= -\frac {\sqrt {4N}}{T_{\boldsymbol{k}}}\widehat {f}^{\,(2)}_i(t,{\boldsymbol{k}})+ \widehat {f}^{\,(1)}_i(t,{\boldsymbol{k}}), \\[-12pt]\nonumber\end{align}

(3.34) \begin{align} {\mathrm{d}}\widehat {f}^{\,(1)}_i(t,{\boldsymbol{k}})&= -\frac {\sqrt {4N}}{T_{\boldsymbol{k}}}\widehat {f}^{\,(1)}_i(t,{\boldsymbol{k}}) {\mathrm{d}} t + \sqrt {q^{(N)}E_\nu ^{\textit{E}}({\boldsymbol{k}})}\widehat {P}_{ip}({\boldsymbol{k}})\widehat {{\mathrm{d}} W}_p(t,{\boldsymbol{k}}), \end{align}

where we introduce on the right-hand side of (3.34) the parameter $q^{(N)}$ defined by

(3.35) \begin{align} q^{(N)}&=\frac {1}{2}\frac {1}{ \int _{\mathbb R}\dfrac {1}{\Big [1+\pi ^2\frac {T_{\boldsymbol{k}}^{2}}{N}\omega ^2\Big ]^N} \textrm{d}{\kern1pt}\omega }\left (\frac {4N}{T_{\boldsymbol{k}}^2}\right )^N =\frac {1}{2}\frac {T_{\boldsymbol{k}}\sqrt {\pi } \varGamma (N)}{\sqrt {N} \varGamma (N-1/2)}\left (\frac {4N}{T_{\boldsymbol{k}}^2}\right )^N. \end{align}

Let us now clarify the Markovian nature of the evolution proposed in the set of equations (3.31)–(3.34). Although the very evolution of the respective Fourier mode $\widehat {u}^{(N),\nu }_i(t,{\boldsymbol{k}})$ cannot be considered as Markovian in itself, the $N$ -dimensional vector made up of that Fourier mode and the $(N-1)$ remaining underlying layers $\widehat {f}^{(m)}_i(t,{\boldsymbol{k}})$ for $N-1\geqslant m\geqslant 1$ , as it is defined in (4.10) considered in the numerical section (§ 4), can be formulated in a Markovian way. The respective matrix evolution is given in (4.13).

We show in Appendix A that, in the statistically stationary regime, we have

(3.36) \begin{align} \lim _{t\to \infty } {\mathbb{E}}\left [\widehat {u}^{(N),\nu }_i(t,{\boldsymbol{k}}_{\boldsymbol{n}})\widehat {u}^{(N),\nu }_{\!j}(t+\tau ,{\boldsymbol{k}}_{\boldsymbol{m}})\right ]= \frac {L^3_{\textit{tot}}}{2}F^{(N)}(\tau /T_{{\boldsymbol{k}}_{\boldsymbol{n}}})E_\nu ^{\textit{E}}(k_{\boldsymbol{n}})\widehat {P}_{\textit{ij}}({\boldsymbol{k}}_{\boldsymbol{n}})\delta ^{(3)}_{{\boldsymbol{n}},-{\boldsymbol{m}}}, \end{align}

where the temporal correlation function $F^{(N)}$ is defined as

(3.37) \begin{align} F^{(N)}(\tau ) = \frac {2 |\sqrt {N}\tau |^{N-1/2} K_{N-1/2}\big ( 2|\sqrt {N}\tau |\big )}{\varGamma (N-1/2)}. \end{align}

In (3.37), $K_n(x)$ is the modified Bessel function of the second kind. We recall that the expression of $F^{(N)}$ (3.37) should be considered for $N\geqslant 2$ . In particular, considering $F^{(N)}$ (3.37) for $N=1$ does not coincide with the expression $F^{(1)}$ given in (3.22). In this situation, the Fourier mode $\widehat {u}^{(N),\nu }_i(t,{\boldsymbol{k}}_{\boldsymbol{n}})$ is differentiable $(N-1)$ times, and we provide the expression of the variance of the $p$ derivatives (A16), which is finite as long as $p\leqslant N-1$ .

Surprisingly, the expression of the correlation function $F^{(N)}$ (3.37) coincides exactly with that of a Matérn Gaussian process, which have been extensively used in machine learning theory, as reviewed by Williams & Rasmussen (Reference Williams and Rasmussen2006). More recently, such Gaussian processes have been exploited in a turbulent modelling context by Lilly et al. (Reference Lilly, Sykulski, Early and Olhede2017). To our understanding, this type of covariance has been called in these formerly cited investigations for the necessity of considering differentiable Gaussian processes. Our approach furthermore shows that this covariance function of Matérn type (3.37) can be realised in a Markovian fashion, following the set of embedded Ornstein–Uhlenbeck evolutions provided in (3.31) to (3.34). In this sense, it can be viewed as a new way to interpret the nature of Matérn Gaussian processes.

Moreover, in the limit of an infinite number of embedded layers $N\to \infty$ , we have

(3.38) \begin{align} F^{(\infty )}(\tau )= \lim _{N\to \infty }F^{(N)}(\tau ) = \textit{e}^{-\tau ^2}. \end{align}

As a consequence, we obtain

(3.39) \begin{align} \lim _{N\to \infty }\lim _{t\to \infty } {\mathbb{E}}\big [\widehat {u}^{\,(N),\nu }_i(t,{\boldsymbol{k}}_{\boldsymbol{n}})\widehat {u}^{\,(N),\nu }_{\!j}(t+\tau ,{\boldsymbol{k}}_{\boldsymbol{m}})\big ]=\frac {L^3_{\textit{tot}}}{2}\textit{e}^{-(\tau /T_{{\boldsymbol{k}}_{\boldsymbol{n}}})^2}E_\nu ^{\textit{E}}(k_{\boldsymbol{n}})\widehat {P}_{\textit{ij}}({\boldsymbol{k}}_{\boldsymbol{n}})\delta ^{(3)}_{{\boldsymbol{n}},-{\boldsymbol{m}}}. \end{align}

Now let us explore the consequences of the particular temporal kernel $F^{(N)}$ given in (3.37) (respectively $F^{(\infty )}$ given in (3.38)) on the regularity in time of the corresponding velocity field $u^{(N),\nu }_i(t,{\boldsymbol{x}})$ (respectively $u^{(\infty ),\nu }_i(t,{\boldsymbol{x}})$ ). We gather in Appendix B the detailed derivation necessary to justify them.

Let us first focus on the temporal structure functions, i.e. the variance of the velocity time increment (2.39). Using the generalised forms of $ E^{\textit{E}}(k)$ (3.6) and $T_k$ (3.7), the expression (2.40) reads

(3.40) \begin{align} S_2^{(N),\textit{temp}}(\tau )&=\lim _{L_{\textit{tot}}\to \infty }\lim _{\nu \to 0}{\mathbb{E}}\Big [ \big |\delta _\tau {\boldsymbol{u}}^{(N),\nu }\big |^2\Big ], \end{align}

which includes the limiting case $N=\infty$ .

Similarly to the temporal structure function (3.23) associated with the exponential kernel (3.22), $S_2^{(N),\textit{temp}}(\tau )$ (3.40) eventually undergoes a transition at small scales, but for a different range of values of $H$ and $\beta$ . We indeed obtain, for any $N$ , including $N=\infty$ ,

(3.41) \begin{align} S_2^{(N),\textit{temp}}&(\tau )\mathrel {\mathop {\propto }\limits _{|\tau |\to 0}^{}}|D_3\tau |^{H/\beta } \quad \mbox{ for } \quad H\lt 2\beta . \end{align}

We provide in (B3) and (B4) an explicit expression of the multiplicative constant. Note that, for a given $\beta$ , the range of accessible spatial regularity $H$ is wider (i.e. $H\lt 2\beta$ ) in (3.41) than what was found (i.e. $H\lt \beta$ ) in (3.24). The behaviour of the structure function for $H\gt 2\beta$ is provided in (B5).

Importantly, for a finite viscosity, the limiting behaviour at small scales $\tau \to 0$ reads, as soon as $N\geqslant 2$ ,

(3.42) \begin{align} S_{2,\nu }^{(N),\textit{temp}}(\tau )\mathrel {\mathop {\propto }\limits _{|\tau |\to 0}^{}} \left | D_3\tau \right |^2. \end{align}

Once again, a more precise expression than (3.42) is given in (B6). Thus, it now behaves proportionally to $\tau ^2$ for any pair of values of $H$ and $\beta$ : the obtained velocity field ${\boldsymbol{u}}^{(N),\nu }(t,{\boldsymbol{x}})$ is now smooth in time as long as $\nu \gt 0$ and $N\geqslant 2$ .

Considering now the associated temporal spectrum $E_\nu ^{(N),\text{T}}(\omega )$ defined by

(3.43) \begin{align} E_\nu ^{(N),\text{T}}(\omega ) &= \lim _{L_{\textit{tot}}\to \infty }\int _{\tau \in {\mathbb{R}}}\textit{e}^{-2i\pi \omega \tau } {\mathbb{E}} \big [ {\boldsymbol{u}}^{(N),\nu }(t,{\boldsymbol{x}})\boldsymbol{\cdot }{\boldsymbol{u}}^{(N),\nu }(t+\tau ,{\boldsymbol{x}})\big ]{\mathrm{d}} \tau , \end{align}

we obtain, for $H\lt (2N-1)\beta$ ,

(3.44) \begin{align} E^{(N),\text{T}}(\omega ) &\equiv \lim _{\nu \to 0} E_\nu ^{(N),\text{T}}(\omega )\mathrel {\mathop {\propto }\limits _{\omega \to \infty }^{}}\omega ^{-\left (1+H/\beta \right )}, \end{align}

with a more precise expression given in (B10), and in (B11) for the limiting case $N\to \infty$ , The case $H\gt (2N-1)\beta$ is treated in (B12).

3.3. A second comparison concerning the temporal structure of DNS and the Gaussian model

In § 2.2 we compared the spatial structure of the model velocity field ${\boldsymbol{u}}^{\nu }(t,{\boldsymbol{x}})$ with the spatial structure of the DNS velocity field ${\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})$ . In this section we would like to compare the temporal structure of these velocity fields. Once again, the DNS data are obtained from the Johns Hopkins turbulence database, where 5028 time steps are stored among the evolution in a statistically stationary regime. We have downloaded the three velocity components of ${\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})$ in a limited region of space defined by three planes of coordinates $x=0$ , $y=0$ and $z=0$ .

We begin by estimating on the DNS field the time correlation of Fourier modes $\widehat {{\boldsymbol{u}}}^{\textit{DNS},\nu }(t,{\boldsymbol{k}}_{\boldsymbol{n}})$ for several values of the wave vector ${\boldsymbol{k}}_{\boldsymbol{n}}$ , as already done by Gorbunova et al. (Reference Gorbunova, Balarac, Canet, Eyink and Rossetto2021). To be more precise, and because we have downloaded data only in a limited region of space, we will be focusing on the correlation in time of the one-dimensional longitudinal Fourier modes. We derive in Appendix C the precise relationship between the time-correlation structure on the one-dimensional longitudinal Fourier modes,

(3.45) \begin{align} C_{\textit{long}}^{(N),\nu }(\tau ,k_n) = \frac {{\mathbb{E}} \Big [ \widehat {u}_{\textit{long}}^{(N),\nu }(t,k_n)\widehat {u}_{\textit{long}}^{(N),\nu }(t+\tau ,-k_n)\Big ]}{{\mathbb{E}} \Big [ \big |\widehat {u}_{\textit{long}}^{(N),\nu }(t,k_n)\big |^2\Big ]} \equiv F^{(N)}_{\textit{long},k_n}\left (\tau \right )\!, \end{align}

and its corresponding three-dimensional counterpart that is fully characterised by the temporal kernel $F$ entering in (3.2). By virtue of statistical homogeneity and isotropy, the one-dimensional longitudinal Fourier mode $\widehat {u}_{\textit{long}}^{(N),\nu }(t,k_n)$ can be defined as

(3.46) \begin{align} \widehat {u}_{\textit{long}}^{(N),\nu }(t,k_n)= \int _{[-L_{\textit{tot}}/2\ ;\ L_{\textit{tot}}/2]}\textit{e}^{-2i\pi k_n x}u^{(N),\nu }_1(t,x,0,0){\mathrm{d}} x. \end{align}

The solid lines in figure 2(b) are the theoretical predictions of $F^{(N)}_{\textit{long},k_n}(\tau )$ for the above mentioned $k_n$ , computed in Appendix C:

(3.47) \begin{align} F_{\textit{long}, k_n}^{(N)}\left (\tau \right ) = L_{\textit{tot}} \int \limits _0^{+\infty } \frac {\rho ^2}{k_n^2+\rho ^2}\frac {E^{\textit{E}}_\nu \left (\sqrt {k_n^2+\rho ^2}\right )}{2} F^{(N)}\left ( \frac {\tau }{T_{\sqrt {k_n^2 + \rho ^2}}}\right )\ 2\pi \rho {\mathrm{d}} \rho. \end{align}

In the case of $F^{(\infty )}$ and for $\beta = 1/2$ , it can be written as

(3.48) \begin{equation} F^{(\infty )}_{\textit{long},k_n} = L_{\textit{tot}}\, \textit{e}^{-\tau ^2 D_3^2 \left ( k_n^2 + L^{-2}\right ) } \int \limits _0^{+\infty } \frac {\rho ^2}{\rho ^2+k_n^2} \frac {E_\nu ^{\textit{E}}\left (\sqrt {\rho ^2+k_n^2}\right )}{2}\textit{e}^{-D_3^2 \tau ^2 \rho ^2} \ 2\pi \rho {\mathrm{d}} \rho. \end{equation}

To estimate $D_3$ from the DNS, we therefore assume the integral to vary slowly with $k_n$ and fit $F^{(\infty )}_{\textit{long},k_n}$ with a Gaussian function. This assumption may be verified by noting that all curves in figure 2(b) collapse while rescaling by $T_k$ . Following this method, the estimated value of $D_3$ is then $D_3 = 3.62$ .

Figure 2. Temporal correlation structure of the longitudinal velocity Fourier modes (3.46). (a) Longitudinal correlation function $C_{\textit{long}}^{\nu }(\tau ,k_n)$ (3.45) of the DNS velocity field $\widehat {u}_{\textit{long}}^{\textit{DNS},\nu }(t,k_n)$ , for $k_n L_{\textit{tot}}=42,\ 75,\ 107,\ 141,\ 173,\ 205$ . The characteristic large time scale $T_E^{\textit{DNS}} = 1.99$ is defined in the readme file of the DNS. Inset: rescaled correlation functions in semi-log scale for the same wavenumbers $k_n$ , with $D_3$ provided in the text and the theoretical Gaussian expression (Gorbunova et al. Reference Gorbunova, Balarac, Canet, Eyink and Rossetto2021) represented by a solid line. (b) A similar representation of the longitudinal correlation function $C_{\textit{long}}^{(4),\nu }(\tau ,k_n)$ (3.45) as in (a) but for the Gaussian model $\widehat {u}_{\textit{long}}^{(4),\nu }(t,k_n)$ . Solid lines represent theoretical predictions while dotted lines are numerical estimations for the same wavenumbers $k_n$ as in (a). The characteristic large time scale $T_E = 1.43$ is defined by $T_E := \sqrt {3}({L_{\textit{int}}}/{\sigma ^\nu })$ , with $L_{\textit{int}}$ being the integral time scale defined in (2.16) Inset: rescaled correlation functions in semi-log scale for the same wavenumbers $k_n$ , with the same $D_3$ as for DNS and provided in the text. (c) Numerical estimation of the characteristic time scale $T_{k_n}$ of the time correlation functions displayed in (a) and (b). We superimpose the numerical estimation of $T_{k_n}$ with a dashed black line given by the large $k$ asymptotic $(D_3 k)^{-1}$ of the theoretical expression (3.7), with relevant free parameters provided in the text. Inset: rescaled estimated $T_{k_n}$ by $(D_3 k)^{-1}$ , the asymptotical behaviour of $T_k$ at large $k$ .

We show in figure 2(a) the results of our numerical estimation of the time correlation function of the DNS $C_{\textit{long}}^{\nu }(\tau ,k_n)$ (3.45) for several wavenumber amplitudes $k_n$ and as a function of the time lag $\tau$ . At vanishing time lag $\tau =0$ , all curves coincide with unity by construction and then decrease toward 0 (they decorrelate) as $\tau$ increases. The characteristic time scale $T_{k_{\boldsymbol{n}}}^{\textit{DNS}}$ of decorrelation is thus clearly a decreasing function of the wavenumber $k_n$ . As observed by Gorbunova et al. (Reference Gorbunova, Balarac, Canet, Eyink and Rossetto2021), the very shape of the correlation function is close to a Gaussian function, and so, to estimate the characteristic time scale $T_{k_n}^{\textit{DNS}}$ , we perform a similar Gaussian fit to Gorbunova et al. (Reference Gorbunova, Balarac, Canet, Eyink and Rossetto2021). In a straightforward manner, for each accessible $1/L_{\textit{tot}}\lt k_n\lt 512/L_{\textit{tot}}$ , we find the best possible $T_{k_n}^{\textit{DNS}}$ such that the temporal correlation function $C_{\textit{long}}^{\nu }(\tau ,k_n)$ (3.45) resembles the explicit Gaussian function $\exp (-(\tau /T_{k_n}^{\textit{DNS}})^2)$ and reproduce the result of this fitting procedure in figure 2(c) using an orange line. We can indeed observe that $T_{k_n}^{\textit{DNS}}$ is a decreasing function of the wavenumber $k_n$ , and follows a power law, as has clearly been evidenced by the linear trends in this doubly logarithmic representation. We can also identify two ranges, below and above $k_nL_{\textit{tot}}\approx 10$ . This transition has already been observed in several studies (Kaneda et al. Reference Kaneda, Ishihara and Gotoh1999; Favier et al. Reference Favier, Godeferd and Cambon2010; Gorbunova et al. Reference Gorbunova, Balarac, Canet, Eyink and Rossetto2021) and remains unexplained, as far as we know. We believe that at large scales, i.e. at low wavenumbers, the flow is sensitive to boundary conditions and, thus, cannot be considered universal. At higher wavenumbers $k_n L_{\textit{tot}}$ , i.e. greater than 10, the power decrease is consistent with the expected behaviour given in (3.7) with $\beta =1/2$ (governing the power-law exponent) and $D_3 = 3.62$ . We superimpose with a dashed black line the proposed function dependence of $T_{k_n}$ (3.7) in the limit of large $k_n$ , with these aforementioned free parameters and clearly reproduce the decreasing behaviour. To check whether such a parameterisation of the characteristic time scale reproduces the behaviours of Fourier mode time correlation functions, we add in the inset of figure 2(a) a representation of $C_{\textit{long}}^{\nu }(\tau ,k_n)$ (3.45) as a function of the rescaled time variable $T_{k_n}^{\textit{DNS}}$ in semi-log scale, and observe that all correlation functions collapse in a satisfactory manner on the expected theoretical curve.

For comparison against our modelling approach, we consider the Gaussian velocity field $u^{{(4)},\nu }_i(t,{\boldsymbol{x}})$ , i.e. for $N=4$ . It is fully characterised in the statistically stationary regime by the covariance structure of its Fourier modes $\widehat {u}^{(4),\nu }_i(t,{\boldsymbol{k}}_{\boldsymbol{n}})$ , which is provided in (3.36). We have considered an evolution based on four layers $N=4$ that ensures a smooth behaviour in time for a finite viscosity, as shown in (B6). This choice is dictated by the necessity of ensuring smooth behaviour in time of the velocity field, corresponding to bounded velocity time derivatives (in this case, $\partial _t u^{(N),\nu }_i(t,{\boldsymbol{x}})$ have a finite variance for all integers $N\gt 1$ ), as is observed in DNS. As we explain in § 4, each layer needs to be stored in memory, which can become prohibitive given the size of the boxes, which corresponds to three components of size $1024^3$ collocation points for each layer. Besides considering $H=1/3$ , the free parameters $D_2$ , $L$ and $\eta _d$ entering in the functional form of the longitudinal PSD (3.5) have the same values as those chosen to describe the spatial structure of the velocity field, as given in § 2.2. The remaining free parameters $D_3$ and $\beta$ that appear in the characteristic correlation time scale $T_k$ of Fourier modes (3.7) are those that we have estimated before, corresponding to $\beta =1/2$ and $D_3 = 3.62$ . We would like to recall that the value of $D_3$ is not universal because it is related, when $\beta =1/2$ , to the standard deviation of the velocity. An efficient and exact in distribution numerical scheme is proposed in § 4 in order to give instances in space and time of the model field $u^{(N),\nu }_i(t,{\boldsymbol{x}})$ for all integers $N\gt 1$ . Overall, the simulation of the vector field over more than five thousand time steps took approximatively two days using 32 CPUs.

As we did for the DNS velocity field, we display in figure 2(b) the temporal longitudinal correlation functions $C_{\textit{long}}^{(N),\nu }(\tau ,k_n)$ (3.45) for the modelled Gaussian velocity field ${\boldsymbol{u}}^{(4),\nu }(t,{\boldsymbol{x}})$ . Although being known to be different from the DNS in a theoretical point of view and provided in (3.36) (once recasted into the longitudinal context (3.47)), the estimations of $C_{\textit{long}}^{(4),\nu }(\tau ,k_n)$ for various wavenumbers $k_n$ displays similar trends as for the DNS. In this case, the collapse is not perfect because we are considering only the longitudinal correlation functions, which necessitate a transformation given in (3.47) when going from the three-dimensional formulation to the one-dimensional setting. We display in figure 2(c) the time $T_k$ as a function of $k$ , for the model for which Fourier mode temporal correlation is prescribed using blue symbols, and the measured Fourier mode correlation time for DNS in orange. We observe a very satisfactory collapse of the estimations of $T_k$ from DNS and the model. We also superimpose with a dashed line the asymptotical limit $ ( D_3 k_n )^{-1}$ using the form provided in (3.7). In the inset, we display the rescaled correlation functions by the asymptotic behaviour of the imposed characteristic time scale $ ( D_3 k_n )^{-1}$ (3.7) at large $k_n$ .

We now present in figure 3 instances of the time evolution of both DNS and Gaussian model velocity fields, with the statistical estimation of the respective temporal spectra and structure functions. We begin in figure 3(a) by displaying the time evolution of the velocity norm $|{\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})|$ , along the line ${\boldsymbol{x}}\in ([-\pi ,\pi ],0,0)$ and across time $t\in [0,5028\Delta t]$ , where the time stepping $\Delta t$ is provided by the Hopkins database. In figure 3 we present a similar representation for the model $|{\boldsymbol{u}}^{(4),\nu }(t,{\boldsymbol{x}})|$ . The same colourbar has been used in both cases. Although the spatial structures of these two velocity fields are very similar, up to some filamentary structures that are not reproduced by the model (see the discussion in § 2.2 and figure 1), the spatio-temporal representation of the velocity norm exhibits a key difference between DNS and the model: for the DNS, we observe that large-scale structures of sizes typical to the integral length scale are swept by the flow, whereas it is not the case in the model. Although not clearly evidenced in this representation, we expect the absence of this sweeping property to be true for eddies at any scale. Thus, while the statistical correlation structures of Fourier modes are well reproduced by the model, as evidenced in figure 2, the model misses one of the important aspects of the sweeping phenomenon, which is the advection of the large scales, as observed in DNS. We invite the reader to the conclusion where this necessary modelling step is proposed as a perspective.

Figure 3. Direct comparison of the temporal structure of the DNS and modelled velocity fields, and estimation of PSDs and second-order structure functions. (a) Spatio-temporal representation of the DNS velocity field $|{\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})|$ along the spatial line ${\boldsymbol{x}}\in ([-\pi ,\pi ],0,0)$ and across time $t\in [0,5027\Delta t]$ , where the time stepping $\Delta t$ is provided by the Hopkins database. We use the same characteristic large time scale $T_E^{\textit{DNS}}$ as in figure 2 to adimensionalise time. (b) Similar representation to (a) but for the model $|{\boldsymbol{u}}^{(4),\nu }(t,{\boldsymbol{x}})|$ with a characteristic large time scale $T_E$ of the order of $T_E^{\textit{DNS}}$ . We use the same colourbar for panels (a) and (b). (c) Estimation of the temporal PSD $E_\nu ^{\text{T}}(\omega )$ (2.43) obtained as the variance of the temporal Fourier modes (see text). We use orange for DNS and blue for the model. The dashed black line corresponds to the power law given in (B10) without any fitting parameters. (d) Estimation of the second-order temporal structure function $S_2^{\textit{T},\nu }(\tau )$ (given by (2.40) in the inviscid limit); same colours as in (c). Dashed black lines correspond to two times the variance at large time lags, to (B3) in the inertial range, and to (B6) in the dissipative range.

Our statistical estimation of the temporal power spectrum density $E_\nu ^{\text{T}}(\omega )$ defined in (2.43) is shown in figure 3(c). To do so, we performed a discrete one-dimensional Fourier transform in time of the profile of the scalar product of velocity fields at a fixed position in space, as required in the definition (2.43). The obtained evolution in time is not periodic and so we introduce a spurious discontinuity in the picture. We have checked that windowing this profile, aimed at minimising the consequences of this discontinuity, does not alter the observed behaviours in the inertial range of scales, but it does slightly changing the results in the respective dissipative range (data not shown). As a result, the temporal spectra correspond to the variance of the (temporal) Fourier modes at each accessible frequency. The variances of the (temporal) Fourier modes of the velocity fields ${\boldsymbol{u}}^{\textit{DNS},\nu }(t,{\boldsymbol{x}})$ and ${\boldsymbol{u}}^{(N),\nu }(t,{\boldsymbol{x}})$ , obtained from the same statistical sample previously described, as functions of the frequency, are displayed doubly logarithmically in figure 3(c). We can see that with previously given free parameters, the temporal spectra superimpose in a very satisfactory manner in the inertial range and follow a typical $\omega ^{-5/3}$ power law decrease.

We finish with the estimation of the second-order temporal structure function $S_2^{\textit{T},\nu }(\tau ):= {\mathbb{E}} [ ( {\boldsymbol{u}}^\nu ( t+\tau ,{\boldsymbol{x}} ) - {\boldsymbol{u}}^\nu ( t,{\boldsymbol{x}} ) )^2 ]$ , as defined in (2.40) in the vanishing viscosity limit, and display our results in figure 3(d). The trends are consistent with the observations made for the temporal spectra. Both DNS and model velocity fields exhibit a power-law behaviour both in the inertial and the dissipative ranges, with respective power-law exponents indicated on the figure. To see this, we superimpose the prediction made in (B3) with a black dashed line. Note that in the inertial range of scales, predictions depend only weakly on the value of the number of layers N to define the model.

To conclude this comparison, let us mention that additional comparisons could be made. In particular, we could also display the full spatio-temporal correlation function, as done by Clark Di Leoni, Cobelli & Mininni (Reference Clark Di Leoni, Cobelli and Mininni2015), and the frequency-wave spectrum as estimated by de Kat & Ganapathisubramani (Reference de Kat and Ganapathisubramani2015), for both the DNS and the model. We are keeping this for further investigations. We would also like to mention that while the statistical structure of the DNS is fairly satisfactory (at least up to the aforementioned second-order statistical characterisation), our model does not capture in a fully appropriate way the advection of the eddies by some large-scale flow, as clear from the inspection of figures 3(a) and 3(b). To our knowledge, there is no model capable of doing this except some recent propositions by Armstrong & Vicol (Reference Armstrong and Vicol2025), where the temporal structure of the advecting velocity field, although simplistic in terms of regularity, is defined recursively through scales using an inverse Lagrangian flow map at each scale of the construction. We will return to these possible perspectives in the final § 5.

4.1. Remarks concerning the DFT

In all simulations we consider a finite number of Fourier modes, and thus, we need to truncate in an appropriate manner the Fourier series and perform an associated DFT. For the full benefit of the fast Fourier transform algorithm, we consider an even number of modes $N_x$ in each direction, which will furthermore be taken as a power of 2, i.e. $N_x=2^p$ , with $p$ being a given integer. In this spirit, over the grid ${\boldsymbol{x}}_{\boldsymbol{n}}={\boldsymbol{n}} L_{\textit{tot}}/N_x$ with $[\![\ -N_x/2 \ ;\ N_x/2-1\ ]\!]^3$ , here and further, we define the DFT as, for ${\boldsymbol{k}}_{\boldsymbol{m}}={\boldsymbol{m}}/L_{\textit{tot}}$ with ${\boldsymbol{m}}\in [\![\ -N_x/2 \ ;\ N_x/2-1\ ]\!]^3$ ,

(4.1) \begin{align} \widehat {u}^\nu _i(t,{\boldsymbol{k}}_{\boldsymbol{m}}) = \left (\frac {L_{\textit{tot}}}{N_x}\right )^3\sum _{ {\boldsymbol{n}}\in [\![\ -N_x/2 \ ;\ N_x/2-1\ ]\!]^3} \textit{e}^{-2i\pi {\boldsymbol{k}}_{\boldsymbol{m}}\boldsymbol{\cdot }{\boldsymbol{x}}_n} u^\nu _i(t,{\boldsymbol{x}}_n), \end{align}

with the corresponding inverse DFT

(4.2) \begin{align} u^\nu _i(t,{\boldsymbol{x}}_{\boldsymbol{n}}) = \frac {1}{L^3_{\textit{tot}}}\sum _{{\boldsymbol{m}}\in [\![\ -N_x/2 \ ;\ N_x/2-1\ ]\!]^3} \textit{e}^{2i\pi {\boldsymbol{k}}_{\boldsymbol{m}}\boldsymbol{\cdot }{\boldsymbol{x}}_{\boldsymbol{n}}} \widehat {u}^\nu _i(t,{\boldsymbol{k}}_{\boldsymbol{m}}). \end{align}

4.2. Numerical schemes

4.2.1. Numerical scheme for the Markovian evolution ( $N=1$ )

Consider first the Ornstein–Uhlenbeck process given by (3.19), and recall that in the statistically stationary regime its solution is given by (3.20). The numerical scheme we use is derived from the solution (3.20) by computing analytically the time evolution of the field $\widehat {u}^{(1),\nu }_i(t,{\boldsymbol{k}})$ at time $t_{m+1} := t_m + \delta t$ as the function of the solution at time $t_m$ :

(4.3) \begin{equation} \widehat {u}^{(1),\nu }_i(t_{m+1},{\boldsymbol{k}}) = \textit{e}^{-\frac {\delta t}{T_{\boldsymbol{k}}}}\widehat {u}^{(1),\nu }_i(t_m,{\boldsymbol{k}}) + \sqrt {\frac {E_\nu ^{\textit{E}}(k)}{T_{\boldsymbol{k}}}}\widehat {P}_{ip}({\boldsymbol{k}})\int _{s=t_{m}}^{t_{m+1}} \textit{e}^{-\frac {t_{m+1}-s}{T_{\boldsymbol{k}}}}\widehat {{\mathrm{d}} W}_p(s,{\boldsymbol{k}}). \end{equation}

The two terms on the right-hand side of (4.3) are independent Gaussian random variables; therefore, the values of $\widehat {u}^{(1),\nu }_i(t_{m},{\boldsymbol{k}})$ can be computed recursively if one is able to sample the sequence of independent Gaussian random variables $ (\int _{s=t_{m}}^{t_{m+1}} \textit{e}^{-({t_{m+1}-s}/{T_{\boldsymbol{k}}})}\widehat {{\mathrm{d}} W}_p(s,{\boldsymbol{k}}) )_{m\geqslant 0}$ . It suffices to compute the expected value

(4.4) \begin{align} {\mathbb{E}}\left [\int _{s=t_{m}}^{t_{m+1}} \textit{e}^{-\frac {t_{m+1}-s}{T_{\boldsymbol{k}}}}\widehat {{\mathrm{d}} W}_p(s,{\boldsymbol{k}})\right ]=0 \end{align}

and, using (3.15) and the Hermitian symmetry, the second-order moments are

(4.5) \begin{align} {\mathbb{E}}\left [\left |\int _{s=t_{m}}^{t_{m+1}} \textit{e}^{-\frac {t_{m+1}-s}{T_{\boldsymbol{k}}}}\widehat {{\mathrm{d}} W}_p(s,{\boldsymbol{k}})\right |^2\right ]=L^3_{\textit{tot}}\int _{s=t_m}^{t_{m+1}} \textit{e}^{-\frac {2(t_{m+1}-s)}{T_{\boldsymbol{k}}}} \text{d}s=\frac {L^3_{\textit{tot}}T_{\boldsymbol{k}}\left (1-\textit{e}^{-\frac {2\delta t}{T_{\boldsymbol{k}}}}\right )}{2}. \end{align}

For any $m$ , consider a Gaussian white noise vector field $\varGamma _{p,m}({\boldsymbol{x}})$ and define the associated Fourier modes $\widehat {\varGamma }_{p,m}({\boldsymbol{k}})$ on a periodic domain:

(4.6) \begin{align} \widehat {\varGamma }_{p,m}({\boldsymbol{k}})=\int _{[-L_{\textit{tot}}/2\ ;\ L_{\textit{tot}}/2]^3}\textit{e}^{-2i\pi {\boldsymbol{k}} \boldsymbol{\cdot }{\boldsymbol{x}}}\varGamma _{p,m}({\boldsymbol{x}}) {\mathrm{d}} {\boldsymbol{x}}. \end{align}

As a consequence of (3.15), one has the equality in distribution

(4.7) \begin{equation} \int _{s=t_{m}}^{t_{m+1}} \textit{e}^{-\frac {t_{m+1}-s}{T_{\boldsymbol{k}}}}\widehat {{\mathrm{d}} W}_p(s,{\boldsymbol{k}})=\sqrt {\frac {T_{\boldsymbol{k}}\left (1-\textit{e}^{-\frac {2\delta t}{T_{\boldsymbol{k}}}}\right )}{2}}\widehat {\varGamma }_{p,m+1}({\boldsymbol{k}}), \end{equation}

meaning that the random variables on the left- and right-hand sides of (4.7) are Gaussian, centred and have the same covariance structure.

This suggests that we introduce the numerical scheme defined by

(4.8) \begin{equation} \widehat {u}^{(1),\nu }_{i,m+1}({\boldsymbol{k}}) = \textit{e}^{-\frac {\delta t}{T_{\boldsymbol{k}}}}\widehat {u}^{(1),\nu }_{i,m}({\boldsymbol{k}}) + \sqrt {\frac {E_\nu ^{\textit{E}}({\boldsymbol{k}})\left (1-\textit{e}^{-\frac {2\delta t}{T_{\boldsymbol{k}}}}\right )}{2}}\widehat {P}_{ip}({\boldsymbol{k}})\widehat {\varGamma }_{p,m+1}({\boldsymbol{k}}), \end{equation}

which can be simulated in practice by sampling a new Gaussian white noise vector field and its Fourier modes at each iteration. This means that the numerical approximation defines a discrete time Markov process.

If the initial value $\widehat {u}^{(1),\nu }_{i,0}({\boldsymbol{k}})$ is equal in distribution to the initial value $\widehat {u}^{(1),\nu }_{i}(0,{\boldsymbol{k}})$ of (4.3) then, for all $m\geqslant 0$ , we have $\widehat {u}^{(1),\nu }_{i,m}({\boldsymbol{k}})=\widehat {u}^{(1),\nu }_{i}(t_m,{\boldsymbol{k}})$ in distribution. As a result, the proposed numerical scheme solves exactly in distribution the Ornstein–Uhlenbeck dynamics (3.19) (at grid times $t_m=m\delta t$ ).

A stationary version of the numerical approximation is obtained by choosing the initial value $\widehat {u}^{(1),\nu }_{i,0}({\boldsymbol{k}})$ distributed as the fractional Gaussian vector field studied in § 3.2.1, which can be sampled following (3.17) as

(4.9) \begin{equation} \widehat {u}^{(1),\nu }_{i,0}({\boldsymbol{k}})=\sqrt {\frac {E_\nu ^{\textit{E}}({\boldsymbol{k}})}{2}}\widehat {P}_{ip}({\boldsymbol{k}})\widehat {\varGamma }_{p,0}({\boldsymbol{k}}). \end{equation}

It is then straightforward to check that by construction the distribution of $\widehat {u}^{(1),\nu }_{i,m}({\boldsymbol{k}})$ is independent of $m$ .

Applying a numerical scheme that is exact in distribution for the Ornstein–Uhlenbeck dynamics naturally preserves its stationary Gaussian distribution, this is fundamental to obtain numerical simulations that exhibit the correct spatial regularity and statistical properties.

4.2.2. Numerical scheme for the differentiable-in-time process ( $N\geqslant 2$ )

The objective of this section is to propose a scheme for the process studied in § 3.2.3 for arbitrary $N\geqslant 2$ . As in the case $N=1$ , the aim is to define a sequence of random variables $ (\widehat {u}_{i,m}^{(N),\nu }({\boldsymbol{k}}) )_{m\geqslant 0}$ such that one has the equality in distribution $\widehat {u}_{i,m}^{(N),\nu }({\boldsymbol{k}})=\widehat {u}_{i}^{(N),\nu }(t_m,{\boldsymbol{k}})$ for any non-negative integer $m\geqslant 0$ . When $N\geqslant 2$ , the construction is more involved since the dynamics is Markovian only when considering the evolution of the additional underlying layers.

A Markovian formulation of the dynamics is obtained by considering the unknown

(4.10) \begin{equation} \widehat {\boldsymbol{F}}_i^{(N),\nu }(t,{\boldsymbol{k}})= \begin{pmatrix} \widehat {u}_i^{(N),\nu }(t,{\boldsymbol{k}})\\[5pt] \widehat {f}_i^{(N-1),\nu }(t,{\boldsymbol{k}})\\[5pt] \vdots \\[5pt] \widehat {f}_i^{(2),\nu }(t,{\boldsymbol{k}})\\[5pt] \widehat {f}_i^{(1),\nu }(t,{\boldsymbol{k}}) \end{pmatrix} \! . \end{equation}

Introducing the square matrix $\boldsymbol{A}^{(N)}({\boldsymbol{k}})$ defined by

(4.11) \begin{equation} \boldsymbol{A}^{(N)}(k)= -\frac {\sqrt {4N}}{T_k} \begin{pmatrix} 1 & 0 & \ldots & 0 & 0 \\[5pt] 0 & 1 & \ddots & \ddots & 0 \\[5pt] \vdots & \ddots & \ddots & \ddots & \vdots \\[5pt] 0 & \ddots & \ddots & 1 & 0 \\[5pt] 0 & 0 & \ldots & 0 & 1 \end{pmatrix} + \begin{pmatrix} 0 & 1 & \ldots & 0 & 0 \\[5pt] 0 & 0 & \ddots & \ddots & 0 \\[5pt] \vdots & \ddots & \ddots & \ddots & \vdots \\[5pt] 0 & \ddots & \ddots & 0 & 1 \\[5pt] 0 & 0 & \ldots & 0 & 0 \end{pmatrix} \end{equation}

and the column vector

(4.12) \begin{equation} \boldsymbol{e}^{(N)}= \begin{pmatrix} 0\\[5pt] \vdots \\[5pt] 0\\[5pt] 1 \end{pmatrix}\!, \end{equation}

the vector-valued unknown $\widehat {\boldsymbol{F}}_i^{(N),\nu }(t,{\boldsymbol{k}})$ can be interpreted as the solution to the multi-dimensional Ornstein–Uhlenbeck dynamics

(4.13) \begin{equation} {\mathrm{d}}\widehat {\boldsymbol{F}}_i^{(N),\nu }(t,{\boldsymbol{k}})=\boldsymbol{A}^{(N)}({\boldsymbol{k}})\widehat {\boldsymbol{F}}_i^{(N),\nu }(t,{\boldsymbol{k}}){\mathrm{d}} t+\sqrt {q^{(N)}E_\nu (k)}\widehat {P}_{ip}({\boldsymbol{k}})\boldsymbol{e}^{(N)}\widehat{\mathrm{d}W}_{p}(t,{\boldsymbol{k}}). \end{equation}

Considering the matrix exponential $\textit{e}^{t\boldsymbol{A}^{(N)}(k)}$ for all $t\geqslant 0$ , the variation of constants formula provides an expression for the solution at any time $t\geqslant 0$ :

(4.14) \begin{equation} \widehat {\boldsymbol{F}}_i^{(N),\nu }(t,{\boldsymbol{k}})=\textit{e}^{t\boldsymbol{A}^{(N)}(k)}\widehat {\boldsymbol{F}}_i^{(N),\nu }(0,{\boldsymbol{k}})+\sqrt {q^{(N)}E_\nu (k)}\widehat {P}_{ip}({\boldsymbol{k}})\!\int \limits _{s=0}^{t}\textit{e}^{(t-s)\boldsymbol{A}^{(N)}(k)}\boldsymbol{e}^{(N)}\widehat {\mathrm{d}W}_{p}(s,{\boldsymbol{k}}). \end{equation}

Using (4.11), it is straightforward to check that, for all $t\geqslant 0$ , one has

(4.15) \begin{equation} \textit{e}^{t\boldsymbol{A}^{(N)}(k)}=\textit{e}^{-\frac {\sqrt {4N} t}{T_{\boldsymbol{k}}}} \begin{pmatrix} 1 & t & \ldots & \dfrac {t^{N-2}}{(N-2)!} & \dfrac {t^{N-1}}{(N-1)!} \\[10pt] 0 & 1 & \ddots & \ddots & \dfrac {t^{N-2}}{(N-2)!} \\[5pt] \vdots & \ddots & \ddots & \ddots & \vdots \\[5pt] 0 & \ddots & \ddots & 1 & t \\[5pt] 0 & 0 & \ldots & 0 & 1 \end{pmatrix} \end{equation}

and, as a result, one has

(4.16) \begin{equation} \textit{e}^{t\boldsymbol{A}^{(N)}(k)}\boldsymbol{e}^{(N)} = \textit{e}^{-\frac {\sqrt {4N} t}{T_k}} \begin{pmatrix} \dfrac {t^{N-1}}{(N-1)!} \\[10pt] \dfrac {t^{N-2}}{(N-2)!} \\[5pt] \vdots \\[5pt] t \\[5pt] 1 \end{pmatrix}\! . \end{equation}

As in the derivation of the numerical scheme when $N=1$ , it is convenient to express the solution at time $t_{m+1}=t_m+\delta t$ depending on the solution at time $t_m$ as

(4.17) \begin{equation} \begin{aligned} \widehat {\boldsymbol{F}}_i^{(N),\nu }(t_{m+1},{\boldsymbol{k}})&=\textit{e}^{\delta t\boldsymbol{A}^{(N)}(k)}\widehat {\boldsymbol{F}}_i^{(N),\nu }(t_m,{\boldsymbol{k}})\\ &+\sqrt {q^{(N)}E_\nu (|{\boldsymbol{k}}|)}\widehat {P}_{ip}({\boldsymbol{k}})\int _{s=t_{m}}^{t_{m+1}}\textit{e}^{(t_{m+1}-s)\boldsymbol{A}^{(N)}(k)}\boldsymbol{e}^{(N)}\widehat {\mathrm{d}W}_{p}(s,{\boldsymbol{k}}). \end{aligned} \end{equation}

It is required to sample at each iteration the Gaussian random vector

(4.18) \begin{equation} \widehat {\boldsymbol{Z}}_{p,m}^{(N)}({\boldsymbol{k}}) = \begin{pmatrix} \widehat {Z}_{p,m}^{(N)}({\boldsymbol{k}})\\[5pt] \widehat {Z}_{p,m}^{(N-1)}({\boldsymbol{k}})\\[5pt] \vdots \\[5pt] \widehat {Z}_{p,m}^{(2)}({\boldsymbol{k}})\\[5pt] \widehat {Z}_{p,m}^{(1)}({\boldsymbol{k}}) \end{pmatrix} = \int _{s=t_{m}}^{t_{m+1}}\textit{e}^{(t_{m+1}-s)\boldsymbol{A}^{(N)}(k)}\boldsymbol{e}^{(N)}\widehat {\mathrm{d}W}_{p}(s,{\boldsymbol{k}}) ,\end{equation}

which has components given by

(4.19) \begin{equation} \widehat {Z}_{p,m}^{(n)}({\boldsymbol{k}})=\int _{s=t_{m}}^{t_{m+1}}\textit{e}^{-\frac {\sqrt {4N}(t_{m+1}-s)}{T_{\boldsymbol{k}}}}\frac {(t_{m+1}-s)^{n-1}}{(n-1)!}\widehat {\mathrm{d}W}_{p}(s,{\boldsymbol{k}}). \end{equation}

The components of $\widehat {Z}_{p,m}^{(n)}({\boldsymbol{k}})$ are not independent since they depend on the same white noise. In order to sample the centred Gaussian random vector $\widehat {\boldsymbol{Z}}_{p,m}^{(N)}({\boldsymbol{k}})$ , one needs to compute its covariance matrix $\boldsymbol{Q}^{(N)}({\boldsymbol{k}})$ with entries given by

(4.20) \begin{align} \big (\boldsymbol{Q}^{(N)}({\boldsymbol{k}})\big )_{n_1,n_2}&={\mathbb{E}}[\widehat {Z}_{p,m}^{(n_1)}({\boldsymbol{k}})\overline {\widehat {Z}_{p,m}^{(n_2)}({\boldsymbol{k}})}]\nonumber \\ &=L^3_{\textit{tot}}\int _{s=t_{m}}^{t_{m+1}}\textit{e}^{-2\frac {\sqrt {4N}(t_{m+1}-s)}{T_{\boldsymbol{k}}}}\frac {(t_{m+1}-s)^{n_1-1}}{(n_1-1)!} \frac {(t_{m+1}-s)^{n_2-1}}{(n_2-1)!}{\rm d}s\nonumber \\ &=L^3_{\textit{tot}}\int _{s=0}^{\delta t}\textit{e}^{-2\frac {\sqrt {4N}s}{T_{\boldsymbol{k}}}}\frac {s^{n_1-1}}{(n_1-1)!} \frac {s^{n_2-1}}{(n_2-1)!}{\rm d}s. \end{align}

In order to sample a centred $N$ -dimensional Gaussian random vector with a given covariance matrix $\boldsymbol{Q}$ , it is sufficient to identify any matrix $\boldsymbol{R}$ such that $\boldsymbol{Q}=\boldsymbol{R}{\kern-1pt}\boldsymbol{R}^T$ . Indeed, if $\varGamma$ is a Gaussian random vector with independent components distributed as standard centred normal random variables, then $\boldsymbol{R}\varGamma$ is a centred Gaussian random vector with covariance matrix $\boldsymbol{R}{\kern-1pt}\boldsymbol{R}^T=\boldsymbol{Q}$ . The identification of such matrices $\boldsymbol{R}$ can be performed in several ways. Recalling that the covariance matrix $\boldsymbol{Q}$ is non-negative symmetric, one may compute its square root (which requires us to first compute its eigenvalues and eigenvectors). In practice, when $\boldsymbol{Q}$ is positive definite, a more convenient approach is to compute the Cholesky decomposition of $\boldsymbol{Q}$ , which amounts to finding a matrix $\boldsymbol{R}$ that is lower triangular. The matrix $\boldsymbol{R}$ is unique if the diagonal entries of $\boldsymbol{Q}$ are positive. Moreover, Cholesky decomposition can be computed using a straightforward algorithm.

Let $\boldsymbol{R}^{(N)}({\boldsymbol{k}})$ denote the matrix obtained by the Cholesky decomposition of the covariance matrix $\boldsymbol{Q}^{(N)}({\boldsymbol{k}})$ whose entries are given above. The Gaussian random vector $\widehat {\boldsymbol{Z}}_{p,m}^{(N)}({\boldsymbol{k}})$ can be sampled as $\boldsymbol{R}^{(N)}({\boldsymbol{k}})\widehat {\boldsymbol{\varGamma }}_{p,m}^{(N)}$ , where $\widehat {\boldsymbol{\varGamma }}_{p,m}^{(N)}$ is a Gaussian random vector given by

(4.21) \begin{align} \widehat {\boldsymbol{\varGamma }}_{p,m}^{(N)}=\begin{pmatrix} \widehat {\varGamma }_{p,m}^{(N)}({\boldsymbol{k}})\\[5pt] \widehat {\varGamma }_{p,m}^{(N-1)}({\boldsymbol{k}})\\[5pt] \vdots \\[5pt] \widehat {\varGamma }_{p,m}^{(2)}({\boldsymbol{k}})\\[5pt] \widehat {\varGamma }_{p,m}^{(1)}({\boldsymbol{k}}) \end{pmatrix}\!, \end{align}

with entries $\widehat {\varGamma }_{p,m}^{(N)}({\boldsymbol{k}}),\ldots ,\widehat {\varGamma }_{p,m}^{(1)}({\boldsymbol{k}})$ given as Fourier modes of independent Gaussian white noise vector fields $\varGamma _{p,m}^{(N)}({\boldsymbol{x}}),\ldots ,\varGamma _{p,m}^{(1)}({\boldsymbol{x}})$ .

Finally, we obtain the numerical scheme defined by

(4.22) \begin{equation} \widehat {\boldsymbol{F}}_{i,m+1}^{(N),\nu }({\boldsymbol{k}})=\textit{e}^{\delta t\boldsymbol{A}^{(N)}(k)}\widehat {\boldsymbol{F}}_{i,m}^{(N),\nu }({\boldsymbol{k}})+\sqrt {q^{(N)}E_\nu (|{\boldsymbol{k}}|)}\widehat {P}_{ip}({\boldsymbol{k}})\boldsymbol{R}^{(N)}({\boldsymbol{k}})\widehat {\boldsymbol{\varGamma }}_{p,m}^{(N)}. \end{equation}

The numerical approximation defines a discrete time Markov process, where at each iteration $N$ new Gaussian white noise vector fields and their Fourier modes need to be computed.

If the initial value $\widehat {\boldsymbol{F}}_{i,0}^{(N),\nu }({\boldsymbol{k}})$ is equal in distribution to the initial value $\widehat {\boldsymbol{F}}_{i}^{(N),\nu }(0,{\boldsymbol{k}})$ of the continuous time dynamics then, for all $m\geqslant 0$ , we have $\widehat {\boldsymbol{F}}_{i,m}^{(N),\nu }({\boldsymbol{k}})=\widehat {\boldsymbol{F}}_{i}^{(N),\nu }(t_m,{\boldsymbol{k}})$ in distribution. As a result, the proposed numerical scheme solves exactly in distribution the Ornstein–Uhlenbeck dynamics (at grid times $t_m=m\delta t$ ).

A stationary version of the numerical approximation is obtained by choosing the initial value $\widehat {\boldsymbol{F}}_{i,0}^{(N),\nu }({\boldsymbol{k}})$ distributed according to the Gaussian stationary distribution described in § 3.2.3. If this holds then the distribution of $\widehat {\boldsymbol{F}}_{i,m}^{(N),\nu }({\boldsymbol{k}})$ is independent of $m$ .