# Relationship between unpredictability and intermittency in shell models of turbulence and experiments

Ewen Frogé, Carlos Granero-Belinchon, Stéphane G. Roux, Thierry Chonavel, Nicolas B. Garnier

► To cite this version:

Ewen Frogé, Carlos Granero-Belinchon, Stéphane G. Roux, Thierry Chonavel, Nicolas B. Garnier. Relationship between unpredictability and intermittency in shell models of turbulence and experiments. Physics of Fluids, 2025, 37 (11), pp.115159. 10.1063/5.0285208. hal-05110160v2

HAL Id: hal-05110160

https://hal.science/hal-05110160v2

Submitted on 15 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Relationship between unpredictability and intermittency in
shell models of turbulence and experiments

Ewen Frogé*¹,², Carlos Granero-Belinchón¹,², Sébastien G. Roux³, Thierry Chonavel¹,²,
and Nicolas B. Garnier³

¹Department of Mathematical and Electrical Engineering, IMT Atlantique, Lab-STICC, UMR CNRS 6285,, 655 Av. du Technopôle, Plouzané, 29280, Bretagne, France

²Odyssey, Inria/IMT Atlantique, 263 Av. Général Leclerc, Rennes, 35042, Bretagne,
France

³ENS de Lyon, CNRS, LPENSL, UMR5672, 69342, Lyon Cedex 07, France

October 15, 2025

Abstract

We study the predictability of turbulent velocity signals using probabilistic analog-forecasting.
Here, predictability is defined by the accuracy of forecasts and the associated uncertainties.
We study the Gledzer–Ohkitani–Yamada (GOY) shell model of turbulence as well as ex-
perimental measurements from a fully developed turbulent flow. In both cases, we identify
the extreme values of velocity at small scales as localized unpredictable events that lead
to a loss of predictability: worse mean predictions and increase of their uncertainties. The
GOY model, with its explicit scale separation, allows to evaluate the prediction performance
at individual scales, and so to better relate the intensity of extreme events and the loss of
forecast performance. Results show that predictability decreases systematically from large
to small scales. These findings establish a statistical connection between predictability loss
across scales and intermittency in turbulent flows.

1 Fully developed turbulence, intermittency and predictability

Fully developed turbulence is a highly nonlinear and multiscale phenomenon in which energy
injected at large scales is transferred through a cascade process to smaller scales, where it is
dissipated by viscosity. This cascade leads to the emergence of statistical scale invariance across
an intermediate range of scales η < l < L, known as the inertial range. The scales smaller than
the Kolmogorov scale η constitute the dissipative domain and the scales larger than the integral
scale L define the integral range. The current multifractal description of turbulence states that
structure functions S<sub>p</sub> behave as:

$$
S_p(l) = \langle |\delta_l v(x)|^p \rangle \sim l^{\zeta(p)} \quad \text{for } \eta < l < L. \qquad (1)
$$

where $\delta_l v(x) = v(x+l) - v(x)$ is the velocity increment with $l$ denoting the scale, $\langle \rangle$ is the spatial average on $x$, and $\zeta(p)$ is the scaling exponent function.

In Kolmogorov’s 1941 (K41) theory, these exponents are linear: $\zeta(p) = p/3$, corresponding
to a self-similar velocity field [19, 14]. However, empirical observations showed that $\zeta(p)$ is a
concave nonlinear function of $p$ [1, 17]. This is known as anomalous scaling and is associated
with intermittency, i.e. the irregular occurrence of intense, localized fluctuations of the kinetic

<footnote>*Corresponding author: ewen.froge@imt-atlantique.fr</footnote>
energy dissipation [23, 20]. These intense fluctuations lead to extreme events in the velocity field at small scales, both in the inertial and dissipative domains [9] that strongly impact high-order statistics.

The multifractal formalism describes the roughness of the velocity field by considering that at each location of the flow, the velocity increment behaves as $|\delta_l v(x)| \sim l^{h(x)}$, where $h(x)$ is the Hölder exponent characterizing the local regularity or roughness of the field [15, 24, 14]. The Hölder exponent also characterizes the order of the singularity at each spatial location. Contrary to the K41 theory which assumes a self-similar turbulent velocity field with $h = 1/3$ everywhere, the multifractal models consider that the turbulent velocity field contains a variety of local singularities of different orders. The set of points with the same $h$ is fractal, and its Hausdorff dimension is given by the singularity spectrum $D(h)$[15]. Having singularities of different orders is another signature of intermittency. The statistical properties of turbulence, including the structure functions, can then be derived from the distribution of local Hölder exponents[14]:

$$ \zeta(p) = \inf_{h} [ph + d - D(h)] \qquad (2) $$

where $d$ is the dimension. This formulation illustrates the relationship between anomalous scaling and the singularity distribution.

In this context, intermittency reflects the heterogeneity of the roughness of the flow, with regions of small $h$ that correspond to intense velocity gradients, and regions of mild gradients associated with large $h$. Each region contributes differently to the overall statistics. This heterogeneity explains the breakdown of self-similarity and the non-Gaussianity of velocity increments at small scales: as the scale decreases, the probability distribution of $\delta_l v$ exhibits heavy tails, reflecting the increasing contribution of rare, intense events. These features —anomalous scaling, heavy-tailed increment distributions, and spatial heterogeneity— are signatures of intermittency that can be described through the multifractal formalism.

Forecasting in turbulent flows, and in particular the prediction of extreme events in velocity gradients remains an active area of investigation [29, 7, 31]. Previous works investigated Lyapunov exponents to connect intermittency and predictability, both in dynamical shell models [11, 2] and in direct numerical simulations of homogeneous and isotropic turbulence [8], highlighting how intense small-scale fluctuations can influence large-scale predictability. In a recent paper, we explored the link between unpredictability and intermittency with a data-driven approach based on analog probabilistic forecasting applied to fully developed turbulent experimental measurements from the Modane wind tunnel [16] and showed that strong velocity fluctuations are linked to a loss of forecast accuracy.

In the present study, we extend this analysis to a simplified dynamical model of turbulence using the Gledzer–Ohkitani–Yamada (GOY) shell model. [32] Using well-separated scales by construction, this model provides a setting in which scale-specific predictability can be analyzed. This separation allows us to better understand the results observed in experimental data, particularly the relationship between intermittency and unpredictability.

We perform analog forecasting on (1) individual GOY shell variables, (2) a pseudo-velocity signal reconstructed from the sum of all GOY signals, and (3) the experimental velocity. In all cases, analog forecasting provides a probability density function (pdf) of possible future states conditioned on the past that we consider as the full probabilistic prediction. In this work, we focus on the study of its mean and variance and we illustrate their dependence on the magnitude of local fluctuations: localized high-amplitude increments are associated with worst predictions, indicating that extreme events cause larger forecasting errors in terms of both incorrect mean prediction and increase of variance. Moreover, when studying individual GOY shell variables, the estimated variance of the prediction, that we call here the uncertainty of the prediction, increases with the GOY index, i.e. when the scale decreases. This trend highlights that small scales are also the most unpredictable, with intermittency and nonlinear interactions being more
important. The statistical scarcity and intensity of extreme events dominate the prediction errors.

