1. Introduction
Over 140 years after Reynolds’ experiment (Reynolds Reference Reynolds1883), the transition to turbulence of Newtonian fluids flowing in a rigid smooth pipe continues to puzzle researchers, (e.g. the recent review by Avila, Barkley & Hof (Reference Avila, Barkley and Hof2023)). The problem depends on the Reynolds number
$\textit{Re}=({\textit{UD}}/{\nu})$
, where
$U$
is the bulk velocity,
$D$
the pipe diameter and
$\nu$
the kinematic viscosity of the fluid, and on the amplitude and shape of the disturbances to the laminar flow. Specifically, pipe flow is linearly stable for at least
$\textit{Re} \leqslant 10^{7}$
(Meseguer & Trefethen Reference Meseguer and Trefethen2003) but still, when sufficiently perturbed, it transitions to turbulence at
$\textit{Re} \approx 2000$
. It first does so in the form of localised turbulent patches known as turbulent puffs (Wygnanski & Champagne Reference Wygnanski and Champagne1973). We show a snapshot of a turbulent puff in figure 1(a). Once triggered, the puff dynamics become independent of the disturbance that created it. Turbulent puffs can either elongate and split, thereby increasing the turbulent fraction, or suddenly decay, resulting in laminar flow, figure 1(c). Decay events are more frequent than split events at
$\textit{Re} \lesssim 2040$
(Avila et al. Reference Avila, Moxey, De Lozar, Avila, Barkley and Hof2011), and obey a memoryless process: the statistics of decay events follow an exponential distribution (Faisst & Eckhardt Reference Faisst and Eckhardt2004) with a mean lifetime
$\tau$
that scales super-exponentially with the Reynolds number (Hof et al. Reference Hof, De, Alberto, Jan and Westerweel2008; Avila, Willis & Hof Reference Avila, Willis and Hof2010). Lifetime statistics at
$\textit{Re}=1850$
are shown in figure 1(d). Goldenfeld, Guttenberg & Gioia (Reference Goldenfeld, Guttenberg and Gioia2010) theorised that the super-exponential scaling of the mean lifetime of decay events can be explained through extreme value theory. Nemoto & Alexakis (Reference Nemoto and Alexakis2021) derived the super-exponential scaling by assuming that relaminarisation is linked to extremely low axial vorticity events throughout the structure of the puff. However, their study raises the following questions: Are extremely low values of the axial vorticity a cause or a consequence of decay events? Is their proposed threshold for decay actually a sufficient, or even a necessary, condition for decay?
Figure 1. Decaying events in transitional shear flows. (a) Snapshot of a turbulent puff at
$\textit{Re}=1850$
. Grey denotes low axial velocity streaks
$u_{x}' \approx -0.4$
. Red denotes regions where
$u_{r}^{2}+u_{\theta }^{2} \geqslant 0.02$
. (b) Cross-sectionally averaged, cross-flow kinetic energy
$q= \langle u_{r}^{2}+u_{\theta }^{2} \rangle _{r,\theta }$
, of the puff in (a). (c) Volume-averaged cross-sectional kinetic energy, (2.2), of a decaying puff. The marker denotes the instant of time shown in the snapshot and
$t_{d}$
stands for the time at which we detect decay, (2.4). (d) Survival probability of the MFE model (
$\textit{Re}=400$
) and puffs (
$\textit{Re}=1850$
). Time is normalised by the corresponding Lyapunov time. The dotted line corresponds to the value of the exponential distribution proposed by Avila et al. (Reference Avila, Moxey, De Lozar, Avila, Barkley and Hof2011) that fits the experimental data of Hof et al. (Reference Hof, De, Alberto, Jan and Westerweel2008). Here,
$S \approx \exp (-t/\tau )$
, where
$\tau = \exp ( \exp ( 5.56\times 10^{-3} \boldsymbol{\cdot }\textit{Re}-8.5) )$
. (e) Decaying trajectory of the MFE model (Moehlis et al. Reference Moehlis, Faisst and Eckhardt2004). Lighter colour means higher
$j$
.
Memoryless decay events are not limited to pipe flow, and are also found in other transitional shear flows such as Taylor–Couette (Borrero-Echeverry, Schatz & Tagg Reference Borrero-Echeverry, Schatz and Tagg2010) and Couette (Bottin & Chaté Reference Bottin and Chaté1998; Shi, Avila & Hof Reference Shi, Avila and Hof2013) flows. In channel flows, turbulent bands also decay following a memoryless process in short spanwise (Shimizu, Kanazawa & Kawahara Reference Shimizu, Kanazawa and Kawahara2019) and tilted (Tuckerman, Chantry & Barkley Reference Tuckerman, Chantry and Barkley2020) domains. In long domains the decay statistics are not memoryless but depend on the history of the turbulent bands (Xu & Song Reference Xu and Song2022). However, when the length of the turbulent bands is fixed, decay events follow a memoryless process too (Wu & Song Reference Wu and Song2025). This supports the idea that the behaviour of transiently chaotic turbulent states is universal among all subcritical shear flows in the transitional regime (Rempel, Lesur & Proctor Reference Rempel, Lesur and Proctor2010; Linkmann & Morozov Reference Linkmann and Morozov2015).
Reduced-order models of shear flows are also able to capture these decay events. Willis & Kerswell (Reference Willis and Kerswell2009) defined a reduced-order model of pipe flow by truncating their Galerkin approximation to a single azimuthal Fourier mode. Their model retained memoryless puff decay. Another example is the reduced-order model developed by Moehlis, Faisst & Eckhardt (Reference Moehlis, Faisst and Eckhardt2004), hereafter referred to as the Moehlis Faisst and Eckhardt (MFE) model. This model is a system of nine Ordinary Differential Equations (ODE) obtained from a truncated Galerkin projection of the flow between two parallel free-slip walls, driven with a wall-normal sinusoidal force (Waleffe Reference Waleffe1997). An example of a trajectory of this model is shown in figure 1(e), and lifetime statistics in figure 1(d). Lellep et al. (Reference Lellep, Prexl, Eckhardt and Linkmann2022) investigated the decay events of the model using an explainable artificial intelligence method. They successfully performed predictions at short times before decay, but observed a quick degradation of their predictions as they considered longer times before the event.
Generally there are two main challenges in predicting relaminarisation. On the one hand, there are many routes to decay in phase space (Chantry & Schneider Reference Chantry and Schneider2014; Budanur, Dogra & Hof Reference Budanur, Dogra and Hof2019). On the other hand, turbulent flows are chaotic (Lorenz Reference Lorenz1965) and the uncertainty in the determination of a flow state can result in radically different decay times. Without a correct assessment of these two issues, the ability of predictive models is severely limited and thresholds of decay are necessarily conservative.
In this paper, we examine the threshold of decay from a new probabilistic perspective: by studying predictability. We propose that the threshold of decay is the first state of a trajectory that becomes fully predictable to decay, given a finite size uncertainty (Palmer Reference Palmer2000). Uncertainties have many different origins (measurements, model, numerical method, etc.) and are unavoidable. We argue that fully characterising the predictability of decay is a prerequisite to successfully finding the threshold and causes of decay. Predictability has historically been characterised with the Lyapunov exponents (Lorenz Reference Lorenz1965). They measure the rate at which trajectories, initially separated by an infinitesimally small uncertainty
$\delta$
, separate exponentially in phase space. However, they fail to correctly measure predictability when
$\delta$
is not infinitesimally small (Boffetta et al. Reference Boffetta, Cencini, Falcioni and Vulpiani2002), or in the case of systems with several time scales (Aurell et al. Reference Aurell, Boffetta, Crisanti, Paladin and Vulpiani1996), like turbulence. They are also general measurements of a chaotic system, and not specific to a single event of interest. Another, more powerful, alternative is to directly study the evolution of the probability distribution function of possible future states of the system. The Liouville conservation equation models the evolution of probability distribution functions and has an analytical solution (Ehrendorfer Reference Ehrendorfer2006). However, it is unfeasible to define this equation for high-dimensional systems like turbulence. Recently, Jiménez (Reference Jiménez2023) integrated the probability density function of a reduced representation of a turbulent channel flow using the Perron–Frobenius operator. However, it is difficult to replicate this analysis for other turbulent flows without knowledge of appropriate reduced-order projections of the flow.
Vela-Martín & Avila (Reference Vela-Martín and Avila2024) proposed an alternative measurement of predictability. Inspired by weather forecasting (Palmer et al. Reference Palmer, Molteni, Mureau, Buizza, Chapelet and Tribbia1993), they performed massive ensembles of simulations of the two-dimensional Kolmogorov flow, and used a metric derived from information theory, the Kullback–Leibler divergence (Kullback & Leibler Reference Kullback and Leibler1951), to characterise the predictability of extreme dissipation events. In a proof of concept (Morón et al. Reference Morón, Vela-Martín and Avila2024), we used this metric to characterise predictability of decay events in the MFE model and pipe flow. In this paper, we extend our previous study and find thresholds of decay both in time and in phase space using massive ensembles and the Kullback–Leibler divergence to assess the predictability of decay. The rest of the paper is structured as follows. In § 2 we describe the MFE and puffs in pipe flow in more detail. In § 3 we describe the method we use to characterise predictability, and discuss results of predictability for the two decay events of interest. In §§ 4 and 5 we analyse MFE and puff trajectories, according to their predictability and in § 6 we draw our main conclusions.
2. Methods
2.1. Puffs in pipe flow
We consider the flow of a viscous Newtonian fluid with constant properties in a straight smooth rigid pipe of circular cross-section. The flow is incompressible and governed by the dimensionless Navier–Stokes equations
\begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + \left (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{u} = -\boldsymbol{\nabla }p + \frac {1}{\textit{Re}}\boldsymbol{\nabla} ^{2}\boldsymbol{u}+f_{p} (t )\boldsymbol{\cdot }\boldsymbol{e}_{x} \quad \text{and}\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0 . \end{align}
Here,
$\boldsymbol{u}$
is the fluid velocity,
$p$
the pressure,
$\boldsymbol{e}_{x}$
the unit vector in the axial (streamwise) direction and
$f_{p} (t )$
the pressure gradient that drives the flow;
$f_{p} (t )$
is adjusted at each time step to enforce a constant bulk velocity. All variables are rendered dimensionless using the pipe diameter (
$D$
) and the bulk velocity (
$U$
). The equations are formulated in cylindrical coordinates
$ (r,\theta ,x )$
, with velocity field components
$ (u_{r},u_{\theta },u_{x} )$
in the radial, azimuthal and axial directions, respectively. Throughout this paper we fix
$\textit{Re}=1850$
. We define
\begin{align} Q (t )= \big \langle u_{r}^{2} + u_{\theta }^{2} \big \rangle _{V} , \end{align}
as the volume-averaged kinetic energy of the cross-sectional velocity, and
\begin{align} U_{c}' (t )= \left \langle u_{\textit{HP}}\left (r=0\right )-u_{x}\left (r=0\right ) \right \rangle _{V} , \end{align}
as the volume-averaged deviation from the Hagen–Poiseuille
$u_{\textit{HP}}$
centre line velocity. We choose these two variables as they have been shown to capture the dynamics of transitional pipe flow very well (Barkley Reference Barkley2011; Barkley et al. Reference Barkley, Song, Mukund, Lemoult, Avila and Hof2015).
As a heuristic threshold, we define puff decay at
$t=t_{d}$
when
\begin{align} Q \left (t \geqslant t_{d}\right ) \leqslant 10^{-7} \text{, for } Q \left (t \lt t_{d}\right ) \gt 10^{-7} . \end{align}
We solve (2.1) numerically using our GPU-CUDA pseudo-spectral code (Morón et al. Reference Morón, Vela-Martín and Avila2024), which is publicly available in this Github. We perform Direct Numerical Simulations in a
$L_{x}=50D$
long pipe and use
$N_{r}=48$
radial points,
$M_{\theta }=96$
azimuthal points and
$M_{x}=768$
axial physical points after de-aliasing. The maximum wall Reynolds number measured is
$\textit{Re}_{\tau } \approx 70$
which results in a grid spacing in wall units of
$0.06 \lesssim \Delta r^{+} \lesssim 2.2$
,
$D \Delta \theta ^{+}/2 \approx 4.5$
and
$\Delta x^{+} \approx 9$
. The time step size is set to
$\Delta t = 0.0025$
. The resolution is optimised to reduce computing cost, while giving an accurate estimate of the lifetime, see figure 1(d).
2.2. A reduced-order model of shear flows
The MFE model (Moehlis et al. Reference Moehlis, Faisst and Eckhardt2004) has nine dimensionless time-dependent modes
$a_{\!j}$
governed by a system of nonlinear ODEs of the form
\begin{align} \frac {\mathrm{d} a_{\!j}}{\mathrm{d} t} = L_{\!j} a_{\!j} + N_{\!j} \left (\boldsymbol{a} \right ) + c_{\!j} , \end{align}
where
$L_{\!j}$
and
$N_{\!j} ( \boldsymbol{a} )$
are respectively linear and quadratic operators,
$c_{\!j}$
model constants and
$\boldsymbol{a}= (a_{1},a_{2},\ldots )^{T}$
. All the variables are normalised by the channel half-height
$h/2$
and the laminar velocity
$U_{0}$
at the walls. (See the full equations in Appendix D.)
Each variable
$a_{\!j}$
is the amplitude of one Fourier mode
$\boldsymbol{v}_{\!j}$
. The dimensionless velocity field is computed as
\begin{align} \boldsymbol{u}\left (\boldsymbol{x},t \right )= \sum _{j=1}^{9} a_{\!j}\left ( t \right )\boldsymbol{v}_{\!j}\left (\boldsymbol{x}\right ). \end{align}
The modes
$\boldsymbol{v}_{\!j}$
have a direct interpretation in the flow. The mode
$j=1$
represents the mean velocity profile. When
$a_{1} \equiv 1$
the flow is laminar, and it is turbulent if
$a_{1}\lt 1$
. The mode
$j=2$
represents velocity streaks;
$j=3$
streamwise vortices;
$j=4$
and
$j=5$
two spanwise flows;
$j=6$
and
$j=7$
two wall-normal vortices; and
$j=9$
represents a correction to the mean velocity profile. The mode
$j=8$
does not have a clear interpretation, but is a three-dimensional mode that interacts with the others.
The MFE model depends on the Reynolds number, defined as
$\textit{Re}= ({U_{0}h}/{2 \nu })$
, and the spanwise/streamwise domain sizes
$L_{z}$
and
$L_{x}$
. We fix
$\textit{Re}=400$
,
$L_{x}=4 \pi$
and
$L_{z}=2 \pi$
. We integrate the MFE model forward in time using a low-storage explicit fourth-order Runge–Kutta method, with a constant time step size
$\Delta t=0.05$
.
We define the energy of the mean profile as
\begin{align} E_{1}= \left (1-a_{1}\right )^{2}, \end{align}
and the energy of the fluctuations as
\begin{align} E_{\!j}=\sum _{j=2}^{9} a_{\!j}^{2}. \end{align}
The flow is laminar if
$E_{1}=E_{\!j}=0$
. We define turbulence decay at
$t=t_{d}$
when
\begin{align} a_{1}\left (t \geqslant t_{d}\right ) \geqslant 0.995 \text{, for } a_{1}\left (t \lt t_{d} \right )\lt 0.995 . \end{align}
3. Forecasting of turbulence decay
3.1. Massive ensembles of simulations
We run
$i=1,\ldots ,N_{b}$
MFE/pipe simulations each initialised with a different chaotic state/puff. These are our base trajectory simulations, and for each of them we save
$N_{t}$
instantaneous states at times
$t_{k}$
(every
$5$
time units for the MFE model and
$8$
for pipe flow). We run each base trajectory,
$i$
, until we observe decay according to our conservative heuristic thresholds, (2.9) for the MFE model, and (2.4) for pipe flow. We save the time at which this threshold is reached for each individual base trajectory
$i$
, and call it the time of decay
$t_{d}$
, see figures 1(c) and 1(e). Note that
$t_{d}$
is different for each base trajectory and it is always larger than the times at which we sample the trajectories
$t_{d}\gt t_{k}$
. We discarded base trajectories if
$t_{d}$
is too short,
$t_{d}\leqslant 3000$
in the MFE model and
$t_{d}\leqslant 250$
in pipe flow.
For each of the
$N_{b} \times N_{t}$
instantaneous states, we launch a massive ensemble of simulations. In the case of puffs, a simulation in an ensemble is initialised as
\begin{align} \boldsymbol{u} \left (t_{k}-t_{d}\right )= \boldsymbol{u}_{t_{k}} + \epsilon _{0} \boldsymbol{\mathcal{N}}\left (0,1\right ) , \end{align}
where
$\boldsymbol{u}_{t_{k}}$
is the velocity field of the base trajectory
$i$
at time
$t_{k}$
,
$\boldsymbol{\mathcal{N}}$
is a Gaussian noise with zero mean and unit standard deviation applied to each velocity component and at each point of the domain and
$\epsilon _{0}$
is the magnitude of the noise.
For the MFE model each simulation in the ensemble is initialised as
\begin{align} \boldsymbol{a} \left (t_{k}-t_{d}\right ) = \boldsymbol{a}_{t_{k}} + \epsilon _{0} \boldsymbol{\gamma } , \end{align}
where the nine-dimensional Gaussian random vector
$\boldsymbol{\gamma }$
has a module equal to 1. Crucially,
$\epsilon _{0}$
is always sufficiently small, so two initial conditions of the same ensemble are, from a macroscopic point of view, identical. In particular
$\epsilon _{0} \leqslant 10^{-3}$
in the MFE model, and
$\epsilon _{0} \leqslant 10^{-2}$
in pipe flow.
For each combination of flow parameters and
$\epsilon _{0}$
, we perform
$N_{{\textit{total}}}=N_{b} \times N_{t} \times N_{e}$
individual simulations, where:
$N_{b}$
is the number of base trajectories,
$N_{t}$
is the number of sampled times and
$N_{e}$
is the members per ensemble (see table 1). Unless stated otherwise we set
$\epsilon _{0} = 10^{-4}$
in the MFE model, and
$\epsilon _{0} = 10^{-2}$
in pipe flow.
Table 1. Number of simulations performed for each system of interest. Here,
$N_{b}$
is the number of base trajectories,
$N_{t}$
the number of sampled times,
$N_{e}$
the members per ensemble and
$N_{{\textit{total}}}=N_{b} \times N_{t} \times N_{e}$
the total number of individual simulations.
3.2. Measuring the predictability of decay with the Kullback–Leibler divergence
We study lifetime statistics of the ensemble members. For this purpose, we compute
$\rho$
as the probability distribution function of the ensemble lifetime,
$P$
the corresponding cumulative distribution function and
$S=1-P$
the survivor function (Lawless Reference Lawless2011). We denote
$\rho _{q}$
as the exponential lifetime distribution of puffs in pipe flow, i.e. computed with random initial conditions across phase space (Faisst & Eckhardt Reference Faisst and Eckhardt2004) and with mean lifetime
$\tau$
, see figure 1(d). While
$\rho _{q}$
is a characteristic of the dynamical system,
$\rho$
depends on the particular ensemble of interest, i.e. is the lifetime distribution conditional on being computed using a known instantaneous state plus a small noise as initial condition. We compute statistics using histograms with 20 equispaced bins. All the distributions (
$\rho$
,
$P$
and
$S$
) are horizontally shifted so the first bin corresponds to the earliest decay event observed in each ensemble.
In figure 2(b), we show the survivor functions of different ensembles, each initiated using a different sampled state at different times
$t_{k}$
of a specific base trajectory
$i=1$
. Note how the distribution of puff decay tends to be clustered close to
$t=0$
if the state used to initialise the ensemble is close to decay
$ (t- t_{d} \gtrsim -40 )$
. Here, all trajectories decay at approximately the same time as the base trajectory does, and the distribution of lifetimes inherits the Gaussian distribution of the initial conditions. For such puffs, decay is inevitable. As one uses states farther back in time as initial conditions, the resulting survivor functions gradually tend towards the expected exponential distribution. For the specific base trajectory shown in figure 2(b),
$S$
(conversely
$\rho$
) becomes indistinguishable from
$S_{q}$
(
$\rho _{q}$
) at
$t -t_{d} \lesssim -80$
.
Figure 2. Example of a puff decay event and the method we use to assess its predictability. At the top: snapshots of a turbulent puff as it decays. Grey denotes low axial velocity streaks
$u_{x}' \approx -0.4$
. Red denotes regions where
$q=u_{r}^{2}+u_{\theta }^{2} \geqslant 0.02$
. The time before decay
$t-t_{d}$
is indicated on top of each snapshot. (a) Zoom in to the puff decay trajectory shown in figure 1(c). The snapshots at the top are indicated as coloured markers in the plot. (b) Survival probability of the ensembles initialised about the instantaneous puffs shown in the snapshots. The colour of the distributions corresponds to the colour of the markers in the snapshots and in the trajectory. With a dotted line, we show the exponential distribution,
$S = 1-P_{q} \approx \exp (-t/\tau )$
. The insect is a zoom in for the ensemble at
$t-t_{d}=-40$
. The black line is the cumulative function of a Gaussian distribution,
$f=0.5-0.5 \text{erf}[b ( t - c )]$
, fitted to the data. Here,
$b$
and
$c$
are fitted parameters and
$\text{erf}$
stands for the error function. (c) Time-dependent
$K$
, (3.3), of this puff trajectory. A higher
$K$
means more predictability to decay. At
$K\rightarrow 0$
the corresponding puff is fully unpredictable. The error bars represent the uncertainty of computing
$K$
from our finite sample. The coloured markers correspond to the instantaneous states shown in the snapshots and indicated in (a).
We measure predictability as the difference between the conditional probability
$\rho$
, and the exponential
$\rho _{q}$
one. Specifically, we use the Kullback–Leibler divergence (Kullback & Leibler Reference Kullback and Leibler1951) defined as
\begin{align} K = \frac {1}{{K}_{\textit{max}}} \sum \rho \log \left ( \frac { \rho }{\rho _{q}} \right ) , \end{align}
where
${K}_{\textit{max}}$
is a parameter that normalises
$K$
to
$K \leqslant 1$
;
$K$
is minimum when
$\rho \equiv \rho _{q}$
, then
$K \rightarrow 0$
. The bigger
$K$
is, the more different the two distributions are, and therefore, the more predictable the members of the ensemble become. Since we compute statistics, and therefore
$K$
, using bins,
$K$
is maximum
$ ( K \equiv 1 )$
when the probability distribution function
$\rho$
is clustered in a single bin. In our case, this happens when all the members of the ensemble decay at nearly the same time as the base trajectory. In figure 2, this happens for the puff at
$t- t_{d} \approx -40$
. The parameter
${K}_{\textit{max}}$
is computed for this case, and depends on the number of bins we use. As we show in Appendix A our method is robust to different noise shapes, noise magnitudes and ensemble sizes. Overall,
$K$
measures the information we gain by knowing the initial state of the system up to an uncertainty
$\epsilon _{0}$
, compared with the information we have by assuming a random initial state.
For each base trajectory,
$i=1,\ldots ,N_{b}$
and sampling time
$k=1,\ldots ,N_{t}$
, we compute
$K$
of the corresponding ensemble of simulations. We assume that
$K$
is continuous in time and study its time evolution. In what follows, we refer to
$K$
as the predictability.
3.3. Temporal evolution of the predictability
In figure 2(c), we show the evolution of the predictability for a specific base trajectory. At sufficiently long times before the decay event, the lifetime distributions are similar to the exponential distribution, resulting in a fully unpredictable decay
$ ( K \rightarrow 0 )$
. As one uses as initial condition puffs that are closer to the decay event,
$K$
increases. We report that this increase does not need to be monotonic, as the decaying trajectory visits regions of phase space with varying predictability. Finally,
$K$
saturates, implying the inevitability of decay for our chosen level of uncertainty.
We repeat this analysis for all our MFE and puff base trajectories, see table 1, and compute statistics of
$K$
with respect to time, see figure 3. To compare
$K$
between different trajectories, here, and in the rest of the paper, we define
\begin{align} \Delta t_{d} = t - t_{d} , \end{align}
as our reference time frame. At the same time before a decay event, we observe
$K$
values that differ in orders of magnitude between trajectories, see figure 3. There is, however, a point in time where all the trajectories saturate,
$ ( K \rightarrow 1 )$
. In the case of the MFE this happens at
$\Delta t_{d} \approx -800$
while for the puffs it happens at
$ \Delta t_{d} \approx -40$
. This means that, for this level of uncertainty (
$\epsilon _{0}$
), this is the earliest time when one can perform a prediction to decay with almost perfect certainty, for all base trajectories studied.
Figure 3. Statistics of
$K$
(as a measurement of predictability) with respect to time, for the two decay events of interest (indicated in the title of the plot). In the two panels, the solid black line denotes the mean predictability, and the dashed yellow line a fit of this mean predictability to (3.5). The shaded region denotes the first and the last deciles of the data. The red line corresponds to the case that was unpredictable for a longer time span, and the blue line to the one that was predictable for the longest time. The error bars stand for the uncertainty in the determination of
$K$
after a bootstrapping analysis.
The mean
$K$
, averaged among all the decaying trajectories, and represented as a thick black line in the figures, is approximated here as
\begin{align} \left \langle K \right \rangle _{N_{b}} \approx A \tanh \left ( \frac {\Delta t_{d}+t_{0}}{t_{c}} \right ) - B , \end{align}
where
$A \approx B \approx 0.5$
are two dimensionless parameters, and
$t_{0}$
and
$t_{c}$
are physically relevant times. The formula fits the mean trends of
$K$
reasonably well, as seen with the dashed yellow lines in figure 3.
The time
$t_{c}$
, in (3.5), is the mean time of predictability loss, see table 2. It measures the time scale at which predictability degrades as one goes back in time with respect to the event of interest. We observe that
$t_{c}$
is approximately one order of magnitude larger than the mean Lyapunov time.
Table 2. Time scales of the two systems of interest.
The time
$t_{0}$
is the mean predictability bias. It represents the past time where, on average, predictability reaches half its maximum, and it is mostly affected by the magnitude of uncertainties
$\epsilon _{0}$
, whereas
$t_{c}$
is not, see Appendix A.3.
4. Predictability of turbulence decay in the reduced-order model
In what follows we classify trajectories and instantaneous states of the MFE model according to their predictability
$ (K )$
. As shown in figure 3,
$K$
varies by several orders of magnitudes. Specifically, unpredictable cases have very low
$K$
and therefore little weight when computing averages of predictability. To better represent the strong variations of
$K$
and appropriately weight unpredictable cases, we consider the metric
\begin{align} \kappa = \log \left ( K \right ) \end{align}
in the subsequent analyses. The presented averages of
$\kappa$
correspond to geometric averages of
$K$
and better capture the natural variation of predictability in the flow.
4.1. Classification of model trajectories according to their predictability
We here perform a conditional analysis of trajectories according to their predictability (low/high). For this purpose, we look for the first moment in time when at least
$N_{p}=1000$
trajectories have maximum
$\kappa$
. We find that this happens at
$t_{p} \approx t_{d}-1200$
. We then classify these
$N_{p}=1000$
trajectories as predictable, and conversely, the
$N_{p}=1000$
MFE trajectories that have the smallest averaged
$\kappa$
in the range
$t_{p} \lt t \lt t_{d}$
as unpredictable.
In figure 4(a), we show predictability statistics of all the MFE trajectories (black), conditioned to either predictable (blue) or unpredictable trajectories (red). We fit each mean predictability with (3.5), and analyse the resulting
$t_{0}$
and
$t_{c}$
. As expected,
$t_{0}$
is on average largest for the predictable trajectories, and smallest for the unpredictable trajectories. Between the two there is a difference of
$ \Delta t_{0}= 452.12$
. We observe that
$t_{c}$
is shorter, both for the predictable and unpredictable trajectories compared with the entire data set. This means that, inside the predictable and unpredictable groups, the mechanisms behind decay are very similar among all the trajectories in the group, and that they take place in a relatively short time.
Figure 4. Classification of MFE trajectories according to their predictability of decay. We consider three trajectory groups: all the trajectories
$N_{b}=14\,000$
(in black), only the
$N_{b}=1000$
most predictable trajectories (in blue) and only the
$N_{b}=1000$
most unpredictable trajectories (in red) at
$t_{p}\lt t\lt t_{d}$
. Thick lines denote averaged quantities among the members of the group, the shaded area denotes the first and the last deciles. (a) Predictability
$(\kappa )$
, (4.1); (b) energy of the mean profile
$E_{1}$
, (2.7); (c) energy of the fluctuations
$E_{\!j}$
, (2.8). The vertical line denotes the threshold in time
$t_{p}$
we use to differentiate between predictable and unpredictable trajectories and R stands for the temporal correlation between the plotted variable and
$\kappa$
, averaged among all the MFE trajectories.
We attempt to explain the two types of MFE decay events (predictable and unpredictable) from a fluid dynamics perspective. In figures 4(b) and 4(c), we show statistics of
$E_{1}$
(energy of the mean profile, (2.7)) and
$E_{\!j}$
(energy of the fluctuations, (2.8)) of the three groups. The temporal correlation between
$E_{1}$
(or
$E_{\!j}$
) with
$K$
, averaged among all the base trajectories, is negative. The predictable cases have a high amplitude
$E_{\!j}$
phenomenon at
$t \approx t_{p}$
. This event is subsequently followed by a flattening of the mean profile,
$E_{1} \rightarrow 0.75$
. In transitional shear flows a flat mean profile is linked with reduced turbulent fluctuations (Hof et al. Reference Hof, De Lozar, Avila, Tu and Schneider2010; Barkley et al. Reference Barkley, Song, Mukund, Lemoult, Avila and Hof2015; Kühnen et al. Reference Kühnen, Song, Scarselli, Budanur, Riedl, Willis, Avila and Hof2018). During the flattening of the profile,
$E_{\!j}$
quickly decreases. At a certain time,
$E_{1}$
starts to quickly decrease. Subsequently
$E_{\!j}$
decreases in a non-monotonic way, but it never reaches a sufficiently high amplitude to re-trigger chaos. Before decaying
$E_{\!j}$
shows some damped oscillations at
$\Delta t_{d} \gtrsim -800$
.
The unpredictable cases have a larger variance at all times
$t_{p} \lt t \lt t_{d}$
than the predictable trajectories. This means that the identification of a clear mechanism of decay is more difficult for this group. Nevertheless, unpredictable trajectories have a peak
$E_{\!j}$
event at
$\Delta t_{d} \approx -1000$
that is followed by a rapid decrease of
$E_{\!j}$
. Although on average
$E_{1}$
slightly increases after the
$E_{\!j}$
peak, it quickly decreases together with
$E_{\!j}$
. Decay events for unpredictable trajectories show a sudden collapse of all the variables at almost the same time. They all quickly decay without showing damped oscillations as
$t \rightarrow t_{d}$
.
4.2. Regions of the reduced-order model phase space according to predictability
We project our MFE database in the
$\log _{10} (E_{1} )$
and
$\log _{10} ( E_{\!j} )$
plane. We then divide the projected phase space into bins of equal sizes. In figure 5(a), we colour the bins with the maximum
$\kappa$
among all the instantaneous states in the bin. The bins coloured in red have at least one member inside that is highly predictable to decay. It is easy to recognise a big red region of phase space. This region occupies almost all the projected phase space where we have data points, showing that almost in any region of phase space there is at least one trajectory that quickly decays. In terms of dynamical systems this demonstrates that the edge is dense in the chaotic saddle (Budanur et al. Reference Budanur, Dogra and Hof2019). However, there is a region in the plot, at high values of
$\log _{10} (E_{\!j} )$
, not coloured in red. The bins in this region of the plot only have instantaneous states that have a small predictability. This means that in this region decay is unpredictable. Cases that fall in this region are expected to (on average) remain chaotic for a relatively long time span. We call this region the mixing region (MIX). We find that the formula
\begin{align} y_{{\textit{MIX}}} \approx -0.96+0.73 \boldsymbol{\cdot }\tanh \left \{ 2.1 \boldsymbol{\cdot }\left [ \log _{10}\left (E_{1}\right ) - 1.3 \right ] \right \} , \end{align}
separates reasonably well this region from the rest of the projected phase space. At a given
$\log _{10} (E_{1} )$
, states with
$\log _{10} ( E_{\!j} )\gt y_{{\textit{MIX}}}$
fall in the MIX Region.
Figure 5. Regions of the projected phase space of the MFE model depending on their predictability with respect to decay events. The sampled states of the base trajectories are projected in the
$\log _{1} (E_{1})$
and
$\log _{10} ( E_{\!j})$
reduced phase space, that is divided in
$100 \times 100$
bins of equal sizes. The colour of each bin depends on the predictability of the states inside of it: red means high predictability (high
$\kappa$
), blue low (low
$\kappa$
). (a) The colour of each bin corresponds to the maximum predictability among all the states inside of it; (b) to the mean and (c) to the minimum. The solid black lines separate regions of phase space depending on their predictability. (d) Quality of the classification of cases in the inevitable decay region as one minus the ratio of successful predictions (RSP) and its dependence on the noise magnitude
$\epsilon _{0}$
.
In figure 5(b), we colour the bins with the mean
$\kappa$
among all the instantaneous states in the bin. There is a big region of the phase space coloured in red where the average predictability of decay is very high. We call it the mean predictable region (MPR) and separate it from the rest of the projected phase space with the curve
\begin{align} y_{{\textit{MPR}}} \approx -22.1+20.4 \boldsymbol{\cdot }\tanh \left \{ -1.5 \boldsymbol{\cdot }\left [ \log _{10}\left (E_{1}\right ) - 0.5 \right ] \right \} . \end{align}
At a given
$\log _{10} (E_{1} )$
, states with
$\log _{10} ( E_{\!j} )\lt y_{{\textit{MPR}}}$
fall in the MPR. In figure 5(c), we colour the bins with the minimum
$\kappa$
among all the instantaneous states classified in that bin. The region in blue represents bins where at least one case classified in that bin is unpredictable to decay. The region in red, referred to as the Inevitable Decay Region (IDR), corresponds to bins where all the cases classified inside are highly predictable to decay. Cases found in this region will almost certainly decay in a short time. We found that the best classifier between the blue and red regions is the line
\begin{align} y_{\textit{IDR}} \approx -0.95 \log _{10}\left (E_{1}\right ) -4.4 , \end{align}
computed using a linear support vector machine (Cortes Reference Cortes1995). At a given
$\log _{10} (E_{1} )$
states with
$\log _{10} ( E_{\!j} )\lt y_{\textit{IDR}}$
are almost guaranteed to decay in a short time. We recompute the regions of phase space according to the minimum
$K$
for several uncertainty magnitudes
$\epsilon _{0}\in [4\times 10^{-5},2 \times 10^{-3} ]$
, and find that the separating line is unaffected by
$\epsilon _{0}$
. We test the quality of the classification by computing a rate of successful predictions (RSP). We define the RSP as the quotient between the number of cases classified in the IDR that have
$K \geqslant 0.99$
with respect to the total number of cases classified in the IDR. In case
$\textrm{RSP}=1$
, it means that all the cases in the IDR are fully predictable to decay. In figure 5(d) we show the behaviour of RSP with respect to
$\epsilon _{0}$
. We observe that, as
$\epsilon _{0}$
decreases, the classification significantly improves. At
$\epsilon _{0}=10^{-3}$
,
$1\,\%$
of the cases are wrongly classified in the IDR, while, at
$\epsilon _{0}=10^{-4}$
only
$0.1\,\%$
are. We note that RSP will never reach
$\textrm{RSP} \rightarrow 1$
. As we show later, there are extremely rare trajectories that can enter the IDR and subsequently remain chaotic for longer times.
4.3. A predictor of decay for the MFE model
We use the different predictability regions discussed above to develop a simple predictive model of decay events in the MFE model. The predictive model only needs two inputs: the region in the
$ ( E_{1},E_{\!j} )$
plane in which the trajectory is currently found in, and how long it has resided in that region. The model returns an indicator
${\textit{Ind}}$
that is maximum
$ ( {\textit{Ind}}=1 )$
if decay is fully predictable, and smaller
$0 \leqslant {\textit{Ind}} \lt 1$
otherwise. See Appendix C for more details. Our model correctly predicts decay (or not decay) more than
$99 \,\%$
of the time.
We tested the predictor’s performance with an extremely rare MFE trajectory (shown in figure 6), which reaches an extremely low value of the kinetic energy
\begin{align} E=E_{1} + E_{\!j} \leqslant 0.004 , \end{align}
at
$t=t_{ex}$
, but remains chaotic for at least 2000 time units thereafter. This event is similar to other extreme events observed in the MFE model, and described in Appendix B. We found this particular trajectory after trying billions of random initial conditions.
Figure 6. Rare MFE trajectory reaching
$E=E_{1}+E_{\!j}\leqslant 0.004$
at
$t=t_{ex}$
and with
$t_{d} \gt t_{ex}+2000$
. Panels show (a)
$E_{1}$
; (b)
$E_{\!j}$
; (c) our decay indicator Ind: Ind
$=0$
, no imminent decay is predicted; Ind
$\gtrsim 0.8$
decay is imminent; (d) projection of the trajectory in the (
$\log _{10} (E_{1} )$
,
$\log _{10} (E_{\!j})$
) plane. The black lines separate the regions of predictability discussed in § 4.2 (MIX with solid line, MPR with dashed, IDR with dotted). In all the plots, red means predictable, blue unpredictable and the markers help to identify the moment in time.
We observe that this rare trajectory initially behaves chaotically in a small region of phase space, and then enters the MPR. It bounces back and forth, entering and exiting the MPR, oscillating about the
$y_{{\textit{MPR}}}$
line, as Ind increases. Even though we obtained the
$y_{{\textit{MPR}}}$
line from the perspective of the mean predictability, it is representative of the dynamics of this and other trajectories (figure 12). However, instead of decaying, it then enters the MIX region. As soon as a trajectory enters the MIX region, no predictions of decay are made for a time span of
$\approx 2t_{c}$
time units. This is because, as discussed above, in this region states are highly unpredictable, and tend to decay after long times. Interestingly, after visiting the MIX region, the trajectory quickly crosses the MPR and IDR, and reaches the extremely low value of
$E \leqslant 0.004$
. Instead of decaying, however, it goes back to the MIX region and it then remains chaotic for a longer time, in agreement with our predictor (Ind
$=0$
).
This example shows that certain, albeit rare, trajectories can cross the IDR and still not subsequently decay. Lellep et al. (Reference Lellep, Prexl, Eckhardt and Linkmann2022) considered decay irreversible when
$E \leqslant 0.005$
. As we show here, trajectories can remain chaotic for long times after going below this threshold. We use this example to stress the importance of determining the threshold of decay for a level of allowed uncertainty, and the importance of characterising regions of phase space according to predictability.
5. Predictability of puff decay events
5.1. Classification of puff trajectories according to predictability
We consider as predictable (unpredictable) the
$N_{p}=10$
trajectories that have the largest (smallest) time-averaged
$\kappa$
, (4.1), in the range
$t_{p} \lt t \lt t_{d}$
, for
$t_{p}=t_{d}-120$
. In figure 7(a), we show the behaviour of
$\kappa$
for the three groups. We fit the averaged
$\kappa$
of each group to (3.5). As expected,
$t_{0}$
is larger for the predictable group, and smaller for the unpredictable group. We report here a difference of
$\Delta t_{0}=23.52$
between the two.
Figure 7. Classification of puff trajectories according to their predictability of decay. We consider three trajectory groups: all the trajectories
$N_{b}=100$
(in black), only the
$N_{b}=10$
most predictable trajectories (in blue) and only the
$N_{b}=10$
most unpredictable trajectories (in red). In all the plots the thick lines denote averaged quantities among the members of the group, and the shaded area denotes the first and the last deciles. (a) Predictability (
$\kappa$
, (4.1)); (b) the volume-averaged deviation from the centre line velocity, (2.3) and (c) the volume-averaged cross-sectional kinetic energy, (2.2). Here, R stands for the time correlation between
$\kappa$
and
$Q$
or
$U_{c}'$
, averaged among all trajectories.
In contrast to the results of the MFE model,
$t_{c}$
is much larger for the predictable group,
$t_{c} \approx 28.96$
, than for the unpredictable group
$t_{c} \approx 8.47$
. A smaller
$t_{c}$
implies that the group of trajectories have similar mechanisms of turbulence decay and that these mechanisms happen in a very short time, in an abrupt way. A large
$t_{c}$
corresponds to trajectories that undergo a more gradual decay process, and whose mechanisms of decay are more different between each other.
In figure 7(b) we show the behaviour of the defect of centre line velocity
$U_{c}'$
, (2.3), with respect to time. As the decay event is approached, the mean profile becomes more similar to the laminar one
$ (U_{c}' \rightarrow 0 )$
. We observe that the predictable trajectories have, on average, a smaller value of
$U_{c}'$
(corresponding to a more laminar-like profile) than the unpredictable ones at all times. In figure 7(c), we show the cross-sectional kinetic energy
$Q$
, (2.2), with respect to time. As the decay event approaches,
$Q$
decreases, and it is always smaller on average in the case of predictable trajectories. Both
$U_{c}'$
and
$Q$
correlate negatively with the predictability, as seen with the parameter R in the plots, and
$Q$
has a better averaged temporal correlation.
We note that, when full predictability is reached at
$\Delta t_{d} \approx -40$
, there is still a big uncertainty in the values of
$Q$
and
$U_{c}'$
a trajectory can have.
5.2. Regions of pipe flow phase space according to predictability
We project our puff database onto the
$Q$
and
$U_{c}'$
plane. We then divide the projected phase space into bins of equal sizes. In figure 8(a) we show the case of colouring each bin with the maximum
$\kappa$
among the states classified in that bin. We observe a big red region that corresponds to bins where at least one case inside is highly predictable to decay. This region covers a huge portion of the sampled phase space, particularly for any
$Q \lesssim 5\times 10^{-4}$
. It is only for large excursions of
$Q$
that puffs become highly unpredictable to decay. In figure 8(b), we colour each bin with the mean
$\kappa$
. We observe a red region that is almost completely defined by the line
$Q \lesssim 2.5\times 10^{-4}$
. Above this line predictability gradually decreases for larger values of
$Q$
, and we refer to the phase space below as the MPR of puffs.
Figure 8. Regions of the projected phase space of puffs, depending on their predictability with respect to decay events. The sampled states of the base trajectories are projected in the
$U_{c}'$
and
$Q$
reduced phase space, that is divided in
$50 \times 50$
bins of equal sizes. The colour of each bin depends on the predictability of the states inside of it: red means high predictability (high
$\kappa$
, (4.1)), blue low. (a) The colour of each bin corresponds to the maximum predictability of all the states inside of it, in (b) to the mean and in (c) to the minimum. The solid black line in (b) indicates the MPR and in (c) the IDR of puffs, (5.1).
In figure 8(c), we colour each bin with the minimum
$\kappa$
. There is a clear red region at the bottom of the plot that represents the IDR of puffs. We find that one can separate this region with a parabola
\begin{align} Q_{\textit{IDR}} \approx 0.17\,U_{c}{\prime }^{2} - 0.031\,U_{c}' + 0.0015 . \end{align}
There is an intermediate value of
$U_{c}' \approx 0.09$
, where
$Q_{\textit{IDR}}$
is minimum
$Q_{\textit{IDR}} \approx 1.51 \times 10^{-4}$
. Accordingly, there are puffs that, despite reaching such a low
$Q$
value, remain subsequently turbulent for longer times.
At high
$U_{c}' \gtrsim 0.1$
unless
$Q$
is high enough, the flow will decay. This region corresponds to flows with a relatively flat mean profile. As we also saw for the MFE model, such flows need a high enough
$Q$
to remain chaotic.
Interestingly, for
$U_{c}' \lesssim 0.07$
,
$Q$
must be high enough so the flow does not decay. At these low
$U_{c}'$
values the mean profile is very similar to the laminar one. One would expect a lower
$Q$
threshold in this region, as fluctuations are more likely to proliferate using the laminar-like mean profile. This, however, is not the case. We speculate that these
$U_{c}'$
, fluctuations, no matter how intense they are, tend to be at wall-normal locations that do not promote turbulence survival, resulting in a high
$Q_{\textit{IDR}}$
threshold.
6. Conclusions
We propose a novel, probabilistic perspective to study the threshold of decay events in transitional shear flows. Instead of directly looking for the causes behind decay events, we first characterise their predictability using massive ensembles of simulations, and then identify the flow configurations that are highly predictable to decay. Our results show that all puff and MFE trajectories saturate to a maximum predictability at a given moment in time before decay (
$\Delta t_{d} \approx -800 \text{ and }-40$
for the MFE model and puffs, respectively). This result answers part of the question we initially pose, as to what is the first moment in time at which decay becomes inevitable.
By analysing the change of predictability with respect to time, we find a characteristic time of predictability loss
$t_{c}$
. Its inverse,
$1/t_{c}$
, is the rate at which predictability degrades, on average, as one goes back in time with respect to decay events. As we show in Appendix A, this measurement is almost unaffected by the level of uncertainty
$\epsilon _{0}$
and, as we show in Appendix B, it does not depend on the (rare) event of interest: it is an intrinsic characteristic of the system. This has profound implications for the development of predictors of decay, and more generally of rare events, as the ability of any predictor will inevitably degrade exponentially on time spans of the order of
$t_{c}$
.
We identify the mechanisms by which two types of MFE trajectories decay. Trajectories that become fully predictable long before decay have a sudden flattening of the mean profile that results in an abrupt collapse of the chaotic behaviour and a decay behaviour similar to a damped oscillator. Trajectories that remain unpredictable for longer times before decay have a sudden high amplitude fluctuation event, after which turbulence quickly collapses and becomes laminar without oscillations. By projecting our data in a two-dimensional phase space, we identify different regions depending on their predictability. We report a region where all the states are highly unpredictable, a region where most cases are predictable and a region where (almost) all cases are fully predictable. The last region answers part of our initial question as to what flow configurations are fully predictable to (and therefore are a threshold of) decay. We use these regions of predictability to develop a simple predictive model of MFE decay that returns a high RSP
${\gt } 0.99$
.
We repeat the above analysis for the case of puffs in pipe flow, and consider two key variables defined by Barkley (Reference Barkley2011) and thereafter widely used to study transitional pipe flow: the kinetic energy
$Q$
of cross-flow fluctuations, and the deviation from the laminar centre line velocity,
$U_{c}'$
. We identify a region of phase space that is fully predictable to decay, and report, for the first time, a threshold of puff decay in terms of
$Q$
and
$U_{c}'$
. We observe that the threshold of decay is mostly set by a low value of
$Q$
, in line with previous studies (Nemoto & Alexakis Reference Nemoto and Alexakis2021), but depends slightly on
$U_{c}'$
. We corroborate that a pipe flow with a flat mean profile needs a high
$Q$
to survive (Kühnen et al. Reference Kühnen, Song, Scarselli, Budanur, Riedl, Willis, Avila and Hof2018). Interestingly we also observe the opposite; flows with a mean profile similar to the laminar one also need a high
$Q$
to survive. We speculate that, for these flows, the fluctuations may be found at locations far from the wall that do not promote turbulence survival.
We argue that the study of predictability can ultimately identify the causes of relaminarisation. Specifically, a flow configuration is a cause of relaminarisation to the same extent the latter is predictable from the former. It is our objective to, in future analyses, look for flow variables that better correlate with predictability and, therefore, with the causes of decay. A promising prospect is the combination of our probabilistic approach with the study of invariant solutions in shear flows (Kawahara et al. Reference Kawahara, Uhlmann and Van Veen2012). In future analyses we aim to use invariant solutions of pipe flow at the edge of chaos (e.g. Duguet, Willis & Kerswell (Reference Duguet, Willis and Kerswell2008)) as initial conditions for our ensembles. We believe that, by doing this, we can identify regions of phase space that are more/less attracting to these invariant solutions and exploit their characteristics to better understand the dynamics of decay events.
Acknowledgements
Part of this work was done during the Fifth Madrid Summer Workshop, funded by the European Research Council (ERC) under the Caust grant ERC-AdG-101018287. The members and organisers of the workshop are here greatly acknowledged. The authors would also like to thank Dr. M.P. Encinar for fruitful discussions.
Data availability statement
All the codes and data are available in Pangaea under the following link: https://doi.pangaea.de/10.1594/PANGAEA.977819. The GPU code to simulate pipe flow can be obtained in the following Github: https://github.com/Mordered/nsPipe-GPU.git.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Robustness of the method to compute predictability
A.1. Effect of the shape of the uncertainties on predictability
Here, we explore the effect of using a different type of perturbation than Gaussian noise to initialise the ensembles of simulations. In the case of puff decay, we performed a simulation with the same discretisation described in § 2.1 but at
$\textit{Re}=3000$
. We trigger turbulence in the simulation and wait until the whole domain is fully turbulent. We then save perturbed velocity fields
$\boldsymbol{u}'=\boldsymbol{u}- (0,0,u_{\textit{HP}} ( r ) )$
at
$N_{e}$
time steps, after every
$1.25$
advective time unit. We initialise each member of the ensemble with one of the scaled turbulent fields as
\begin{align} \boldsymbol{u} \left (t=0\right )= \boldsymbol{u}_{t_{k}}^{i} + \epsilon _{T} \frac {\boldsymbol{u}'}{\left | \boldsymbol{u}' \right |} . \end{align}
Note that this type of perturbation satisfies continuity of the flow. We fix
$\epsilon _{T}=2.5\times 10^{-2}$
to make the energy of the perturbation equivalent to the energy of the Gaussian perturbation at
$\epsilon _{0} \approx 10^{-2}$
. We re-compute ensembles of simulations using
$N_{t}=10$
instantaneous states of only
$N_{b}=4$
of our base trajectories, see table 1. In figure 9(a) we show predictability with respect to time depending on the type of perturbation used for the ensembles. We observe that the type of perturbation has little effect on the predictability of individual trajectories. Even after averaging over only
$N_{b}=4$
trajectories, the mean trends of predictability remain unaffected.
Figure 9. Robustness of the predictability measurement. (a) Predictability of puff decay events depending on the type of initial condition used in the ensembles of simulations. Blue corresponds to random Gaussian noise with
$\epsilon _{0} \approx 10^{-2}$
, red to a scaled fully turbulent field. The thick lines correspond to the mean predictability: averaged over
$N_{b}=100$
base trajectories, for the case of random noise, and over only
$N_{b}=4$
base trajectories for the case of scaled turbulent fields. The thin lines correspond to the same base trajectory: in red with a predictability characterised with the scaled turbulent field, in blue with random noise. (b) Predictability of decay events in the MFE model, for ensembles initialised with Gaussian noise (blue), and ensembles initialised with a scaled MFE chaotic state (red). (c) Predictability with respect to time of a single MFE decay trajectory, characterised with ensembles of different sizes
$N_{e}$
.
For the case of the MFE decay events, we save
$N_{e}$
instantaneous chaotic states of the MFE model at
$t\lt t_{d}-2000$
. We then re-scale the states to a magnitude of
$\epsilon _{0}=10^{-4}$
and use them to initialise ensembles of simulations for
$N_{b}=2000$
trajectories. In figure 9(b) we show predictability with respect to time depending on the type of perturbation used for the ensembles. We observe that the shape of the initial condition does not have an important impact on the mean trends of predictability, on the statistics or even on the predictability of individual trajectories.
A.2. Effect of using fewer members per ensemble
We also consider here the effect of using fewer members per ensemble
$N_{e}$
on the determination of the predictability of MFE decay events. We use the case
$N_{e}=1000$
as the base case scenario, and then reduce
$N_{e}$
to study the changes on
$K$
, see figure 9(c). We observe that, by decreasing
$N_{e}$
, we obtain more noisy predictability estimations, especially for small values of
$K$
. This is expected as we have a larger uncertainty in our statistics. Nevertheless we report that the mean trends of predictability are well captured for all the
$N_{e}$
considered here. Even the more limiting case of
$N_{e}=20$
is able to capture the high amplitude oscillations of predictability observed at
$N_{e}=1000$
. We repeat this analysis for additional
$N_{b}$
trajectories, compute the ensemble-averaged predictability and fit it to (3.5). We do this for all
$N_{e}$
. Although we do not show it here, we observe that, between
$N_{e}=20$
and
$N_{e}=1000$
, the fitted
$t_{c}$
only changes by
$15 \,\%$
, while between
$N_{e}=200$
and
$N_{e}=1000$
it changes by less than
$2 \,\%$
. Thus, ensembles with fewer members result in almost identical predictability characteristics as those at
$N_{e}=1000$
.
A.3. The effect of the magnitude of uncertainties
$\epsilon _{0}$
on predictability
We recompute the predictability of the
$N_{b}=14\,000$
decaying MFE trajectories, sampled at
$N_{t}=400$
time steps and perturbed with a Gaussian noise, see table 1, but with different noise magnitudes:
$\epsilon _{0}\in [4\times 10^{-5}, 2 \times 10^{-3} ]$
. In figure 10 we show the mean predictability for each
$\epsilon _{0}$
and observe that, despite a certain horizontal shift, they are almost identical.
Figure 10. Effect of the magnitude of uncertainties
$\epsilon _{0}$
on predictability of MFE decay. (a) Mean
$K$
with respect to time among
$N_{b}=14\,000$
MFE base trajectories, computed with different magnitudes of the Gaussian noise. Lighter colour means smaller
$\epsilon _{0}$
. (b) Results of fitting mean
$K$
of different
$\epsilon _{0}$
with (3.5).
We fit the mean predictability to (3.5). We report that, by decreasing
$\epsilon _{0}$
, the bias
$t_{0}$
increases
$t_{0} \propto - \log (\epsilon _{0} )$
. In our analysis,
$\epsilon _{0}$
models the level of uncertainty one has (e.g. from measurement device, measurement method, numerical method, …). As it decreases, the bias increases as expected: with smaller uncertainties, one can predict decay earlier. Note that at
$\epsilon _{0}=0$
,
$t_{0}\rightarrow \infty$
.
Although we observe a slight increase of
$t_{c}$
as
$\epsilon _{0}$
decreases, the characteristic time of predictability loss is almost unaffected by the size of
$\epsilon _{0}$
. Here, we report a difference of
$\approx 13 \,\%$
between the two limiting
$\epsilon _{0}$
cases
$\epsilon _{0}=2 \times 10^{-3}$
and
$\epsilon _{0}=4\times 10^{-5}$
. This further suggests that
$t_{c}$
is an intrinsic measurement of the dynamical system: it does not considerably depend on the coarse graining that one chooses.
Appendix B. Rare events in the reduced-order model
At
$\textit{Re}=400$
, at random times, the mode
$a_{1}$
of the MFE model becomes very large
$a_{1} \geqslant 0.8$
, giving the impression that chaos may be about to decay. However, instead of subsequently decaying, the model then remains chaotic for relatively long times. See an example of a trajectory with such an event in figure 11(a) and also in figures 1(c), 6 and 12.
Figure 11. Description and predictability of rare
$a_{1}$
events in the MFE model. (a) An MFE trajectory with an event of high
$a_{1} \geqslant 0.8$
at time
$t=t_{ex}$
and with
$t_{d} \gt t_{ex}+2000$
. (b) Statistics of this type of event in the MFE model. At
$\textit{Re}=400$
these events follow a Poisson distribution, with a mean waiting time
$\tau _{ex}\approx 2.5 \times 10^{4}$
. (c) Predictability of the MFE model to rare
$a_{1}$
events. Thin lines correspond to the more (blue) and less (red) predictable trajectories, the shaded area denotes the first and last deciles of the statistics and the thick line denotes the mean predictability. The dashed yellow line corresponds to the resultant fit of the mean
$K$
to the formula in (3.5).
We call these events ‘extreme events’ of the MFE and identify them at time
$t=t_{ex}$
according to two conditions:
-
(i) at
$t=t_{ex}$
,
$a_{1}\geqslant 0.8$
; -
(ii) and
$t_{d} \geqslant t_{ex}+2000$
.
This means that these events have a very high
$a_{1}$
magnitude, but do not quickly decay afterwards. The mean waiting time between these events follows an exponential distribution, as seen in figure 11(b). This means that, like the decay events, these events also follow a memoryless process.
As we did for the case of decay events, we compute the predictability of
$N_{b}=10\,000$
base trajectories that have an extreme event at time
$t_{ex}$
. We sample the trajectories at
$N_{t}=200$
time steps before the extreme event and launch ensembles of
$N_{e}=1000$
simulations with
$\epsilon _{0}=10^{-4}$
for each sampled state. We gather statistics of the times of extreme event recurrences for each ensemble, and compute
$K$
as in (3.3), but comparing the conditional distributions with the expected exponential distribution in figure 11(b).
We show the resultant predictability statistics in figure 11(c). As for the case of predictability of decay, predictability increases on average as one approaches the extreme event in time. We also observe differences in predictability of orders of magnitude between different trajectories. We fit the mean predictability to (3.5). We report that the resultant
$t_{c}=238.28$
is almost equal to the
$t_{c} = 264.08$
computed for the predictability of decay events, see table 2. This suggests that the predictability loss time scale is a general measurement of a chaotic system, and does not depend on the particular memoryless process of interest.
Appendix C. Description of the predictor of MFE decays
Figure 12. The MFE trajectory, and our time-dependent decay predictor. Panels show (a)
$E_{1}$
; (b)
$E_{\!j}$
; (c) our decay indicator Ind: Ind
$=0$
, no imminent decay is predicted; Ind
$\gtrsim 0.8$
decay is imminent; (d) projection of the trajectory in the (
$\log _{10} (E_{1})$
,
$\log _{10} (E_{\!j})$
) plane. The black lines separate the regions of predictability discussed in § 4.2 (MIX with solid line, MPR with dashed, IDR with dotted). In all the plots, red means predictable, blue unpredictable and the markers help to identify the moment in time.
The indicator of decay or not decay in our predictor model is computed as
\begin{align} {\textit{Ind}}=\frac {\sigma }{t_{{\textit{Ind}}}} \int ^t_{t-t_{{\textit{Ind}}}} \frac {1}{2} \left \{\text{sgn}\left [y_{{\textit{MPR}}}-\log _{10}\left (E_{\!j}\right )\right ]+1\right \} \left ( \frac {2t}{t_{{\textit{Ind}}}} \right ) \mathrm{d}t , \end{align}
where
\begin{align} \sigma & =0 \text{ if } \int ^t_{t-t_{{\textit{MIX}}}} \frac {1}{2} \left \{\text{sgn}\left [\log _{10}\left (E_{\!j}\right ) - y_{{\textit{MIX}}}\right ]+1\right \} \mathrm{d}t \gt 0 \text{, and } \end{align}
\begin{align} \sigma &=1 \text{ otherwise,} \end{align}
with
$\text{sgn}$
being the signum function. The parameter
$t_{{\textit{Ind}}}$
is computed as
\begin{align} t_{{\textit{Ind}}}= \min \left (t_{{\textit{MPR}}},t-t_{\sigma }\right ) , \end{align}
with
$t_{\sigma }$
being the last time step where
$\sigma \lt 1$
.
Here,
$y_{{\textit{MIX}}}$
and
$y_{{\textit{MPR}}}$
are the delimiters of the mixing, (4.2), and MPRs, (4.3); and
$t_{{\textit{MPR}}}$
and
$t_{{\textit{MIX}}}$
are two parameters of the predictor. Here, we tune
$t_{{\textit{MPR}}}=100$
and
$t_{{\textit{MIX}}}=400$
so they have the same order of magnitude as
$t_{c}$
.
When the MFE enters the MPR the indicator
${\textit{Ind}}\gt 0$
, and increases in magnitude the longer the time the trajectory spends there. As soon as the time spent in the MPR is
$t \geqslant t_{{\textit{MPR}}}$
, the indicator saturates to
$1$
. The predictor also checks if the trajectory enters the MIX region. In this region decay is highly unpredictable. As soon as the trajectory enters the MIX region,
${\textit{Ind}}=0$
and no predictions are made for
$t_{{\textit{MIX}}}$
time units.
As long as
${\textit{Ind}}\gt 0$
, the predictor also outputs a forecast decay time
\begin{align} \Delta t_{{\textit{decay}}} \left (t \right )=\frac {4{\textit{Re}}}{\pi ^{2}}\log \left ( \frac {1-a_{1}\left (t \right )}{1-0.995} \right ) . \end{align}
This represents the forecasted remaining lifetime before decay. This decay time is derived from the equation of the variable
$a_{1}$
(Moehlis et al. Reference Moehlis, Faisst and Eckhardt2004) by neglecting the nonlinear terms
\begin{align} \frac {\mathrm{d}a_{1}}{\mathrm{d}t}= \frac {\pi ^{2}}{4{\textit{Re}}}\left (1 -a_{1}\right ) , \end{align}
integrating and assuming that decay is completed at
$a_{1} \approx 0.995$
(our heuristic threshold).
See in figure 12(c) an example of the time-dependent Ind for the MFE base trajectory shown in figure 1(a). In figure 12(d) we show the position of the trajectory in the projected phase space. We observe how the trajectory starts in the right top corner of the plot, where we observe most trajectories spend the majority of their time. The trajectory then starts to approach the MPR. As soon as the trajectory enters the MPR, Ind increases (represented by a change in colour of the line).
The trajectory then enters the MIX region, and all predictions are discarded for a certain time span. Finally, the trajectory enters again the MPR, and Ind increases monotonically to 1, as the trajectory finally decays.
C.1. Ability of the predictor
To assess the quality of our predictor, we define the function
\begin{align} g\left ( t \right )=\exp \left [- \left (\frac {\Delta t_{{\textit{decay}}} (t ) - t_{d} + t }{t_{{\textit{tol}}}} \right )^{6} \right ] , \end{align}
where
$t_{{\textit{decay}}}$
is the forecasted remaining lifetime,
$t_{d}$
the actual time to decay and
$t_{{\textit{tol}}}=800$
. We set this tolerance to alleviate the strong assumption we make in (C6), where we assume a viscous decay. The above function is equal to
$\approx 1$
at
$\left | \Delta t_{{\textit{decay}}} (t ) - t_{d} + t \right | \leqslant 400$
and smaller for any other time span:
$g \approx 0.5$
at
$\left | \Delta t_{{\textit{decay}}} (t ) - t_{d} + t \right | \leqslant 750$
. Note that, as shown in figure 3(b), all MFE trajectories reach maximum predictability at
$t \approx 800$
time units before decay
$t_{d}$
.
We compute the error of the predictor with respect to time as the product
\begin{align} \epsilon _{1}={\textit{Ind}} \boldsymbol{\cdot }\left (1-g\right ) . \end{align}
We compute the false positives of a base trajectory as
\begin{align} {FP} = \frac {1}{\Delta t_{{\textit{Ind}}}} \int _{0}^{\Delta t_{{\textit{Ind}}}} \epsilon _{1} \mathrm{d} t , \end{align}
where
$\Delta t_{{\textit{Ind}}}$
are all the time steps where
${\textit{Ind}}\gt 0$
. We compute the ratio of successful predictions as
\begin{align} \textrm{RSP} = 1- \frac {1}{T} \int _{0}^{T} \epsilon _{1} \mathrm{d} t , \end{align}
to account for the amount of time the predictor is successfully predicting decay (or not decay). Here,
$T=2100$
is the total time span of each base trajectory. We compute the error and the RSP for all our
$N_{b}=14\,000$
MFE base trajectories.
We report that the model returns a mean
$\left \langle {FP} \right \rangle _{b}=7.16 \,\%$
, out of which only
$0.77 \,\%$
are full false positives. Full false positives are cases where
$\epsilon _{1} \geqslant 0.95$
. We report that the model returns a mean
$\left \langle \textrm{RSP} \right \rangle _{b,t}=99.73 \,\%$
, as it predicts correctly most of the time that the MFE is, or is not, about to decay.
Appendix D. The MFE equations
We include in this appendix the MFE equations, as originally derived by Moehlis et al. (Reference Moehlis, Faisst and Eckhardt2004). Let
\begin{align} \beta = \frac {\pi }{2}{, }\alpha =\frac {2 \pi }{L_{x}} {, }\gamma =\frac {2 \pi }{L_{z}}{, }k_{\alpha \gamma }=\sqrt {\alpha ^{2} + \gamma ^{2}}{, }k_{\beta \gamma }=\sqrt {\beta ^{2} + \gamma ^{2}}\text{ and }k_{\alpha \beta \gamma }=\sqrt {\alpha ^{2} + \beta ^{2} + \gamma ^{2}} . \end{align}
The equations for each mode read
\begin{align} \frac {\mathrm{d}a_{1}}{\mathrm{d}t} &= \left ( \frac {\beta ^2}{\textit{Re}} \right ) (1 - a_1) + \sqrt {1.5} \beta \gamma \left ( \frac {a_2 a_3}{k_{\beta \gamma }} - \frac {a_6 a_8}{k_{\alpha \beta \gamma }} \right ) ,\end{align}
\begin{align} \frac {\mathrm{d}a_{2}}{\mathrm{d}t} &=- \frac {a_2}{\textit{Re}} \left ( \frac {4\beta ^2}{3} + \gamma ^2 \right ) + a_4 a_6 \frac {\sqrt {50/27} \gamma ^2}{k_{\alpha \gamma }} - a_5 a_7 \frac {\gamma ^2}{\sqrt {6} k_{\alpha \gamma }}\nonumber \\ &\quad - a_5 a_8 \frac {\alpha \beta \gamma }{\sqrt {6} k_{\alpha \gamma } k_{\alpha \beta \gamma }} - a_1 a_3 \frac {\sqrt {1.5} \beta \gamma }{k_{\beta \gamma }} - a_3 a_9 \frac {\sqrt {1.5} \beta \gamma }{k_{\beta \gamma }}, \end{align}
\begin{align} \frac {\mathrm{d}a_{3}}{\mathrm{d}t} &= - a_3 \frac {\beta ^2 + \gamma ^2}{\textit{Re}} + a_4 a_7 \frac {2 \alpha \beta \gamma }{\sqrt {6} k_{\alpha \gamma } k_{\beta \gamma }} + a_5 a_6 \frac {2 \alpha \beta \gamma }{\sqrt {6} k_{\alpha \gamma } k_{\beta \gamma }}\nonumber \\ &\quad + a_4 a_8 \frac {\beta ^2 (3 \alpha ^2 + \gamma ^2) - 3 \gamma ^2 (\alpha ^2 + \gamma ^2)}{\sqrt {6} k_{\alpha \gamma } k_{\beta \gamma } k_{\alpha \beta \gamma }} ,\end{align}
\begin{align} \frac {\mathrm{d}a_{4}}{\mathrm{d}t} &= - a_4 \frac {3\alpha ^2 + 4\beta ^2}{3{\textit{Re}}} - a_1 a_5 \frac {\alpha }{\sqrt {6}} - a_2 a_6 \frac {10 \alpha ^2}{3 \sqrt {6} k_{\alpha \gamma }} - a_3 a_7 \frac {\sqrt {1.5} \alpha \beta \gamma }{k_{\alpha \gamma } k_{\beta \gamma }}\nonumber \\ &\quad - a_3 a_8 \frac {\sqrt {1.5} \alpha ^2 \beta ^2}{k_{\alpha \gamma } k_{\beta \gamma } k_{\alpha \beta \gamma }} - a_5 a_9 \frac {\alpha }{\sqrt {6}} ,\end{align}
\begin{align} \frac {\mathrm{d}a_{5}}{\mathrm{d}t} &= - a_5 \frac {\alpha ^2 + \beta ^2}{\textit{Re}} + a_1 a_4 \frac {\alpha }{\sqrt {6}} + a_2 a_7 \frac {\alpha ^2}{\sqrt {6} k_{\alpha \gamma }} - a_2 a_8 \frac {\alpha \beta \gamma }{\sqrt {6} k_{\alpha \gamma } k_{\alpha \beta \gamma }} + a_4 a_9 \frac {\alpha }{\sqrt {6}}\nonumber \\ &\quad + a_3 a_6 \frac {2 \alpha \beta \gamma }{\sqrt {6} k_{\alpha \gamma } k_{\beta \gamma }} ,\end{align}
\begin{align} \frac {\mathrm{d}a_{6}}{\mathrm{d}t} &= - a_6 \frac {3\alpha ^2 + 4\beta ^2 + 3\gamma ^2}{3{\textit{Re}}} + a_1 a_7 \frac {\alpha }{\sqrt {6}} + a_1 a_8 \frac {\sqrt {1.5} \beta \gamma }{k_{\alpha \beta \gamma }} + a_2 a_4 \frac {10\bigl(\alpha ^2 - \gamma ^2\bigr)}{3 \sqrt {6} k_{\alpha \gamma }} \nonumber \\ &\quad - a_3 a_5 \frac {2 \alpha \beta \gamma }{\sqrt {1.5} k_{\alpha \gamma } k_{\beta \gamma }} + a_7 a_9 \frac {\alpha }{\sqrt {6}} + a_8 a_9 \frac {\sqrt {1.5} \beta \gamma }{k_{\alpha \beta \gamma }}, \end{align}
\begin{align} \frac {\mathrm{d}a_{7}}{\mathrm{d}t} &= - a_7 \frac {\alpha ^2 + \beta ^2 + \gamma ^2}{\textit{Re}} - a_1 a_6 \frac {\alpha }{\sqrt {6}} - a_6 a_9 \frac {\alpha }{\sqrt {6}} + a_2 a_5 \frac {\gamma ^2 - \alpha ^2}{\sqrt {6} k_{\alpha \gamma }} + a_3 a_4 \frac {\alpha \beta \gamma }{\sqrt {6} k_{\alpha \gamma } k_{\beta \gamma }} ,\end{align}
\begin{align} \frac {\mathrm{d}a_{8}}{\mathrm{d}t} &= - a_8 \frac {\alpha ^2 + \beta ^2 + \gamma ^2}{\textit{Re}} + a_2 a_5 \frac {2 \alpha \beta \gamma }{\sqrt {6} k_{\alpha \gamma } k_{\alpha \beta \gamma }} + a_3 a_4 \frac {\gamma ^2 \bigl(3 \alpha ^2 - \beta ^2 + 3 \gamma ^2\bigr)}{\sqrt {6} k_{\alpha \gamma } k_{\beta \gamma } k_{\alpha \beta \gamma }} ,\end{align}
\begin{align} \frac {\mathrm{d}a_{9}}{\mathrm{d}t} &= - a_9 \frac {9 \beta ^2}{\textit{Re}} + a_2 a_3 \frac {\sqrt {1.5} \beta \gamma }{k_{\beta \gamma }} - a_6 a_8 \frac {\sqrt {1.5} \beta \gamma }{k_{\alpha \beta \gamma }}. \end{align}
Open access