Very importantly, the variance of the prediction is estimated independently of the values we aim to forecast, and only depends on the pertinence of 1) the found ensemble of analogs and 2) the forecast model on this ensemble. So, localized peaks of this variance imply regions where the analogs are not behaving in a similar way. The simultaneity of the peaks of the forecast variance and the extreme velocity gradients suggests the existence of precursors of these extreme events.

The paper is organized as follows. Section II introduces the GOY shell model and the experimental turbulent velocity. Section III presents the analog forecasting method. Section IV reports results for GOY variables, pseudo-velocity, and experimental data. Section V summarizes the findings and outlines future directions.

# 2 Shell models and experimental turbulent velocity

## 2.1 The GOY shell Model

Shell models provide a simplified representation of turbulence inspired by the Fourier representation of Navier-Stokes equations. They are dynamical models of turbulence consisting of a finite set of coupled ordinary differential equations, each equation describing the dynamics of the turbulent velocity on a shell of wavenumbers. [12, 6, 10, 23] The shells, or scales, represent eddies of different sizes and are chosen to have logarithmically spaced wavenumbers. Coupling between neighboring shells is based on nonlinear triadic interactions. Shell models are built to reproduce important properties such as the nonlinear energy transfer between neighboring scales and intermittency, while enforcing conservation of global quantities such as energy and helicity.

The Gledzer-Ohkitani-Yamada (GOY) model [32] consists of a hierarchy of complex variables $\tilde{u}_n$, each associated with a wavenumber in the form:

$$k_n = k_0 \lambda^n, \quad n \in \{0, \dots, N-1\} \quad n \in \mathbb{N}, \qquad (3)$$

where $k_0$ is the smallest wavenumber, $\lambda$ sets the ratio between adjacent shells and $N$ is the total number of shells. A larger index $n$ corresponds to a larger wavenumber and hence a smaller scale. The time evolution of the system is governed by the following differential equations:

$$\begin{aligned} \frac{d\tilde{u}_n}{dt} = & ik_n \left(a\tilde{u}_{n+1}^*\tilde{u}_{n+2}^* + b\tilde{u}_{n-1}^*\tilde{u}_{n+1}^* + c\tilde{u}_{n-2}^*\tilde{u}_{n-1}^*\right) - \\ & -\nu k_n^2 \tilde{u}_n + f_n, \end{aligned} \qquad (4)$$

where $\tilde{u}_n \in \mathbb{C}$, the asterisk indicates the complex conjugate, and $a = 1.0$, $b = -\frac{\alpha}{\lambda}$ and $c = -\frac{(1-\alpha)}{\lambda^2}$ are constants that depend on the parameter $\alpha$ and are chosen to enforce energy and helicity conservation. These parameters together with Eq.(3) control the nonlinear energy transfer between shells. The viscosity $\nu$ dissipates energy at small scales and $f_n$ is an external forcing term applied at large scales.

The GOY model reproduces the Kolmogorov scale invariance for shells in the inertial range. So, the structure functions follow the scaling:

$$S_p(k_n) = \langle |\tilde{u}_n|^p \rangle \sim A_p k_n^{-\zeta(p)}, \qquad (5)$$

where $A_p$ is a constant and $\zeta(p) \sim p/3$ up to weak non-linear corrections. Consequently, the GOY model recovers intermittency, with extreme events in both inertial and dissipative domains.

For simulations, we used the following parameters:

These parameters were chosen to ensure a developed inertial range. Forcing at $n = 4$ injects energy at large scales, maintaining a statistically steady state through a cascade to smaller scales.
<figure><img src="image_2.png" /><figcaption>Figure 1: Normalized time series of the real part $u_i$ of the GOY complex variables a) near integral $u_6$, b) mid inertial $u_{10}$, and c) near dissipative $u_{14}$.</figcaption></figure>

We concentrate our study on the steady state to ensure statistical stationarity, so first tenth of the total simulation duration is discarded.

For simplicity, we restrict our analysis to the real parts of $\tilde{u}_n$. This choice is motivated by the statistical isotropy and homogeneity of the GOY model in wavenumber space, which ensure that the real and imaginary components of shell variables are statistically equivalent when averaged over long times. In particular, their power spectra and structure functions are indistinguishable when averaged over sufficiently long trajectories. Throughout this work, $u_n = \Re(\tilde{u}_n)$ will denote the real part of the shell variable. Figure 1 shows the time evolution of $u_6$, $u_{10}$ and $u_{14}$ which correspond to slow, intermediate and fast dynamics, respectively.

## 2.2 Construction of a multiscale field

The different shells of the GOY model are supposed to represent disjoint ensembles of wavenumbers in the Fourier space. Thanks to this disjointness and the linearity of the Fourier transform, we can define a multiscale field by summing the $u_n$:

$$v_{\text{GOY}}(t) = \sum_{n=0}^{N} u_n(t). \qquad (6)$$

This definition reconstructs a multiscale signal that integrates contributions from all scales, see Appendix A. We call this process GOY pseudo-velocity. Figure 2 a) illustrates the time evolution
<table><caption>Table 1: Simulation parameters used in the GOY model.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Number of shells</td><td>N = 22</td></tr><tr><td>Initial wavenumber</td><td>k<sub>0</sub> = 0.125</td></tr><tr><td>Scale ratio</td><td>λ = 2.0</td></tr><tr><td>Viscosity</td><td>ν = 10<sup>-7</sup></td></tr><tr><td>Nonlinear parameter</td><td>α = 0.5</td></tr><tr><td>Forcing strength</td><td>f<sub>4</sub> = 0.005, f<sub>i</sub> = 0, ∀i ≠ 4</td></tr><tr><td>Simulation time step</td><td>dt<sub>sim</sub> = 10<sup>-5</sup></td></tr><tr><td>Sampling time step</td><td>dt = 10<sup>-3</sup></td></tr><tr><td>Total simulation time</td><td>T<sub>max</sub> = 10<sup>5</sup></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Normalized time series of a) the GOY pseudo-velocity $v_{GOY}(t)$, b) the experimental velocity $v_{Exp}(t)$ from Modane wind tunnel.</figcaption></figure>

of this GOY pseudo-velocity which presents multiscale dynamics.

Figure 3 presents the power spectral density (PSD) of $u_n$ (colored lines) and the summed pseudo-velocity (gray line). Using the eddy turnover time $\tau(k) \sim (kv_k)^{-1}$ and the characteristic velocity at scale $k$, $v_k \sim k^{-1/3}$, temporal frequencies $\omega$ are mapped to an effective wavenumber via $k \sim \omega^{3/2}$, see [12]. This allows the power spectral density to be plotted on a wavenumber-like axis and compared to spatial spectra from experimental turbulence. The resulting spectrum of the pseudo-velocity constructed from all shells exhibits a scaling close to $k^{-5/3}$ over a finite range of scales, consistent with Kolmogorov's inertial range scaling. This illustrates the multiscale nature of the process defined in Eq. (6).

## 2.3 Experimental turbulent velocity

Turbulent velocity was measured via hot-wire anemometry in a turbulent flow in the Modane wind tunnel [18]. With a Taylor-scale Reynolds number $R_λ \approx 2500$, the flow exhibits fully developed turbulence. Using Taylor's frozen field hypothesis [30], temporal variations can be interpreted as spatial variations advected by the mean flow. This gives access to the spatial fluctuations of turbulent velocity from single-point Eulerian velocity measurements. The sampling frequency of the experiment is $f_s = \frac{1}{dt} = 25$ kHz and allows to probe scales smaller than the Kolmogorov scale that is of the order of $5dt$. Figure 2 b) presents the time evolution of the
<figure><img src="image_4.png" /><figcaption>Figure 3: Power spectral density (PSD) of individual $u_n$ (colored lines) and of the pseudo-velocity (gray line). The PSD is plotted against a surrogate wavenumber axis $k$, normalized by the forcing scale $k_{\text{forcing}}$. The theoretical $k^{-5/3}$ inertial range scaling is shown for reference (dashed black line).</figcaption></figure>

experimental turbulent velocity $v_{\text{exp}}$ over one integral scale $T = 2350dt$.

# 3 Forecast & Analogs

## 3.1 Probabilistic forecast and Innovation

For a stochastic process $X = \{x_t\}$ uniformly sampled at $dt$, future values $x_t$ of $X$ can be derived statistically from previous samples $x_{-\infty:t-dt}$. A probability density function of possible future states conditioned on the past, $q(x_t|x_{-\infty:t-dt})$, can be estimated. In this framework, the mean prediction and its uncertainty are given by the expected value of $x_t$ conditioned on its complete past $\hat{x}_t = \mathbb{E}[x_t|x_{-\infty:t-dt}]$ and the variance of $x_t$ conditioned on its complete past $\sigma_t^2 = \mathbb{V}[x_t|x_{-\infty:t-dt}]$. Consequently, the probabilistic prediction is the full pdf $q(x_t|x_{-\infty:t-dt})$. In this work, we will study it through its mean $\hat{x}_t$ and variance $\sigma_t^2$ and for simplicity of language we call $\hat{x}_t$ the prediction and $\sigma_t^2$ its uncertainty.

The innovation $\epsilon_t$ quantifies the unpredictable component of $x_t$. It is defined as:

$$ \epsilon_t = x_t - \mathbb{E}[x_t|x_{-\infty:t-dt}] \quad (7) $$

In practice, conditioning is restricted to a finite past window of length $p$, leading to the approximations:

$$ \mathbb{E}[x_t|x_{t-pdt:t-dt}] = \int x_t q(x_t|x_{t-pdt:t-dt}) dx_t \quad (8) $$

$$ \mathbb{V}[x_t|x_{t-pdt:t-dt}] = \int x_t^2 q(x_t|x_{t-pdt:t-dt}) dx_t - \left( \int x_t q(x_t|x_{t-pdt:t-dt}) dx_t \right)^2 \quad (9) $$

$$ \epsilon_t^{(p)} = x_t - \mathbb{E}[x_t|x_{t-pdt:t-dt}] \quad (10) $$
## 3.2 Analog-Based Estimation

For a given time $t$, we consider the sequence of size $p$ preceding $x_t$:

$$ \vec{x}_t^{(p)} = \begin{pmatrix} x_{t-dt} \\ \vdots \\ x_{t-pdt} \end{pmatrix} \qquad (11) $$

We say that $x_t$ is the successor of $\vec{x}_t^{(p)}$. Two sequences $\vec{x}_t^{(p)}$ and $\vec{x}_{t'}^{(p)}$ are analogs if they are close in terms of a given norm [22]. Based on Poincaré's recurrence theorem, we assume that in a sufficiently long time-series, close analogs have similar successors [27]. Consequently, for a given sequence $\vec{x}_t^{(p)}$, it is possible to perform a statistical prediction of $x_t$ based on regressions between its past analogs and their successors. This approach yields an estimate of the conditional probability density function $q^{(p)}(x_t|x_{t-pdt:t-dt})$ and in particular its mean and variance.

Given $\vec{x}_t^{(p)}$, the $k$ closest historical sequences $\{\vec{x}_{t_i}^{(p)}\}_{1\le i\le k}$ with $t_i < t$ serve as analogs [25, 21]. Following [16], in order to estimate $E[x_t|x_{t-pdt:t-dt}]$ and $W[x_t|x_{t-pdt:t-dt}]$, we use a linear regression approach, where each analog sequence is used to predict its successor. We denote by $\hat{x}_t^{(p)}$ the linear regression estimator of $x_t$ and $\hat{\sigma}_t^{(p)}$ its standard deviation. Given $k$ analogs, the regression minimizes $||\vec{y} - X_p\beta_t||_W^2$ over $\beta_t$, where:

* $\vec{y} = (x_{t'_1}, \dots, x_{t'_k})^T$ is the vector of analog successors,

* $X_p$ is the $k \times (p+1)$ matrix of past sequences augmented with a column of ones to include a bias term,

* $W$ is a diagonal weight matrix.

The least-squares estimator is:

$$ \hat{\beta}_t = (X_p^T W X_p)^{-1} X_p^T W \vec{y}, \qquad (12) $$

leading to:

$$ \hat{x}_t^{(p)} = \hat{\beta}_t^T \vec{x}_t^{(p)}, \quad \hat{\sigma}_t^{(p)} = \sqrt{\sum_{i=1}^k w_{ii} (x_{t'_i} - \hat{x}_{t'_i})^2} \qquad (13) $$

The weights $w_{ii}$ are defined as

$$ w_{ii} = \frac{\exp(-d_i/m)}{\sum_{j=1}^{k} \exp(-d_j/m)}, \qquad (14) $$

where $d_i = ||\vec{x}_t^{(p)} - \vec{x}_{t'}^{(p)}||$ is the distance to the $i$-th analog and $m$ is the median of all $d_i$, giving more weight to closer analogs.

The innovation estimate is given by

$$ \hat{\epsilon}_t^{(p)} = x_t - \hat{x}_t^{(p)} \qquad (15) $$

Very importantly, while the innovation depends on the specific value $x_t$ we aim to predict, this is not the case for $\hat{\sigma}_t^{(p)}$ that is estimated as the difference between the successors and the predictions made on the analogs. Consequently, $\hat{\sigma}_t^{(p)}$ yields information of the relevance of 1) the found analogs and 2) the fitted model on this ensemble of analogs.

To focus on the relative structure of fluctuations, each sequence of size $p$ is centered prior to analog search. This preprocessing step ensures that similarity between sequences is evaluated based on their variation rather than their mean value [16].
The choice of the forecasting parameters —the analog size $p = 3$, the number of analogs $k = 70$, and the size $N = 2^{22}$ of the reference analog dictionary— follows the settings established in our previous study[16]. These values were selected based on empirical performance and robustness criteria evaluated on experimental turbulence data. For GOY time series, the sampling time of the GOY model $dt = 10^{-3}$ simulation units, is used as prediction time step in the analog forecasting method. For the experimental turbulent velocity, the sampling time $dt = 0.4$ ms is used. In both cases, these times correspond to time scales in the dissipative domain.

# 4 Results

The GOY dataset consists of 22 shell variables $u_n \in \mathbb{R}$. For each one, the time-series is partitioned into twelve non-overlapping subsections with $N = 2^{22}$ samples in each one. The first subsection is used as a reference database for analog searching, and the probabilistic forecast is performed on the last ten. To ensure the independence between samples used in analog search and samples used in forecasting, the second partition is not used. The same procedure is applied to the pseudo-velocity field derived constructed following Eq.(6), with segments of $2^{22}$ samples, and to the experimental turbulent velocity with segments of $2^{21}$ samples. The statistics presented in the following are computed over the last ten partitions, while time series plots show only one integral scale.

## 4.1 Prediction and uncertainty in GOY

Figure 4 presents the empirical Mean Squared Error (MSE) $\frac{1}{N} \sum_{t=1}^{N} (\hat{\epsilon}_{t}^{(p)})^2$ of the forecast across all $u_n$. It corresponds to the time-averaged variance of the innovation signal and characterizes the overall predictability of the shell component.

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

Figure 4: Logarithm of the MSE of the forecast $\frac{1}{N} \sum_{t=1}^{N} (x_t - \hat{x}_t^{(p)})^2$ as a function of GOY shell index. Variables are normalized prior to any computation as $x_t \mapsto \frac{x_t - \mu}{\sigma}$, where $\mu$ and $\sigma$ are the mean and standard deviation, respectively. The color scale encodes the index for consistency with other figures.

For shells in the inertial and dissipative domains, the innovation variance clearly increases with the index, i.e. large scales are more predictable than small scales. The systematic reduction in predictability as one moves toward smaller scales is consistent with the picture of turbulence
<figure><img src="image_6.png" /><figcaption>Figure 5: Time series of the increments $\delta_{dt}u_i = u_i(t + dt) - u_i(t)$ (blue), innovation $\hat{\epsilon}_t^{(p)}$ (green), and prediction variance $\hat{\sigma}_t^{2(p)}$ (black) for the real parts of GOY shells 6, 10, and 14.</figcaption></figure>

as an energy cascade where predictability is progressively lost due to the increasing role of intermittency.

To better understand how predictability evolves across scales, we analyze three representative shells of the GOY model: one near the integral scale ($u_6$), one in the mid-inertial range ($u_{10}$), and one in the inertial domain close to the dissipative range ($u_{14}$). Figure 5 presents for the three chosen shell variables: 1) the time serie $u_i$, 2) the smallest increments $\delta_{dt}u_i$ which quantify the intensity of fluctuations over a time step $dt$, 3) the innovation $\hat{\epsilon}_t^{(p)}$ that measures the difference between the actual value and the prediction based on analog selection and 4) the prediction variance $\hat{\sigma}_t^{(p)}$ that estimates the uncertainty of the prediction.

Shell variable $u_6$ depicts a large scale and so varies smoothly over time. The increments $\delta_{dt}u_6$ remain small and evolve slowly, reflecting the variable smoothness. Its innovation $\hat{\epsilon}_t^{(p)}$ is several orders of magnitude smaller than the increments, indicating high predictability. The prediction variance $\hat{\sigma}_t^{2(p)}$ is of the order of $10^{-8}$ and is therefore negligible when compared to the variance of the increment. This indicates stable predictions and confirms the high predictability at this scale. The dynamics at this scale are largely dominated by the external large-scale forcing that is smooth, resulting in structures that evolve smoothly.

The shell 10 corresponds to an intermediate scale in the inertial domain, between the smooth evolution at larger scales and the strongly intermittent fluctuations at smaller scales in the inertial and dissipative ranges. $u_{10}$ is rougher than $u_6$, and its increment $\delta_{dt}u_{10}$ fluctuates more, as expected for the small scales of turbulence. The innovation is larger than in shell 6 but two orders of magnitude smaller than the amplitude of the increment. The prediction variance $\hat{\sigma}_t^{2(p)}$ is of the order of $10^{-7}$, increasing slightly compared to shell 6, but still negligible compared to the amplitude of the increment. These results suggest that while prediction is less effective than
at larger scales, it still approximates future values well.

Shell variable $u_{14}$, which characterizes the smallest scales of the inertial domain, shows sharper and stronger fluctuations compared to $u_6$ and $u_{10}$, and the increment $\delta_{dt}u_{14}$ is visually more intermittent. The innovation is now of the order of the amplitude of the increment and sometime even larger. This indicates that predicting the dynamics of $u_{14}$ is complicated. The prediction variance $\hat{\sigma}_t^2(p)$ also approaches values close to the variance of the increment, illustrating the increased uncertainty in the predictions. This loss of predictability observed for $u_{14}$ is consistent with the increasing significance of extreme events at small scales, a characteristic feature of intermittency in turbulence.

To better understand the loss of predictability, we analyze the variance of the innovation conditioned on the square of the increment:

$$ \langle \hat{\sigma}_t^2(p) | (\delta_{dt}u_i)^2 \rangle = \int P(\hat{\sigma}_t^2(p) | (\delta_{dt}u_i)^2) \hat{\sigma}_t^2(p) d\hat{\sigma}_t^2(p). \quad (16) $$

where $\hat{\sigma}_t^2(p)$ is the variance of the prediction and $P(\hat{\sigma}_t^2(p) | (\delta_{dt}u_i)^2)$ is the distribution of variances conditioned on the values of the square of the increments. This is done via a joint histogram of $\log[(\delta_{dt}u_i)^2]$ and $\log[\hat{\sigma}_t^2(p)]$, from which the conditional mean is computed for each increment bin. This quantity highlights whether forecasting errors are uniformly distributed across different fluctuation magnitudes or whether extreme events contribute disproportionately to the loss of predictability.

Figure 6 presents the conditioned variance of innovation for shells 6, 10, and 14. For $u_6$, the conditioned variance remains relatively stable across all increment magnitudes. Also, increments of largest amplitudes are of the order of $\log((\delta_{dt}u_i)^2) \sim -10$. This suggests that there are no extreme events in this shell variable impacting the prediction performance. Thus, large-scale dynamics are primarily driven by external forcing, leading to relatively predictable evolution. For $u_{10}$, increments with the largest amplitudes are of the order of $\log((\delta_{dt}u_i)^2) \sim -5$ and a slight increase in conditioned variance is observed for squared increment magnitude between $\log((\delta_{dt}u_i)^2) \sim -15$ and $\log((\delta_{dt}u_i)^2) \sim -5$. This suggests that while most prediction errors remain moderate, large increments lead to slightly larger forecasting errors, indicating the growing role of intermittency in shaping predictability. For $u_{14}$, the conditioned variance exhibits a steep increase with the squared increment amplitude, confirming that extreme events dominate the loss of predictability at small scales. The sharp rise in conditioned variance highlights how intermittency-driven extreme events cause unpredictability at small scales. This aligns with results of Figure 4: at small scales extreme events are becoming increasingly frequent and intense leading to a strong loss in predictability.

## 4.2 Statistical characterization of increments and innovations in GOY

To analyze the performance of the analog-based forecasting, we estimate the autocorrelation function $R(\tau)$ of the innovation $\hat{\epsilon}_t^{(p)}$, defined for a generic centered process $X(t)$ as $R(\tau) = \frac{\langle X(t)X(t+\tau) \rangle}{\langle X^2(t) \rangle}$, and its kurtosis $\frac{\langle X^4(t) \rangle}{\langle X^2(t) \rangle^2}$. The autocorrelation of the innovation quantifies the remaining linear dependencies while its kurtosis characterizes the existence of rare, large magnitude prediction errors. In theory, the innovation would be a white noise, but in practice, using a finite past of length $p$, some higher-order dependencies may remain. We compare the results to the autocorrelation function and the kurtosis of the smallest-scale increment $\delta_{dt}u_i$.

Figure 7 displays the autocorrelation of the increment (a) and the innovation (b) for $u_6$, $u_{10}$ and $u_{14}$. All shell variables present an autocorrelation of the increment that decreases and reaches zero at a time scale consistent with the eddy turnover time of the corresponding shell variable. Lower-index shell variables, such as $u_6$, exhibit broader ranges of positive autocorrelation, while higher indices, like $u_{14}$, lose correlation more quickly. This is consistent with faster dynamics at smaller scales. In contrast, the innovations of the shell variables show autocorrelation curves
<figure><img src="image_7.png" /><figcaption>Figure 6: Logarithm of the variance of innovation conditioned on the squared increment magnitude, $\langle \hat{\sigma}_t^{(p)} | (\delta_{dt}x)^2 \rangle$, for GOY $u_6$ (orange), $u_{10}$ (green), and $u_{14}$ (cyan), as well as for the pseudo velocity $v_{GOY}$ (grey) and the experimental turbulent velocity $v_{Exp}$ (black).</figcaption></figure>

that are almost zero when we move away from the origin. This absence of correlation structure indicates that the innovation is close to a white noise, even when the increment exhibits temporal coherence. The analog forecasting procedure thus removes the linearly correlated part of the signal.

Table 2 reports the kurtosis of the innovations and the increments of $u_6$, $u_{10}$, and $u_{14}$. While the increments already display high kurtosis, increasing with GOY index and indicating non-Gaussianity and growing intermittency, the kurtosis of the innovations is substantially larger. The difference between the kurtosis of increments and innovations is stronger at small scales: for $u_{14}$, the innovation kurtosis exceeds that of increments by two orders of magnitude. These values confirm that although temporal correlations are removed, the innovation retains the non-Gaussian features of turbulence — specifically, the intermittent, burst-like events that remain beyond the reach of prediction. The innovation thus concentrates the statistically extreme components of the signal that the analog method fails to forecast.

<table><caption>Table 2: Logarithm of the kurtosis of increments and innovations for GOY $u_6$, $u_{10}$, and $u_{14}$, as well as for the GOY pseudo-velocity and the experimental velocity.</caption><thead><tr><th>Signal</th><th>log Kurt[δ<sub>dt</sub>u<sub>i</sub>]</th><th>log Kurt[ε̂<sub>t</sub><sup>(p)</sup>]</th></tr></thead><tbody><tr><td>u<sub>6</sub></td><td>2.38 ± 0.11</td><td>3.92 ± 0.10</td></tr><tr><td>u<sub>10</sub></td><td>3.46 ± 0.07</td><td>8.37 ± 1.05</td></tr><tr><td>u<sub>14</sub></td><td>4.75 ± 0.09</td><td>8.21 ± 0.20</td></tr><tr><td>v<sub>GOY</sub></td><td>5.50 ± 0.25</td><td>9.23 ± 0.38</td></tr><tr><td>v<sub>Exp</sub></td><td>2.24 ± 0.03</td><td>2.79 ± 0.08</td></tr></tbody></table>

## 4.3 Prediction and uncertainty in GOY pseudo-velocity field and experimental velocity

In this section, we evaluate the performance of analogs forecast in more complex processes: the GOY pseudo-velocity field generated from Eq. (6) and an experimental turbulent velocity signal. Both processes exhibit multiscale dynamics, an energy cascade, and intermittency. In both cases,
<figure><img src="image_8.png" /></figure>

Figure 7: Normalized autocorrelation $R(\tau)$ of the increment (a) and the innovation (b) for $u_6$, $u_{10}$, and $u_{14}$ in function of $\log(\tau/T)$. Vertical dashed lines indicate the eddy turnover time of each shell variable.

we are able to link extreme events in velocity increments to unpredictable events.

In Figure 8, for the GOY pseudo-velocity $v_{GOY}(t)$ and for the experimental velocity $v_{Exp}(t)$, we present the smallest increment of the velocity $\delta_{dt}v(t)$, the innovation $\hat{\epsilon}_t^{(p)}$ and the uncertainty of the prediction $\hat{\sigma}_t^{2(p)}$, all plotted over time.

The increment of the GOY pseudo-velocity exhibits lower variance but higher intermittency with very strong peaks appearing at isolated times. The increment of $v_{Exp}$ is also intermittent but the ratio of the amplitudes of extreme and usual events is smaller. The same behavior is observed for the innovation. In the case of the GOY pseudo-velocity, most of the velocity values are well predicted. Most bursts in the innovation coincide with extreme events of the increment. In the case of $v_{Exp}$, average prediction is less accurate and isolated innovation peaks also correspond to unpredictable events linked to extreme events of the increment. Finally, the prediction uncertainty is also intermittent for both velocities, especially for the GOY pseudo-velocity, with strong peaks aligned with those in innovation, i.e. extreme events of the increment lead to particularly unpredictable events with a high uncertainty in the prediction.

Figure 6 shows the prediction variance conditioned on the squared velocity increment, for both GOY pseudo-velocity and experimental velocity. In both cases, the conditioned variance remains relatively stable for small-amplitude increments, but increases significantly for larger increments. This confirms that extreme events are the dominant contributors to predictability loss, as forecasting errors grow with the increment magnitude. Despite quantitative differences, we see that for both studied processes the predictability loss is driven mainly by rare, large-amplitude events of the velocity increment.

Figure 9 shows the autocorrelation function $R(\tau)$ of the increment and of the innovation for $v_{GOY}$ and $v_{Exp}$. In both cases, the autocorrelation of the increment decays with increasing $\tau$.
<figure><img src="image_9.png" /><figcaption>Figure 8: Time series of (a) the GOY pseudo-velocity $v_{GOY}(t)$ and (b) the experimental velocity $v_{Exp}(t)$, along with their respective increment (1), innovation (2) and variance of the prediction (3).</figcaption></figure>

For $\delta_{dt}v_{GOY}$, the decay is slower and correlation persists over a wider range of time scales. The autocorrelation of $\delta_{dt}v_{Exp}$ decreases more rapidly. In contrast, the autocorrelations of the innovation of both velocities remain almost flat and close to zero at all $\tau$, indicating that innovations are decorrelated. The absence of autocorrelation characterizes the innovation as a residual that does not retain linear temporal structure.

Figure 10 shows the second-order structure function $\log(S_2(\tau))$ of the velocity (red) and the cumulative sum of the innovation (orange), for both the GOY pseudo-velocity (a) and the experimental velocity (b), as a function of the scale of analysis $\log(\tau/T)$. In both cases, the structure function of the velocity follows a scaling close to $S_2(\tau) \sim \tau^{2/3}$ in a full range of scales, consistent with Kolmogorov's theory for fully developed turbulence. The second order structure function of the cumulative sum of the innovation shows power-law behaviors with different exponents. In the GOY case, the slope is close to 1 across all scales, which suggests that the innovation is decorrelated: the cumulative sum follows a random walk, as expected for a white-noise residual. In contrast, the cumulative sum of the innovation in the experimental velocity exhibits a two-regime behavior: at small scales, the slope is approximately 0.76, while at larger scales it transitions to a slope close to 1. This deviation from ideal white-noise behavior at small scales indicates that the innovation of experimental turbulence may retain weak correlations that are not removed by the analog forecasting method.

Table 2 reports the logarithm of the kurtosis of the smallest-scale increment and of the innovation for both the GOY pseudo-velocity and the experimental velocity. The increments of the GOY signal display a high logarithm of the kurtosis of 5.50, while the increments of $v_{Exp}$ have a significantly lower value around 2.24. This difference indicates that the GOY increments are more intermittent, with more pronounced extreme deviations from the mean. For both
<figure><img src="image_10.png" /><figcaption>Figure 9: Normalized autocorrelation function $R(\tau)$ of the smallest-scale increment (a) and the innovation (b), for the GOY pseudo-velocity $v_{\text{GOY}}$ and the experimental velocity $v_{\text{Exp}}$. Vertical dashed lines indicate the eddy turnover time of forced shell variable $u_4$ for GOY and integral scale for experimental velocity.</figcaption></figure>

signals, the kurtosis of the innovation is higher than that of the increment. In the GOY pseudo-velocity, this effect is particularly strong. In the experimental signal, the increase is smaller, but still present. This suggests that in both cases, the unpredictable residual is characterized by heavy-tailed statistics and is dominated by rare, large-amplitude deviations.

# 5 Conclusion

In this study, we analyzed the scale-dependent structure of predictability in turbulent signals using analog-based probabilistic forecasting. We studied three different systems: individual shells of the Gledzer–Ohkitani–Yamada (GOY) model, a GOY pseudo-velocity defined as the sum of all the GOY shell variables, and a high-Reynolds number experimental velocity signal from the Modane wind tunnel.

The results from the individual GOY shell variables reveal a clear increase in the mean prediction error as the scale decreases. In addition, localized extreme events in the increment of the shell variables correspond to high values of the innovation and the prediction uncertainty. This trend is also observed in the GOY pseudo-velocity and in the experimental velocity. This behavior shows that strong local fluctuations linked to intermittent events are harder to predict, and contribute significantly to the overall forecasting error.

To further characterize the statistical structure of the innovation, we analyzed its autocorrelation and kurtosis. For $u_6$, $u_{10}$, and $u_{14}$ as well as for the GOY pseudo-velocity and the experimental velocity, the autocorrelation of the innovation is nearly flat and close to zero, con-
<figure><img src="image_11.png" /></figure>

Figure 10: Second-order structure function $log(S_2(\tau))$ of the velocity and the cumulative sum of the innovation (orange) for the GOY pseudo-velocity $v_{GOY}$ (a) and experimental velocity $v_{Exp}$ (b) in function of the scale of analysis $log(\tau/T)$. For the GOY pseudo-velocity $T$ is the eddy turnover time of forced shell variable $u_4$, and for the experimental velocity $T$ is the integral scale of the flow.

firming that the analog-based method effectively removes linear dependencies regardless of the signal. Moreover, the innovation kurtosis is higher than that of the smallest scale increment, indicating the persistence of heavy tails in the innovation distribution. This suggests that extreme fluctuations remain largely unresolved and continue to dominate the prediction error statistics.

In the case of the GOY pseudo-velocity and the experimental velocity, we study the second order structure function of the velocity and the cumulative sum of the innovation. The second order structure function of both velocities presents the expected behavior for fully developed turbulence. In contrast, the structure function of the innovation of $v_{GOY}$ presents a slope 1 on a log-log plot. This is the expected behavior of a random walk, and hence corresponds to an innovation signal behaving as a white noise. In the case of experimental turbulence, this linear scaling appears only for scales larger than the integral scale, indicating that the analog-based forecast does not fully remove all the linear dependencies.

Our results show that unpredictability is not uniformly distributed across scales but increases systematically toward smaller scales. It is intuitive to consider that large scale dynamics may be easily predicted at small forecast times compared to small scale dynamics. Each GOY shell has a characteristic eddy turnover time, with low frequency shells presenting larger turnover times. So, we expect the low frequency shells to be more predictable at small forecast times than the high frequency shells. However, the variation of unpredictability across scales closely mirrors the amplitude and distribution of extreme events across scales, establishing a possible link between the two: the stronger the intermittent behavior at a given scale, the higher the associated forecast error. Small-scale shells, where intermittent dynamics are strongest, are also the most difficult
to forecast. The appearance of very localized peaks in both the innovation and the variance of the prediction that are concomitant to extreme values of the increments of the shell time series support this conclusion since the unpredictability of high frequency shells is not uniformly distributed in time. In this sense, unpredictability is not an external artifact, but an intrinsic feature of the cascade, tightly connected to the amplitude and distribution of extremes at each scale. This interpretation is consistent with a recent study showing that unresolved small-scale dynamics systematically increase the uncertainty of analog-based predictions [26], highlighting the fundamental limitations of forecasting in the absence of fine-scale information, and with the idea that microscopic fluctuations may govern the selection of flow trajectories in turbulence [28].

By applying a unified forecasting framework to both a synthetic model and experimental data, we provide evidence that the structure of forecast errors is defined by the physics of turbulence, and in particular its intermittent nature. This approach make it possible to assess predictability across scales, and to isolate the statistical contributions of extreme events to forecast uncertainty in turbulent systems. This uncertainty, defined in Eq.(13), depends on the pertinence of the found ensemble of analogs and of the forecast model applied on this ensemble. So, localized peaks of uncertainty imply heterogeneity of the process in terms of predictability. Concomitant extreme values of the forecast uncertainty and of the velocity increments indicate the presence of precursors of extreme events in the velocity field as well as the major unpredictability of these extreme events. These behaviors may be connected to the concept of spontaneous stochasticity, recently proposed as a mechanism for the emergence of intrinsic randomness in deterministic systems under singular limits [3, 4, 29]. In the context of turbulence, such randomness is hypothesized to arise from singularities in the velocity field [13], and in the high-Reynolds-number limit, it would imply both the non-uniqueness of solutions to the Navier-Stokes equations and the fundamental unpredictability of individual flow realizations [3]. In this context, analog-based ensembles may be a promising approach to improve our understanding of spontaneous stochasticity.

Finally, despite their similarity, see figure 2, we conclude that the experimental turbulence exhibits a significantly more complex structure than GOY by analyzing the average forecast performance on the experimental and GOY velocities. The pseudo-velocity signal proposed here was constructed as a real-valued combination of GOY shell variables. This reconstruction reproduces key features of turbulent velocity fields, mainly a second-order structure function consistent with the expected $S_2(\tau) \sim \tau^{2/3}$ scaling. However, the increment also exhibits extremely high kurtosis, pointing to an over-amplification of rare, extreme fluctuations. These results indicate that while the pseudo-velocity captures the global scaling behavior of turbulence, it may distort the statistical distribution of fluctuations, particularly through nonphysical intermittency.

In this work, we decided to study the mean and variance of the forecast pdf $q(x_t|x_{t-pdt:t-dt})$. However, from this conditional probability, other characterizations are possible, for example, predictions can be randomly sampled following $q(x_t|x_{t-pdt:t-dt})$, see [5]. Also, in this work, we studied forecast at a fixed very small time step $dt$ corresponding to a time scale in the dissipative domain. A detailed study at different forecast times $q(x_t|x_{t-ps:t-s})$ could illustrate how unpredictability depends on $s$. This is specially interesting in the case of the GOY shells since each shell presents a different characteristic eddy turnover time and so we expect the shells to become completely unpredictable at time lengths of the order of the turnover times.

# A Construction of a pseudo velocity field from shell Models

Shell models such as GOY provide a simplified representation of turbulence through a set of coupled ordinary differential equations of complex-valued variables $\tilde{u}_n(t)$, each associated with a shell of wavenumber $k_n$. Despite originally being referred to as “Fourier”, “spectral” or “velocity” components, these variables are not physical velocities nor do they represent direct Fourier components of a spatial velocity field.

The velocity-shell variables $\tilde{u}_n(t)$ can be understood as filtered turbulent velocities in the
<figure><img src="image_12.png" /><figcaption>Figure 11: Normalized structure functions of the reconstructed pseudo-velocity signal v<sub>GOY</sub> from GOY shell variables.</figcaption></figure>

Fourier domain, each containing information about a limited range of scales or shell. Importantly,
the GOY model is constructed to emulate the behavior of the spectral Navier-Stokes equation,
but is not derived from it directly. As a result, no evident inverse transformation exists to recover
a true velocity field from the *tilde* u<sub>n</sub>(t).

For the purpose of studying predictability and multiscale statistics, we constructed a surrogate
one-dimensional pseudo-velocity signal. We define v<sub>GOY</sub>(t) by taking the real part of the sum
of all GOY shell variables, see Eq.(6). This construction relies on the assumption that shell
variables have disjoint spectral support, ensuring that each scale is not redundantly counted in
the summation. Though heuristic, this approach yields a signal with multiscale characteristics
and temporal dynamics resembling those of real turbulent velocity fields.

Alternative reconstructions were tested, including phase-associating the shells as Fourier com-
ponents, i.e., assuming v<sub>GOY</sub>(t) = ℜ(∑<sub>n</sub> t̃<sub>n</sub>(t)e<sup>iω<sub>n</sub>t</sup>), with ω<sub>n</sub> = k<sub>n</sub><sup>2/3</sup>. However, this produced
strongly oscillatory signals with spectral spikes at each ω<sub>n</sub>, and did not reproduce a realistic
velocity spectrum.

The reconstructed signal v<sub>GOY</sub>, obtained by summing the real parts of the GOY variables,
exhibits statistical signatures that partly align with those expected in fully developed turbulence.
In particular, the second-order structure function S<sub>2</sub>(τ) displays a power-law scaling over a range
of time scales. As shown in Figure 11, panel (1), the slope of log S<sub>2</sub>(τ) is approximately 2/3
over intermediate scales, consistent with Kolmogorov's prediction. However, a closer inspection
reveals the presence of two distinct scaling regimes, with slopes ≈ 0.80 and ≈ 0.25, indicating a
deviation from ideal inertial range behavior.

Higher-order structure functions further confirm that the pseudo-velocity is not a perfect
turbulent surrogate. Figure 11, panels (2) and (3), show the skewness S<sub>3</sub>/S<sub>2</sub><sup>3/2</sup> and logarithm of
the flatness log(S<sub>4</sub>/3S<sub>2</sub><sup>2</sup>). While both quantities exhibit trends qualitatively similar to turbulence,
a negative skewness and a flatness larger than 3 in both the inertial and dissipative ranges,
their magnitudes and scaling behaviors differ significantly from those observed in experimental
data. This discrepancy may be explained by the multiscale nature of each individual GOY
shell variable. Although each u<sub>n</sub> is supposed to live in a single shell of wavenumbers, its time
series contains energetic contributions across a broad range of temporal scales, as evidenced by
their spectra (Figure 3). Even at a fixed shell, the u<sub>n</sub>(t) exhibits cascade-like dynamics, with
distinguishable integral, inertial, and dissipative regimes. Consequently, the pseudo-velocity
signal v<sub>GOY</sub>(t) aggregates overlapping multiscale processes. This goes against the assumption of
disjoint scale contributions in spectral space, which was used to justify the way the signal was
built, and may introduce artifacts in the resulting statistics.

sectionSentitivity analysis of the analog research parameters
<figure><img src="image_13.png" /><figcaption>Figure 12: Mean Squared Error $\frac{1}{N} \sum_{j=1}^{N} (x_j - \hat{x}_j^{(p)})^2$, where $x_j$ and $\hat{x}_j^{(p)}$ as a function of the number $k$ of neighbors used in the kNN algorithm. The performance is shown for $p = 3, 10, 20$.</figcaption></figure>

Figure 12 illustrates the global performance of the analog forecast in function of $k$ and $p$ for the GOY shells $u_6$, $u_{10}$ and $u_{14}$, the GOY pseudo-velocity $v_{GOY}$ and the experimental turbulent velocity $v_{Exp}$. The global performance is evaluated independently for each process with the mean squared error $\frac{1}{N} \sum_{j=1}^{N} (x_j - \hat{x}_j^{(p)})^2$, where $x_j$ and $\hat{x}_j^{(p)}$ are respectively the value we want to predict and the prediction, both at time $j$. The global performance is evaluated for $p = \{3, 10, 20\}$ and $k \in [10, 500]$, using a reference database of size $N = 2^{22}$ for GOY signals and $N = 2^{21}$ for the experimental velocity. Independently of the studied process the mean squared error reaches a plateau of best performance at large $k$ indicating the minimum value of analogs needed to provide forecasts. This plateau is reached very fast for $p = 3$ and slowly for larger $p$. For the GOY shell $u_6$ and the experimental velocity $v_{Exp}$, the best performance is obtained with $p = 20$, while for the other processes, which are rougher and more intermittent, it is with $p = 3$.

## References

[1] F. Anselmet, Y. Gagne, E. J. Hopfinger, and R. A. Antonia. High-order velocity structure functions in turbulent shear flows. *Journal of Fluid Mechanics*, 140:63–89, 1984.

[2] E. Aurell, G. Boffetta, A. Crisanti, G. Paladin, and A. Vulpiani. Predictability in systems with many characteristic times: The case of turbulence. *Physical Review E*, 53(3):2337–2349, March 1996. doi: 10.1103/PhysRevE.53.2337.

[3] Dmytro Bandak, Alexei A. Mailybaev, Gregory L. Eyink, and Nigel Goldenfeld. Spontaneous Stochasticity Amplifies Even Thermal Noise to the Largest Scales of Turbulence in a Few Eddy Turnover Times. *Physical Review Letters*, 132(10), March 2024.
[4] A. Barlet, A. Cheminet, B. Dubrulle, and A. A. Mailybaev. Spontaneous stochasticity in a 3d Weierstrass-ABC flow. *ArXiv preprint*, 2025.

[5] L. Biferale, A. A. Mailybaev, and G. Parisi. Optimal subgrid scheme for shell models of turbulence. *Physical Review E*, 95(4):043108, 2017.

[6] Luca Biferale. Shell model of energy cascade in turbulence. *Annual Review of Fluid Mechanics*, 35(Volume 35, 2003):441–468, January 2003. ISSN 0066-4189, 1545-4479. doi: 10.1146/annurev.fluid.35.101101.161122.

[7] Patrick J. Blonigan, Mohammad Farazmand, and Themistoklis P. Sapsis. Are extreme dissipation events predictable in turbulent fluid flows? *Physical Review Fluids*, 4(4):044606, April 2019. doi: 10.1103/PhysRevFluids.4.044606.

[8] G. Boffetta and S. Musacchio. Chaos and Predictability of Homogeneous-Isotropic Turbulence. *Physical Review Letters*, 119(5):054102, August 2017. doi: 10.1103/PhysRevLett.119.054102.

[9] L. Chevillard, B. Castaing, A. Arneodo, E. Lévêque, J. F. Pinton, and S. G. Roux. A phenomenological theory of Eulerian and Lagrangian velocity fluctuations in turbulent flows. *Comptes Rendus Physique*, 13(9):899–928, 2012. ISSN 1631-0705.

[10] Peter Constantin, Boris Levant, and Edriss S. Titi. Analytic Study of Shell Models of Turbulence. *Physica D: Nonlinear Phenomena*, 219(2):120–141, July 2006. ISSN 01672789. doi: 10.1016/j.physd.2006.05.015.

[11] A. Crisanti, M. H. Jensen, A. Vulpiani, and G. Paladin. Intermittency and predictability in turbulence. *Physical Review Letters*, 70(2):166–169, January 1993. doi: 10.1103/PhysRevLett.70.166.

[12] P. D. Divletsen. *Turbulence and Shell Models*. Cambridge University Press, 2010.

[13] T. D. Drivas and A. A. Mailybaev. 'Life after death' in ordinary differential equations with a non-Lipschitz singularity. *Nonlinearity*, 34:2296, 2021.

[14] U. Frisch. *Turbulence: The Legacy of A.N. Kolmogorov*. Cambridge University Press, 1995.

[15] U. Frisch and G. Parisi. On the singularity structure of fully developed turbulence. *Turbulence and Predictability in Geophysical Fluid Dynamics and Climate Dynamics*, 01:71–88, 1985.

[16] E. Frogé, C. Granero-Belinchón, S. G. Roux, N. B. Garnier, and T. Chonavel. Analog-based forecasting of turbulent velocity: Relationship between unpredictability and intermittency. *Europhysics Letters*, 149(1):13001, January 2025. ISSN 0295-5075. doi: 10.1209/0295-5075/ada4d3.

[17] Yves Gagne, Emil J. Hopfinger, and Uriel Frisch. *A New Universal Scaling for Fully Developed Turbulence: The Distribution of Velocity Increments*, pages 315–319. Springer US, New York, NY, 1990.

[18] H. Kahalerras, Y. Malecot, Y. Gagne, and B. Castaing. Intermittency and Reynolds number. *Physics of Fluids*, 10:910–921, 1998.

[19] A. Kolmogorov. The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers. *Akademiia Nauk SSSR Doklady*, 30:301–305, January 1941. ISSN 0002-3264.
[20] A. N. Kolmogorov. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. *Journal of Fluid Mechanics*, 13:82-85, 1962.

[21] Redouane Lguensat, Pierre Tandeo, Pierre Ailliot, Manuel Pulido, and Ronan Fablet. The Analog Data Assimilation. *Monthly Weather Review*, 145:4093-4107, October 2017. doi: 10.1175/MWR-D-16-0441.1.

[22] Edward N. Lorenz. Atmospheric Predictability as Revealed by Naturally Occurring Analogues. *Journal of the Atmospheric Sciences*, 26(4):636-646, July 1969. ISSN 0022-4928, 1520-0469. doi: 10.1175/1520-0469(1969)26<636:APARBN>2.0.CO;2.

[23] A. M. Obukhov. Some specific features of atmospheric turbulence. *Journal of Fluid Mechanics*, 13:77-81, 1962.

[24] Giovanni Paladin and Angelo Vulpiani. Anomalous scaling laws in multifractal objects. *Physics Reports*, 156(4):147-225, December 1987. ISSN 0370-1573. doi: 10.1016/0370-1573(87)90110-4.

[25] P. Platzer, P. Yiou, P. Naveau, P. Tandeo, Y. Zhen, P. Ailliot, and J.-F. Filipot. Using local dynamics to explain analog forecasting of chaotic systems. *Journal of the Atmospheric Sciences*, 78(7):2117-2133, July 2021. ISSN 0022-4928, 1520-0469. doi: 10.1175/JAS-D-20-0204.1.

[26] Paul Platzer and Bertrand Chapron. The Effects of Unresolved Scales on Analogue Forecasting Ensembles. pages 223-241. August 2024. ISBN 978-3-031-70659-2. doi: 10.1007/978-3-031-70660-8_10.

[27] Henri Poincaré. Sur le problème des trois corps et les équations de la dynamique. *Acta mathematica*, 13(1):A3-A270, 1890.

[28] David Ruelle. Microscopic fluctuations and turbulence. *Physics Letters A*, 72(2):81-82, June 1979. ISSN 0375-9601. doi: 10.1016/0375-9601(79)90653-4.

[29] W. Ruffenach, E. Simonnet, and N. Valade. Spontaneous stochasticity and the Armstrong-Vicol passive scalar. *ArXiv preprint*, 2025.

[30] G. I. Taylor. The Spectrum of Turbulence. *Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 164(919):476-490, 1938. ISSN 0080-4630.

[31] Alberto Vela-Martín. Predictability of isotropic turbulence by massive ensemble forecasting. *Physical Review Fluids*, 9(12):L122601, December 2024. doi: 10.1103/PhysRevFluids.9.L122601.

[32] M. Yamada and K. Ohkitani. Lyapunov spectrum of a chaotic model of three-dimensional turbulence. *Journal of the Physical Society of Japan*, 56(12):4210-4213, 1987.