2.1. Time evolution of uncertainty production using Navier–Stokes equations

Ge et al. (Reference Ge, Rolland and Vassilicos2023) had derived the governing equation of the velocity field difference $\boldsymbol{w}$ from the incompressible Navier–Stokes equations:

(2.1) \begin{align} \partial _{t}w_{i}+\left (u_{\!j}+w_{\!j}\right )\partial _{\!j}w_{i}+w_{\!j}\partial _{\!j}u_{i} &=-\partial _{i}h+\nu \partial _{\!j}\partial _{\!j}w_{i}+g_{i}, \nonumber\\ \partial _{i}w_{i} &= 0, \end{align}

where $\boldsymbol{u}=\boldsymbol{u}^{(1)}$ is the reference flow velocity. In the following discussion, we omit the superscripts indicating the flow field to which quantities belong and fluid quantities are assumed to be taken from the reference flow field unless otherwise stated. Here, $\partial _{t}$ denotes the partial derivative with respect to time and $\partial _{i}$ denotes the partial derivative with respect to the $i$ th spatial coordinate. In (2.1), $h=p_{r}^{(2)}-p_{r}^{(1)}$ denotes the pressure-difference field (implicitly divided by the fluid density) and $\boldsymbol{g}=\boldsymbol{f}^{(2)}-\boldsymbol{f}^{(1)}$ represents the external forcing difference. In the present work, we assume that there is no external input from the forcing, i.e. $\boldsymbol{g}=0$ .

We first derive the time evolution of the stress $\unicode{x1D64D}$ , with each component defined as $R_{\textit{ij}}\equiv w_{i}w_{\!j}$ . By calculating the outer product of (2.1) with $\boldsymbol{w}$ and then add the transpose of this outer product, we obtain

(2.2) \begin{equation} \begin{aligned} \partial _{t}R_{\textit{ij}} &+\left (u_{k}+w_{k}\right )\partial _{k}R_{\textit{ij}} +R_{\textit{jk}}\partial _{k}u_{i} +R_{\textit{ik}}\partial _{k}u_{\!j} \\ &= -\left (w_{\!j}\partial _{i}+w_{i}\partial _{\!j}\right )h +\nu \partial _{k}\partial _{k}R_{\textit{ij}} -2\nu \left (\partial _{k}w_{i}\right )\left (\partial _{k}w_{\!j}\right )\!. \end{aligned} \end{equation}

Considering that the local uncertainty-production term according to (1.1) is $P_{\varDelta }=-w_{i}S_{\textit{ij}}w_{\!j}=-S_{\textit{ij}}R_{\textit{ij}}$ (Ge et al. Reference Ge, Rolland and Vassilicos2023), where the strain-rate tensor $\unicode{x1D64E}$ is symmetrical with each component equal to $S_{\textit{ij}}= (\partial _{\!j}u_{i}+\partial _{i}u_{\!j} )/2$ , we know the time evolution of the strain rate tensor is

(2.3) \begin{equation} \partial _{t}S_{\textit{ij}}+u_{k}\partial _{k}S_{\textit{ij}}+S_{\textit{ik}}S_{\textit{kj}}+\varOmega _{\textit{ik}}\varOmega _{\textit{kj}}=-P_{\textit{ij}}+\nu \partial _{k}\partial _{k}S_{\textit{ij}}+F_{\textit{ij}}, \end{equation}

where $\boldsymbol{\varOmega }$ is the anti-symmetrical rotation-rate tensor with each component defined as $\varOmega _{\textit{ij}}= (\partial _{\!j}u_{i}-\partial _{i}u_{\!j} )/2$ . Here, $P_{\textit{ij}} \equiv \partial _{i}\partial _{\!j}p_{r}$ is the pressure Hessian tensor and $F_{\textit{ij}} \equiv (\partial _{i}f_{\!j}+\partial _{\!j}f_{i} )/2$ is the effect of the forcing on the strain rate variation. Equations (2.2) and (2.3) give the tempo-spatial evolution of local uncertainty-production term:

(2.4) \begin{eqnarray} \partial _{t}P_{\varDelta }+\left (u_{k}+w_{k}\right )\partial _{k}P_{\varDelta }=&-&R_{\textit{ij}}w_{k}\partial _{k}S_{\textit{ij}}+\underbrace {R_{\textit{ij}}\left [\left (S_{\textit{ik}}\varOmega _{\textit{kj}}-\varOmega _{\textit{ik}}S_{\textit{kj}}\right )+\varOmega _{\textit{ik}}\varOmega _{\textit{kj}}+3S_{\textit{ik}}S_{\textit{kj}}\right ]}_{\text{(i)}}\nonumber \\ &+&\underbrace {S_{\textit{ij}}\left (w_{\!j}\partial _{i}+w_{i}\partial _{\!j}\right )h+R_{\textit{ij}}P_{\textit{ij}}}_{\text{(ii)}}\nonumber \\ &+&\underbrace {\nu \left [\partial _{k}\partial _{k}P_{\varDelta }+2S_{\textit{ij}}\left (\partial _{k}w_{i}\right )\left (\partial _{k}w_{\!j}\right )+2\left (\partial _{k}S_{\textit{ij}}\right )\left (\partial _{k}R_{\textit{ij}}\right )\right ]}_{\text{(iii)}}\nonumber \\ &-&\underbrace {R_{\textit{ij}}F_{\textit{ij}}}_{\text{(iv)}}. \end{eqnarray}

While the uncertainty energy $E_{\varDelta }$ quantifies the similarity or correlation between two flow fields at a given instant, the uncertainty-production term $P_{\varDelta }$ instead characterises the temporal ‘stability’ of points in the flow field. Equation (2.4), which describes the space and time evolution of the uncertainty-production term, is complex and provides limited intuitive insight. In Navier–Stokes turbulence, $P_{\varDelta }$ evolves due to four factors: (i) inertial effects caused locally by strain and vorticity in the flow, as shown in the first line on the right-hand-side of (2.4); (ii) pressure effects, including the global difference between two flows (characterised by $h$ ) and the global influence of the reference flow (characterised by $P_{\textit{ij}}$ ), as shown in the second line on the right-hand-side of (2.4); (iii) viscosity effects; and (iv) external forcing. Here, we focus on the exponential growth regime of the uncertainty energy $\langle E_{\varDelta }\rangle$ . This regime is observed only when the perturbation norm is infinitesimal, i.e. $|\boldsymbol{w}|\ll |\boldsymbol{u}|$ (Ruelle Reference Ruelle1981). In this limit, the contribution $-R_{\textit{ij}}w_k\partial _k S_{\textit{ij}}$ on the right-hand side of (2.4) is much smaller than the retained production terms and is therefore neglected in the following analysis.

Given the complexity of (2.4), in the present work, we focus primarily on the local inertial contribution to $P_{\varDelta }$ . This contribution originates from the coupling between the perturbation field $\boldsymbol{w}$ and the local velocity-gradient field. It therefore provides a direct connection between $P_{\varDelta }$ and local kinematic quantities, such as strain, rotation and the associated flow topology. This connection offers an intuitive way to relate uncertainty production to the local structure of the flow, and provides a basis for identifying and visualising the regions that are most susceptible to subsequent uncertainty growth.

However, the pressure-induced and viscosity-induced contributions involve different mechanisms, and are not analysed in detail here. The pressure is intrinsically non-local in incompressible turbulence, since it is determined by the global velocity gradient through the Poisson equation. As a result, its contribution at a spatial point cannot be interpreted from the local velocity-gradient tensor or from the local sampling used in the following sections. The viscosity-induced contribution is associated with small-scale gradients and dissipative regularisation of the perturbation field. A detailed interpretation of this term would require a separate analysis of dissipative and scale-dependent effects. For these reasons, the present work concentrates mainly on the inertial mechanism, with the aim of clarifying how local flow structures contribute to the growth of uncertainty production.

2.2. Inertial effect on the uncertainty production

According to (2.4), the inertial effect on uncertainty production can be decomposed into the interaction of velocity-difference stress with three different symmetrical tensors. The tensor $[\unicode{x1D64E},\boldsymbol{\varOmega }]$ (with components $[\unicode{x1D64E},\boldsymbol{\varOmega }]_{\textit{ij}}=S_{\textit{ik}}\varOmega _{\textit{kj}}-\varOmega _{\textit{ik}}S_{\textit{kj}}$ ) represents the vortex stretching rate with zero trace. Using the principal axes of $[\unicode{x1D64E},\boldsymbol{\varOmega }]$ as a local orthonormal reference frame, we can write

(2.5) \begin{equation} R_{\textit{ij}}\left (S_{\textit{ik}}\varOmega _{\textit{kj}}-\varOmega _{\textit{ik}}S_{\textit{kj}}\right )=\varLambda ^{a}_{1}(w_{1}^{a})^{2}+\varLambda ^{a}_{2}(w_{2}^{a})^{2}+\varLambda ^{a}_{3}(w_{3}^{a})^{2}, \end{equation}

where $\varLambda ^{a}_{1}$ , $\varLambda ^{a}_{2}$ , $\varLambda ^{a}_{3}$ are the eigenvalues of $[\unicode{x1D64E},\boldsymbol{\varOmega }]$ and $w_{1}^{a}$ , $w_{2}^{a}$ , $w_{3}^{a}$ are the components of the velocity-difference vector projected on the corresponding principal axes. Incompressibility forces $[\unicode{x1D64E},\boldsymbol{\varOmega }]$ are to be traceless, i.e. $\varLambda ^{a}_{1}+\varLambda ^{a}_{2}+\varLambda ^{a}_{3}=0$ . Defining the order of eigenvalues as $\varLambda ^{a}_{1}\leq \varLambda ^{a}_{2}\leq \varLambda ^{a}_{3}$ , we must have $\varLambda ^{a}_{1}\lt 0$ representing local vortex compression and $\varLambda ^{a}_{3}\gt 0$ representing local vortex stretching, while the sign of the intermediate eigenvalue is uncertain. The important point which can now be made based on (2.5) is that uncertainty-production is increased by vortex stretching ( $\varLambda ^{a}_{3}\gt 0$ ) and is decreased by vortex compression ( $\varLambda ^{a}_{1}\lt 0$ ). The total contribution of $R_{\textit{ij}} (S_{\textit{ik}}\varOmega _{\textit{kj}}-\varOmega _{\textit{ik}}S_{\textit{kj}} )$ to $P_{\varDelta }$ depends on the alignment of uncertainty energy with the principal axes of $[\unicode{x1D64E},\boldsymbol{\varOmega }]$ .

Vortex stretching increases uncertainty production. At the same time, it also enhances vorticity. The vorticity term $R_{\textit{ij}}\varOmega _{\textit{ik}}\varOmega _{\textit{kj}}$ , however, reduces uncertainty production. Since $\boldsymbol{\mathsf{\varOmega }}$ is an anti-symmetrical tensor, $\boldsymbol{\varOmega }\boldsymbol{\cdot }\boldsymbol{\varOmega }$ (with components $\varOmega _{\textit{ik}}\varOmega _{\textit{kj}}$ ) is symmetrical with three non-positive eigenvalues, defined to be $\varLambda ^{b}_{1}\leq \varLambda ^{b}_{2}\leq \varLambda ^{b}_{3}\leq 0$ . Following (2.5), we rewrite $R_{\textit{ij}}\varOmega _{\textit{ik}}\varOmega _{\textit{kj}}$ on the principal axes of $\boldsymbol{\varOmega }\boldsymbol{\cdot }\boldsymbol{\varOmega }$ as

(2.6) \begin{equation} R_{\textit{ij}}\varOmega _{\textit{ik}}\varOmega _{\textit{kj}}=\varLambda ^{b}_{1}(w_{1}^{b})^{2}+\varLambda ^{b}_{2}(w_{2}^{b})^{2}+\varLambda ^{b}_{3}(w_{3}^{b})^{2}, \end{equation}

$w_{1}^{b}$ , $w_{2}^{b}$ , $w_{3}^{b}$ are the components of the velocity-difference vector projected on the principal axes of $\boldsymbol{\varOmega }\boldsymbol{\cdot }\boldsymbol{\varOmega }$ . We have $R_{\textit{ij}}\varOmega _{\textit{ik}}\varOmega _{\textit{kj}}\leq 0$ , meaning that the enstrophy decreases uncertainty production.

We now check the tensor $\unicode{x1D64E}\boldsymbol{\cdot }\unicode{x1D64E}$ (with components $S_{\textit{ik}}S_{\textit{kj}}$ ), representing the magnitude of strain rate. Given that $\unicode{x1D64E}$ is symmetrical, the eigenvalues of $\unicode{x1D64E}\boldsymbol{\cdot }\unicode{x1D64E}$ , whose order is defined as $0\leq \varLambda ^{c}_{1}\leq \varLambda ^{c}_{2}\leq \varLambda ^{c}_{3}$ , are all positive. Following (2.5), we have

(2.7) \begin{equation} R_{\textit{ij}}S_{\textit{ik}}S_{\textit{kj}}=\varLambda ^{c}_{1}(w_{1}^{c})^{2}+\varLambda ^{c}_{2}(w_{2}^{c})^{2}+\varLambda ^{c}_{3}(w_{3}^{c})^{2}, \end{equation}

where $\boldsymbol{w}^{c}$ are the velocity-difference vector rewritten on the principal axes of $\unicode{x1D64E}\boldsymbol{\cdot }\unicode{x1D64E}$ . We have that $R_{\textit{ij}}S_{\textit{ik}}S_{\textit{kj}}\geq 0$ , suggesting all the strains of flow, no matter compression or stretching, increase the uncertainty production.

Actually, (2.4) has such a complex form that we cannot use it to localise or predict the extreme events of uncertainty production. However, the connection between uncertainty production and the strain rate or vorticity of the reference field, as shown in (2.4), suggests the possibility of relating the uncertainty-production term to the flow field’s topographical structures, which are defined by strain rate and vorticity. This approach is presented in § 5, and will help us analyse the DNS data to intuitively understand and localise extreme events of uncertainty production. To predict the extreme uncertainty-production events, considering the complexity of the deterministic governing equation (2.4), we will use the committor function for probabilistic forecasting. The committor function is estimated using DNS data and the method proposed by Lucente et al. (Reference Lucente, Rolland, Herbert and Bouchet2022b ), as shown in § 3.

3.1. Position of the problem

The central motivation for our study comes from the observation by Ge et al. (Reference Ge, Rolland and Vassilicos2023) that rare extreme production events contribute to a large proportion of the spatially averaged uncertainty production. Our goal is therefore not only to isolate these high-production events, but also to diagnose which flow properties are most relevant to their occurrence and to the increase of uncertainty.

For this purpose, one could follow an a posteriori approach. By this, we mean sampling histograms of the uncertainty production (figures 1 and 4), then determining the weight of tails where $P_\Delta$ is larger than its average by several standard deviations $\sigma _{P_\Delta }$ . From that point, it is sometimes possible to analyse the flow states corresponding to these tails and, in some cases, trace the time evolution backward to determine what flow configuration, and at what time, led to these events. While conceptually simple, this approach has several limitations related to the rarity of events and the difficulty to sample them. Moreover, following this approach, it is not straightforward to forecast the probability of an extreme event in the future given the current flow state.

Figure 1.

An example of probability density function (p.d.f.) of the local uncertainty-production term $P_{\varDelta }$ obtained as explained in § 4. The sets $\mathcal{B}$ and $\mathcal{A}$ are schematically represented in this plot: $\mathcal{B}$ contains rare extreme events in the positive $P_{\varDelta }$ tail of the p.d.f. and $\mathcal{A}$ contains all the negative production events. The red line is the p.d.f. of $P_{\varDelta }$ normalised by the probability at $P_{\varDelta }=0$ and the standard deviation $\sigma _{P_{\varDelta }}$ of $P_{\varDelta }$ estimated during the exponential growth of $\langle E_{\varDelta }\rangle$ . The working definition of $\sigma _{P_{\varDelta }}$ is given in figure 4.

To address this predictive problem, we introduce the committor. We write a realisation of the dynamics as $z(t)$ , with initial condition $z_0=z(t=0)$ . We define $\mathcal{B}$ as the set of states $z$ corresponding to an extreme event and $\mathcal{A}$ as the set of states $z$ corresponding to normal conditions, for instance, using thresholds on one p.d.f. tail (shaded in purple in figure 1) and around the mode (shaded in blue in figure 1), respectively. These two sets should be disjoint. We then define the first-hitting times of the sets $\mathcal{A}$ and $\mathcal{B}$ for a realisation $z$ evolving in phase space as

(3.1) \begin{equation} \begin{aligned} T_{\mathcal{A}} &\equiv \inf \left \{ t \ge 0 : z(t)\in \mathcal{A} \right \}\!, \\ T_{\mathcal{B}} &\equiv \inf \left \{ t \ge 0 : z(t)\in \mathcal{B} \right \}\!. \end{aligned} \end{equation}

These are simply the first times at which the system reaches either of the two sets.

Although the underlying Navier–Stokes dynamics considered in this paper is deterministic, the phase-space variable $z$ used here is a coarse-grained representation built from only a finite number of physically relevant observables. Consequently, one point $z_0$ in this reduced phase space may correspond to multiple underlying flow states, which can evolve along different trajectories and reach $\mathcal{A}$ or $\mathcal{B}$ at different times. The probability in the committor therefore arises from this coarse-grained description. The committor is then defined as the probability of reaching set $\mathcal{B}$ before set $\mathcal{A}$ , given the starting point $z_0$ :

(3.2) \begin{equation} q(z_0)=P\left (T_{\mathcal{B}}\lt T_{\mathcal{A}} \mid z(t=0)=z_0\right )\!. \end{equation}

In the context of rare and extreme events, $z_0$ is often close to $\mathcal{A}$ , and the probability to reach the extreme set $\mathcal{B}$ is small compared with $1$ . This probability has been an object of study for more than eight decades in chemical kinetics (Onsager Reference Onsager1938) and the term committor was coined approximately twenty years ago (E et al. Reference E, Ren and Vanden-Eijnden2005; E & Vanden-Eijnden Reference E and Vanden-Eijnden2006; Prinz et al. Reference Prinz, Held, Smith and Noé2011). This function can be defined in a deterministic (where it is either $0$ or $1$ ) or in a probabilistic context (where it can span all values between $0$ and $1$ ) (Lucente et al. Reference Lucente, Herbert and Bouchet2022a ). If the committor is known, we can then study extreme events a priori. This probability depends on the initial condition $z_0$ , it can be used to perform forecasts, and it is possible to construct new composites and to perform averages of states ‘ $z_0$ ’ conditioned on their probability to evolve into an extreme (Miloshevich et al. Reference Miloshevich, Cozian, Abry, Borgnat and Bouchet2023). The features of these flow states can be studied to try to highlight precursors of given extremes.

While the system we study is entirely deterministic, taking a probabilistic approach by using the committor function is practically simpler. First, choosing to estimate a committor from a deterministic point of view requires a fine scanning of all points in phase space with minimal error bars on event properties; this is in contradiction with our object of study, where there is by default a perturbation to velocity fields. Second, if we make some additional reasonable assumptions on our process, the committor equation (3.2) becomes strongly structured because it follows a given algebraic or differential equation.

Indeed, this is the case if we assume that the time evolution can be modelled by a Markov chain, that is to say, a stochastic process where each random change between $t$ and $t+\Delta t$ is solely determined by the state of the system at $t$ . This requires two first reasonable approximations: the representation of the time evolution of the system by a process discrete in time with step $\Delta t$ (larger or equal to $\textrm {d}t$ , the time step of our simulations) as well as rapid enough decorrelations in time. In that case, the starting point for the definition and computation of the committor is the conjugate probability $p(z(t)\,|\,z(t=0)=z_0)$ of being in a state $z(t)$ at some time $t\gt 0$ conditioned on having been in state $z_0$ at time $t=0$ , i.e. $z(t=0)=z_0$ . Note that while $z$ is the independent variable for $p$ , $z_0$ is the independent variable for $q$ . This probability $p$ contains the information on all possible realisations of the process at a given time.

In our discrete-in-time Markov chain, the advance in time of $p$ is given by the application of a linear operator $G^\dagger$ , i.e.

(3.3) \begin{equation} p(t+\Delta t)=G^\dagger (p(t)), \end{equation}

which is known as a Chapman–Kolmogorov equation (Gardiner Reference Gardiner2009; Friedrich et al. Reference Friedrich, Peinke, Sahimi and Reza Rahimi Tabar2011). The operator $G^\dagger$ contains the transition probabilities from any state $z$ to any other state $z'$ during a time $\Delta t$ , and its action is expressed as the sum of transition probabilities times probabilities (an actual sum for discrete systems or an integral for continuous ones). We chose this notation because the committor is then a steady solution of the affine adjoint problem

(3.4) \begin{equation} q=G(q),\quad q(z_0)=0\quad \text{if} \quad z_0 \in \mathcal{A},\quad q(z_0)=1\quad \text{if}\ z_0 \in \mathcal{B}. \end{equation}

The two conditions directly follow from the definitions of the hitting times equation (3.1) and of the committor equation (3.2) (Prinz et al. Reference Prinz, Held, Smith and Noé2011): if the process starts in $\mathcal{A}$ (respectively $\mathcal{B}$ ), it stops at the first step, $T_{\mathcal{A}}=0$ , $T_{\mathcal{B}}\gt 0$ (respectively $T_{\mathcal{B}}=0$ , $T_{\mathcal{A}}\gt 0$ ) and $\mathbb{P}(T_{\mathcal{B}}(z_{0})\lt T_{\mathcal{A}}(z_{0}))=0$ (respectively $\mathbb{P}(T_{\mathcal{B}}(z_{0})\lt T_{\mathcal{A}}(z_{0}))=1$ ).

The proxies $\tilde {q}$ and $G^\dagger$ are interesting because it is simple and moderately expensive to estimate them from data with enough precision given an unconstrained, representative, sample of the dynamics. Note that this adds the additional assumption that our system is ergodic: sample means converge to ensemble averages and, among other conditions, there are no isolated areas in phase space. In particular, transition probabilities contained in $G$ can be estimated from unconstrained realisations of the process. Conversely, a direct sampling of the committor from data would require sampling many realisations of the dynamics (until they reach $\mathcal{A}$ or $\mathcal{B}$ ) for each value of $z_0$ (Lucente et al. Reference Lucente, Herbert and Bouchet2022a ). Using that method, the estimate is the proportion of realisations that reach $\mathcal{B}$ before $\mathcal{A}$ , and the variance of the estimated committor is at least the inverse of the number of realisations started from $z_0$ .

The final stage of approximation consists in discretising the phase space in which our process $z(t)$ evolves. In studies such as that by Lucente et al. (Reference Lucente, Rolland, Herbert and Bouchet2022b ), the main reason for this discretisation was numerical cost reduction. Indeed, the purpose of this study was to estimate the committor for a system described by an explicit, analytic stochastic differential equation. In that case, the direct and adjoint of the exact linear operator $G$ were known. In the limit $\Delta t \rightarrow 0$ , the probabilities $\tilde {q}$ and $q$ were respectively solutions of the direct Fokker–Planck equation and the steady adjoint Fokker–Planck equation (Gardiner Reference Gardiner2009; Friedrich et al. Reference Friedrich, Peinke, Sahimi and Reza Rahimi Tabar2011), which are the time-continuous equivalents of (3.3) and (3.4). However, the random variable had up to six degrees of freedom: the cost of the numerical resolution of the Fokker–Planck equation was prohibitive, while the analytical resolution was intractable. In the present case, the original system, a flow governed by the incompressible Navier–Stokes equations, does not follow a Markov chain. This Markov chain has to be created and creating this discrete process is one of the possible approaches. This is our last stage of modelling: describing our system in a probabilistic manner, assuming that it is Markovian with a time step $\Delta t$ , and finally reducing its dimension to a finite number of states.

After that stage of discretisation, our process can now visit only a finite number $N$ of discrete states. As a consequence, the committor $\boldsymbol{q}$ and the probability $\boldsymbol{p}$ are now vectors in $\mathbb{R}^N$ , and $\boldsymbol{G}$ and $\boldsymbol{G}^\dagger$ are $N\times N$ matrices. Each component $q_{\!j}$ and $p_{\!j}$ corresponds to a state $z_{\!j}$ , $1\le j \le N$ . Meanwhile, the entries $G_{nj}$ of $\boldsymbol{G}^\dagger$ are the transition probabilities from state $n$ to state $j$ . The Chapman–Kolmogorov equation (3.3) for advancing $\boldsymbol{p}$ in time then reads

(3.5) \begin{equation} p_{\!j}(t+\Delta t)=\sum _{n=1}^N G_{nj}p_{n}(t),\quad G_{nj}\ge 0,\quad \sum _{n=1}^N G_{nj}=1,\quad 1\le j \le N, \end{equation}

where the transition probability follows a normalisation condition, because the system necessarily transits somewhere. Thus, if we are able to estimate from data the transition probabilities contained in $G_{nj}$ , through one representation of the process or another, we can then obtain an estimate of the committor on the discrete states by solving (3.4). Using that estimate, it may then be possible to interpolate the committor on any state, for instance, using an approach derived from the method used to discretise the phase space.

3.2. Construction of the analogue Markov chain

In this section, we explain how the phase space is discretised in practice to construct the $N$ states on which our analogue Markov chain evolves and how the corresponding transition probabilities are estimated from DNS data. The resulting process is termed an analogue Markov chain and it provides the basis for the committor computation presented in the next section.

The method that we follow to construct our effective Markov dynamics on discrete states is based on the method of analogues proposed by Lorenz (Reference Lorenz1969b ,Reference Lorenz c ) for simple atmospheric predictability. The original approach is deterministic and is based on the observation that, up to a predictability time horizon, similar initial conditions will remain close as they evolve in time. Thus, if one has sampled a long enough time series of a process $\{ x_a(t)\}$ , the most basic formulation of the analogue method to forecast the future $x(t+\Delta t)$ of a new observation $x(t)$ is to search for the closest state, its nearest neighbour, according to a norm in phase space, among the so-called analogues $x_a(t)$ . The prediction is then $x_a(t+\Delta t)\simeq x(t+\Delta t)$ (figure 2 a). To account for the chaotic nature of the turbulent flow and to propose a process from which we can estimate a committor, we use a variant of the analogue method, termed the analogue Markov chain. In this method, the process does not deterministically evolve along the same path as its nearest neighbour among the analogues. Instead, at each time step, it is allowed to jump at random to one out of $K$ states: the jump is chosen among the $K$ time evolutions of the $K$ nearest neighbours of our state (figure 2 b). This creates a Markov chain in the state space comprising the observed analogues.

To construct this analogue Markov chain, we first need to define the state space where it takes place, which is constructed from data observed in our turbulent flow. For this matter, we sample predictors: vectors $\boldsymbol{Z}$ in the phase space $\boldsymbol{\chi }\subset \mathbb{R}^{D}$ . Each of the $D$ components of $\boldsymbol{Z}$ is given by the value of a predictor field at a given spatial point $\boldsymbol {x}$ , at a given instant $t$ . Actually, the main challenge in this paper is the choice of predictor fields in the first place. Thus, in this paper, we consider more than one spatial point $\boldsymbol {x}$ , say a number $N_x$ of such points, where the predictor fields are calculated. At each one of these points, we sample uninterrupted time series of $N_t$ samples of $\boldsymbol{Z}$ with time separation $\Delta t$ , i.e. for $t=0, \Delta t, 2\Delta t, 3\Delta t, \ldots , (N_{t}-1)\Delta t$ (where $t=0$ stands for an arbitrary initial time). Our state space is constituted by the resulting $\{\boldsymbol{Z}_{n}\}_{1\leq n \leq N_{t} N_x}$ , concatenated as $N_x$ successive time series. For all $1\le m\le N_x$ , $\boldsymbol{Z}_{l+1+(m-1) N_t}$ comprises the $D$ predictors sampled at points $\boldsymbol {x}_m$ and $\boldsymbol{Z}_{n}$ at instants $l\Delta t$ for $0\le l\le N_{t}-1$ .

We now present the construction of the transition matrix $\boldsymbol{G}^\dagger$ of the analogue Markov chain from the series of observed events $\boldsymbol{Z}_{n}\in \left \{\boldsymbol{Z}_{n}\right \}_{1\leq n \leq N_{t} N_x}$ at the $N_x$ spatial points. In our analogue Markov chain, each state $\boldsymbol{Z}_{n}$ can follow the transition of each of its $K$ nearest neighbours with equal probability (see figure 2 b, for $K=5$ ). The nearest neighbours are determined according to their closeness in Euclidean norm to $\boldsymbol{Z}_{n}$ in phase space among all $\boldsymbol{Z}_{l+1 + (m-1)N_t}$ with $1\le m\le N_x$ and $0 \le l \le N_t-1$ (note that in the nearest-neighbour search, we exclude the $\boldsymbol{Z}$ at the final instant of every time series corresponding to every point $\boldsymbol {x}_m$ so as to ensure that each of the $K$ neighbours has a recorded evolution in state space). We denote by $\mathcal{T}_{nk}$ the indices of the $K$ nearest neighbours ( $1\le k\le K$ ) of $\boldsymbol{Z}_{n}$ , which means that the $\boldsymbol{Z}_{n'}$ with $n'=\mathcal{T}_{nk}$ realise the $K$ smallest values of the Euclidean norm $\parallel \boldsymbol{Z}_{\!n}-\boldsymbol{Z}_{n'}\parallel$ . Since each $\boldsymbol{Z}_{\mathcal{T}_{nk}}$ evolves into $\boldsymbol{Z}_{\mathcal{T}_{nk}+1}$ in the recorded time series, in our analogue Markov chain, $\boldsymbol{Z}_{n}$ can transit to any of the $\boldsymbol{Z}_{\mathcal{T}_{nk}+1}$ with probability $1/K$ and does not transit to any other state. Given that rule, we can write all the entries of the $(N_tN_x)\times (N_tN_x)$ transition matrix $\unicode{x1D642}^\dagger$ as follows:

(3.6) \begin{equation} G_{nj}^\dagger =\left \{ \begin{array}{ccl} 1/K & &\text{if } \exists \, k\in [1,K] \text{ such that } j=\mathcal{T}_{nk}+1,\\[5pt] 0 & &\text{otherwise.} \end{array} \right . \end{equation}

In this article, we use this transition matrix as a basis to compute the committor in our state space, as we detail in the next section.

Figure 2.

(a) Schematic of the analogue method to predict a deterministic evolution (figure reproduced from Miloshevich et al. Reference Miloshevich, Lucente, Yiou and Bouchet2024). (b) Schematic of the analogue Markov chain method taking $K=5$ as an example. On the left-hand side, the event $\boldsymbol{Z}_{n}$ (shown in red and not corresponding to the event at the final instant of any time series) is surrounded by its $K$ nearest events $\{\boldsymbol{Z}_{\mathcal{T}_{nk}}\}_{1\leq \mathcal{T}_{nk} \leq N_{t}-1}, 1\leq k\leq K$ including itself, as indicated within the circle. On the right-hand side are their subsequent observation events and the transition to these subsequent observations is denoted by black lines. The red lines are the possible transitions in the analogue Markov chain. The schematic is adapted from figure 3 of Lucente et al. (Reference Lucente, Rolland, Herbert and Bouchet2022b ). The observations $\boldsymbol{Z}_{m}$ and $\boldsymbol{Z}_{m+1}$ are added here to represent events outside of the neighborhood of $\boldsymbol{Z}_{n}$ .

3.3. Computing the committor function from the analogue Markov chain

Now that we have presented the construction of the transition matrix $\tilde {\unicode{x1D642}}$ of our analogue Markov chain, we first explain how $\tilde {\unicode{x1D642}}$ is modified compared with its original form $\unicode{x1D642}^\dagger$ to get the committor. We then detail how we use this transition matrix $\tilde {\unicode{x1D642}}$ of analogue Markov chain to compute the committor (Prinz et al. Reference Prinz, Held, Smith and Noé2011; Lucente et al. Reference Lucente, Rolland, Herbert and Bouchet2022b ).

A series of $\boldsymbol{Z}$ that are representative enough for a committor computation by our method necessarily contain a number $N_{\mathcal{A}}\gt 0$ of instants $n$ where $\boldsymbol{Z}_{\!n}$ is in $\mathcal{A}$ and a number $N_{\mathcal{B}}\gt 0$ of instants where $\boldsymbol{Z}_{n'}$ is in $\mathcal{B}$ . In the initial Markov chain, the system can enter these states by following recorded transitions, but can also exit them in the same way. However, if we wish to compute a committor $\boldsymbol{q}$ , the dynamics must be allowed to enter either of these two sets, but not to exit them: they must be absorbing. We therefore modify the transition matrix and create an auxiliary Markov chain where the two sets are strictly absorbing. A simple way to modify our process consists in changing the rows of $\boldsymbol{G}$ corresponding to the states $\boldsymbol{Z}_{\!n} \in \mathcal{A}\cup \mathcal{B}$ so that each of these states transits only to itself. We thus create a new transition matrix $\tilde {\boldsymbol{G}}$ such that

(3.7) \begin{align} \notag &\tilde {G}_{nj}=G_{nj},\quad 1\le j \le N_{x} N_t \quad \text{if} \quad \boldsymbol{Z}_{\!n} \notin \mathcal{A}\cup \mathcal{B},\\ \tilde {G}_{nn}=1\quad \text{and}\quad \tilde {G}_{nj} &=0 ,\quad 1 \le j \le N_{x} N_t,\, j\ne n \quad \text{if} \quad \boldsymbol{Z}_{\!n} \in \mathcal{A}\cup \mathcal{B}. \end{align}

Note that this modification respects the normalisation condition of (3.5).

For our Markov chain, we can define the committor as a vector $\tilde {\boldsymbol{q}}\in \mathbb{R}^{N_{x} N_t}$ , each of its components $q_n$ being the committor value associated with state $\boldsymbol{Z}_{\!n}$ . For our finite-dimensional analogue Markov chain process, a time-invariant committor satisfying the backward Chapman–Kolmogorov equation (3.5) then reads

(3.8) \begin{equation} \sum _{j=1}^{N_tN_x}\tilde {G}_{\textit{ij}}q_{\!j}=q_i,\quad 1\le i\le N_xN_t,\quad q_n=0\,\quad \text{if} \quad \boldsymbol{Z}_{\!n} \in {\mathcal{A}},\quad q_{n'}=1\quad \text{if} \quad \boldsymbol{Z}_{n'} \in {\mathcal{B}}. \end{equation}

Following Prinz et al. (Reference Prinz, Held, Smith and Noé2011), Lucente et al. (Reference Lucente, Rolland, Herbert and Bouchet2022b ), we solve (3.8) in two steps: first, by computing the eigenvectors of $\tilde {G}$ for its largest eigenvalue and second, by expressing $q$ as a linear combination of these eigenvectors. Perron–Frobenius theorems indicate that the matrix $G$ , which is non-negative (zero or positive entries) and irreducible (it conserves probability), see (3.5), has a single eigenvalue equal to $1$ , its spectral radius. The corresponding eigenspace therefore has dimension one. The same properties are found for $G^\dagger$ . Because of the modifications applied to turn $G$ into $\tilde {G}$ (turning $\mathcal{A}$ and $\mathcal{B}$ into absorbing states), the spectral radius of $\tilde {G}$ remains $1$ ; however, the eigenvalue $1$ is no longer simple: its corresponding eigenspace has dimension $N_{\mathcal{A}}+N_{\mathcal{B}}$ (Prinz et al. Reference Prinz, Held, Smith and Noé2011; Lucente et al. Reference Lucente, Rolland, Herbert and Bouchet2022b ). Since $\tilde {G}q=q$ , the committor vector lies in that eigenspace. If we term $\{\boldsymbol{v}_{\!m}\}_{1\le m\le N_{\mathcal{A}}+N_{\mathcal{B}}}$ a basis of this eigenspace, there exists a unique collection of $M$ scalars $\alpha _m$ such that

(3.9) \begin{equation} \boldsymbol{q}=\sum _{m=1}^M \alpha _m\boldsymbol{v}_{\!m}. \end{equation}

These scalars are in turn the solution of the system of equations given by the $N_{\mathcal{A}}$ conditions

(3.10) \begin{equation} 0=q_n=\sum _{m=1}^M \alpha _m v_{n,m} \quad \text{for } \quad \boldsymbol{Z}_{\!n} \in \mathcal{A}, \end{equation}

and the $N_{\mathcal{B}}$ conditions

(3.11) \begin{equation} 1=q_{n'}=\sum _{m=1}^M \alpha _m v_{n',m} \quad \text{for } \quad \boldsymbol{Z}_{n'} \in \mathcal{B}, \end{equation}

where $v_{n,m}$ is the component $n$ of eigenvector $\boldsymbol{v}_{\!m}$ . Thus, in our python implementation, we first compute a number of eigenvectors for the largest eigenvalues, strictly larger than $M$ , so as to be able to span the eigenspace of $\tilde {\boldsymbol{G}}$ for eigenvalue $1$ . From the eigenvectors associated with eigenvalue $1$ , we construct the system of $N_{\mathcal{A}}+N_{\mathcal{B}}$ equations ((3.10) and (3.11)) and solve it. Using the computed scalars, we then construct our committor vector equation (3.9) defined on the states of the analogue Markov chain.

Note that, in principle, all states in $\mathcal{A}$ and all states in $\mathcal{B}$ could be further aggregated into two internally mixed macrostates. In that case, a reduced transition matrix compared with (3.7) could be constructed. The multiplicity associated with the eigenvalue $1$ would then decrease from $N_{\mathcal{A}}+N_{\mathcal{B}}$ to $2$ , which would lower the dimension of the corresponding linear system and simplify the computation. Since the committor only depends on which set is reached first, such an aggregation would preserve the same committor. In the present work, however, we retain the original formulation to remain consistent with the current implementation of the transition matrix and the associated event indexing.

Note that this result also provides a method to estimate the value of the committor for any new state $\boldsymbol{Y}\in \chi$ sampled in the flow. According to the principle of the representation of the dynamics by an auxiliary analogue Markov chain, we can start a new realisation of the chain starting from $\boldsymbol{Y}$ . It can transit towards the future step of any of its $K$ nearest neighbours among the $N_{t} N_x$ sampled analogues $\boldsymbol{Z}_{\!n}$ , drawn uniformly. The chain then proceeds until it is absorbed by either $\mathcal{A}$ or $\mathcal{B}$ . Thus, the probability that a trajectory starting from $\boldsymbol{Y}$ reaches $\mathcal{B}$ before $\mathcal{A}$ is given by the arithmetic average of the committor values over the $K$ nearest neighbours of $\boldsymbol{Y}$ . Here, the out-of-sample event $\boldsymbol{Y}$ is an isolated query event, so transitions from in-sample events to $\boldsymbol{Y}$ are not admissible, and the original transition structure and committor values remain unchanged. If we denote the indices of these $K$ states by $\{\mathcal{T}_{k}(\boldsymbol{Y})\}_{1 \leq k \leq K}$ , an estimate for the committor on $\boldsymbol{Y}$ is thus given by the arithmetic average

(3.12) \begin{equation} q(\boldsymbol{Y}) = \frac {1}{K} \sum _{k=1}^{K} q_{\mathcal{T}_{k}(\boldsymbol{Y})+1}. \end{equation}

In the present work, we set $K=100$ , but we also tried $K=50$ and $K=150$ without qualitative change of conclusions. The set $\{\boldsymbol{Z}_{n}\}_{1\leq n \leq N_{x}N_{t}}$ is sampled from the DNS data described in § 4. The sampling details are presented in § 4.3, where the definitions of the sets $\mathcal{B}$ and $\mathcal{A}$ are also provided.

3.4. Evaluation of the quality of an approximate committor function

Following Lucente et al. (Reference Lucente, Rolland, Herbert and Bouchet2022b ), we use the Brier score to evaluate the precision of the obtained approximate committor function. It was originally proposed by Brier (Reference Brier1950) to estimate the precision of weather forecasts. First, we define $q_{B}$ as a binary verification function such that $q_{B}(\boldsymbol{Y})=1$ if $T_{\mathcal{B}}(\boldsymbol{Y})\lt T_{\mathcal{A}}(\boldsymbol{Y})$ , and $q_{B}(\boldsymbol{Y})=0$ if $T_{\mathcal{B}}(\boldsymbol{Y})\gt T_{\mathcal{A}}(\boldsymbol{Y})$ . The accurate committor function is represented by $\hat {q}$ . The outcomes of $q_{B}$ have a Bernoulli distribution, where $\mathbb{P}(q_{B}(\boldsymbol{Y}) = 1) = \hat {q}(\boldsymbol{Y})$ and $\mathbb{P}(q_{B}(\boldsymbol{Y}) = 0) = 1 - \hat {q}(\boldsymbol{Y})$ . The Brier score is defined as

(3.13) \begin{equation} \textit{BT}_{N}=\frac {1}{N}\sum _{i=1}^{N}\left (q(\boldsymbol{Y}_{i})-q_{B}(\boldsymbol{Y}_{i})\right )^{2}, \end{equation}

where $N$ is the number of events used for testing the performance of the approximate committor function. The Brier score therefore takes values between $0$ and $1$ . Assuming that the sampled process in the set $\{\boldsymbol{Y}_{i}\}$ is ergodic in phase space, we have the equivalence between the arithmetic mean and the ensemble mean, i.e.

(3.14) \begin{equation} \lim _{N\to \infty }\textit{BT}_{N}=\mathbb{E}(\textit{BT}_{N})=\left \Vert \sqrt {\hat {q}(1-\hat {q})}\right \Vert _{\rho }^{2}+\left \Vert \hat {q}-q\right \Vert _{\rho }^{2}, \end{equation}

where $\Vert f\Vert _{\rho }^{2}=\int _{\chi }f^{2}(\boldsymbol{Y})\rho (\boldsymbol{Y})\,\mathrm{d}\boldsymbol{Y}$ is the norm weighted according to the stationary distribution function $\rho (\boldsymbol{Y})$ . The first term in (3.14) depends only on the nature of the system and is independent of the quality of the approximate committor function; it is a fixed lower bound resulting from the probabilistic nature of the prediction. The second term in (3.14) represents the quadratic deviation of the approximate committor function from the accurate one. The closer $q$ is to $\hat {q}$ , the lower is the Brier score.

In the present work, this point is particularly important because each choice of predictor defines a different coarse-grained phase space. Therefore, the Brier score does not only evaluate the quality of the approximate committor itself, but also the amount of relevant information retained by the chosen predictor. A lower Brier score means that the selected predictor, or equivalently the reduced phase-space representation built from it, is more informative about the future occurrence of extreme uncertainty-production events. In this sense, comparing the Brier scores obtained from different predictors provides a quantitative way to assess how closely the corresponding coarse-grained phase spaces capture the dynamically relevant structure of the original Navier–Stokes phase space.

4.1. DNS numerical set-up

To study the growth of average uncertainty energy in three-dimensional periodic/homogeneous turbulence, we use a fully de-aliased pseudo-spectral code to perform DNS of forced incompressible Navier–Stokes turbulence in a periodic box of size $\mathcal{L}=2\pi$ in all three normal directions. Time advancement is achieved with a second-order Runge–Kutta scheme. We use the same adaptive time step $\textrm {d}t$ for both $\boldsymbol{u}^{(1)}$ and $\boldsymbol{u}^{(2)}$ as follows. The Courant–Friedrichs–Lewy (CFL) condition, with CFL number equal to $0.4$ , gives a time step $\textrm {d}t^{(1)}$ when applied to $\boldsymbol{u}^{(1)}$ and $\textrm {d}t^{(2)}$ when applied to $\boldsymbol{u}^{(2)}$ . We set $\textrm {d} t=\min (\textrm {d}t^{(1)},\textrm { d}t^{(2)})$ . The code strategy is detailed by Vincent & Meneguzzi (Reference Vincent and Meneguzzi1991) and has been used in our previous study (Ge et al. Reference Ge, Rolland and Vassilicos2023, Reference Ge, Rolland and Vassilicos2025a , Reference Ge, Rolland and Vassilicosb ). We first generate a statistically stationary reference flow. To generate the reference flow, we use a von Kármán initial energy spectrum with the same coefficients as Yoffe (Reference Yoffe2012) and random initial Fourier phases. We integrate the reference flow until it reaches a statistically stationary state and then seed it with small perturbations to create the perturbed flow at a time which we refer to as $t_{0}=0$ . The reference flow and the perturbed flow are simulated in parallel using the same time step. The key parameters of the reference flow are given in table 1.

Table 1.

Parameters of the reference flows, where all the parameter here are presented after averaged in time when the flows are statistically stationary. The large eddy turnover time $T_{0}^{(1)}$ is $T_{0}\equiv L/U$ for the reference flow where $U=\sqrt {2\langle E\rangle /3}$ is the root mean square velocity and $L=(3\pi /4\langle E\rangle )\int\!k^{-1}\hat {E}(k)\,\textrm {d}k$ is the integral length scale. $N$ is the resolution of the simulations, $\nu$ is the kinematic viscosity, $\varepsilon$ is the turbulence dissipation rate. ${Re}= UL/\nu$ is the Reynolds number based on $L$ and ${Re}_{\lambda }\equiv U\lambda /\nu$ is the Reynolds number based on the Taylor length $\lambda \equiv \sqrt {10\langle E\rangle \nu /\varepsilon }$ . $k_{max }=N/3$ is the maximum resolvable wavenumber and $\eta =(\nu ^{3}/\varepsilon )^{1/4}$ is the Kolmogorov length.

Table 1. Long description

A table comparing parameters of reference flows in turbulence simulations. The table has 4 rows and 10 columns. Column headers are Case, N cubed, nu, epsilon, U, L, T 0, Re, Re lambda, and k max eta. Row 1: Case, F 1 128; N cubed, 128 cubed; nu, 0.0378; epsilon, 0.0999; U, 0.5153; L, 1.6806; T 0, 3.2611; Re, 22.9487; Re lambda, 16.7970; k max eta, 6.41. Row 2: Case, F 1 512; N cubed, 512 cubed; nu, 0.0010; epsilon, 0.0971; U, 0.6038; L, 1.066; T 0, 1.7646; Re, 643.6905; Re lambda, 143.7088; k max eta, 1.71. Row 3: Case, F 2 128; N cubed, 128 cubed; nu, 0.0300; epsilon, 0.1466; U, 0.4748; L, 1.087; T 0, 2.2924; Re, 17.1952; Re lambda, 13.2071; k max eta, 4.90. Row 4: Case, F 2 512; N cubed, 512 cubed; nu, 0.0015; epsilon, 0.2171; U, 0.6501; L, 0.6473; T 0, 0.9976; Re, 280.2203; Re lambda, 90.7673; k max eta, 1.90.

The generation of the perturbed flow is detailed in our previous work (Ge et al. Reference Ge, Rolland and Vassilicos2023), where

(4.1) \begin{equation} \hat {\boldsymbol{u}}^{(2)}(\boldsymbol{k},t_0)=\left \{ \begin{array}{ccc} \hat {\boldsymbol{u}}^{(1)}(\boldsymbol{k},t_0)& &\text{if } \left |\boldsymbol{k}\right |\lt k_{0}, \\[5pt] \text{randomly generated}& &\text{if } \left |\boldsymbol{k}\right |\geq k_{0}. \end{array}\right . \end{equation}

We denote respectively by $\hat {\boldsymbol{u}}^{(1)}$ and $\hat {\boldsymbol{u}}^{(2)}$ the velocity fields of the reference flow and the perturbed flow in Fourier space. In the present work, we set $k_{0}\eta =1.5$ , the same as in our previous work (Ge et al. Reference Ge, Rolland and Vassilicos2023). We use this perturbation strategy to keep the perturbed flow as close as possible to the reference flow. Specifically, the perturbed field is constructed so that the one-dimensional energy spectrum is preserved at the initial time $\hat {E}^{(1)}(k,t_{0})=\hat {E}^{(2)}(k,t_{0})$ for all resolved wavenumbers. In our previous works (Ge et al. Reference Ge, Rolland and Vassilicos2025a ), we have shown that this way of constructing the perturbation does not affect the subsequent uncertainty-growth rate when compared with the more standard approach of directly adding a white-noise perturbation. We have also shown that the uncertainty-growth dynamics is insensitive to both the amplitude of the initial perturbation and the range of perturbed wavenumbers (Ge et al. Reference Ge, Rolland and Vassilicos2023).

In the present work, we consider two different types of external forcing. For each forcing type, we test two Reynolds-number configurations. The first one corresponds to a Reynolds number chosen to be as low as possible within the present framework. This allows us to maintain a comparatively high spatial resolution and, at the same time, to obtain a relatively long chaotic exponential-growth regime of the uncertainty energy. Such a long growth regime is advantageous for the present statistical analysis because it provides long continuous time series and therefore a large number of temporal samples $N_t$ . The second configuration corresponds to a Reynolds number chosen to be as high as possible within our computational limits. The purpose of these higher-Reynolds-number cases is to reduce the direct influence of the forcing on the smaller scales and thereby to better isolate the inertial effects of interest. Because increasing the Reynolds number also shortens the chaotic exponential-growth regime, these higher-Reynolds-number simulations do not provide time series as long as those of the lower-Reynolds-number cases. To compensate for this reduced temporal extent and still obtain sufficient sampling, we increase the number of observation windows $N_x$ . The corresponding sampling details are discussed in § 4.3.

One type of forcing, referred to as case F1, is a negative-damping force applied exclusively at low wavenumbers, a method commonly used in previous studies (Boffetta & Musacchio Reference Boffetta and Musacchio2017; Berera & Ho Reference Berera and Ho2018; Ge et al. Reference Ge, Rolland and Vassilicos2023). Here, ‘negative damping’ refers to the fact that the forcing is aligned with the large-scale velocity modes by injecting energy into them, whereas a damping force is anti-aligned (opposed) to the velocity (Linkmann & Morozov Reference Linkmann and Morozov2015; McComb et al. Reference McComb, Linkmann, Berera, Yoffe and Jankauskas2015). Therefore, the forcing is expressed as

(4.2) \begin{equation} \hat {\boldsymbol{f}}^{(2)}\left (\boldsymbol{k},t\right )=\hat {\boldsymbol{f}}^{(1)}\left (\boldsymbol{k},t\right )=\left \{ \begin{array}{ccc} \frac {\varepsilon _{0}}{2E_{\!f}^{(1)}}\hat {\boldsymbol{u}}^{(1)}\left (\boldsymbol{k},t\right )& &\text{if } 0\lt \left |\boldsymbol{k}\right |\leq k_{\!f},\\[7pt] 0& &\text{otherwise,} \end{array}\right . \end{equation}

where $\varepsilon _{0}=0.1$ and $k_{\!f}=2.5$ .

Figure 3.

Time evolution of (a) the average uncertainty energy in a semilogarithmic plot and (b) the corresponding growth exponent in a linear plot. $\gamma \equiv ({1}/{2})\,\mathrm{d}\ln \langle E_{\varDelta }\rangle /\mathrm{d}t$ , $\tau _{\eta }$ denotes the Kolmogorov time scale and $\langle E_{tot}\rangle =\langle |\boldsymbol{u}^{(1)}|^{2}/2+|\boldsymbol{u}^{(2)}|^{2}/2 \rangle$ is the total mean kinetic energy.

The other forcing type, referred to as case F2, is a single-mode divergence-free force that remains constant in time:

(4.3) \begin{equation} \boldsymbol{f}^{(2)}\left (\boldsymbol{x},t\right )=\boldsymbol{f}^{(1)}\left (\boldsymbol{x},t\right )=f_{0}\left ( \begin{array}{ccc} \cos \left (k_{0}y\right )\sin \left (k_{0}z\right )\\[5pt] \cos \left (k_{0}z\right )\sin \left (k_{0}x\right )\\[5pt] \cos \left (k_{0}x\right )\sin \left (k_{0}y\right )\\ \end{array}\right )\!, \end{equation}

where $f_{0}=0.5$ and $k_{0}=2$ .

4.2. DNS results

Figure 3 presents the time evolution of $\langle E_{\varDelta } \rangle$ , where $\tau _{\eta }\equiv (\nu /\varepsilon )^{1/2}$ denotes the Kolmogorov time scale. It is observed that the average uncertainty energy $\langle E_{\varDelta } \rangle$ enters an exponential-growth regime after a sharp initial decrease, as already noted by Ge et al. (Reference Ge, Rolland and Vassilicos2023). As mentioned in § 4.1, the lower-Reynolds-number cases exhibit a much longer exponential-growth interval. To characterise this interval more clearly, figure 3(b) also shows the time evolution of the instantaneous growth exponent, defined as $\gamma \equiv ({1}/{2})\,\mathrm{d}\ln \langle E_{\varDelta }\rangle /\mathrm{d}t,$ whose nearly constant value indicates exponential growth. In the present work, we focus exclusively on the regime of exponential growth of $\langle E_{\varDelta } \rangle$ , namely the plateau interval identified in figure 3(b). The time-sampling intervals used in the following analysis are all chosen within this range, as listed in table 2. Figure 4 shows the p.d.f. of the uncertainty production at different times during this regime. The p.d.f.s of $P_{\varDelta }$ at different times during this time range approximately collapse if normalised by the p.d.f.’s maximum value and standard deviation. The average uncertainty production rate normalised by its standard deviation is therefore an approximately statistically stationary process as figure 5 in fact demonstrates. It is therefore reasonable to disregard the chronological order of the observations, after normalising them by their standard deviation, when constructing the analogue Markov chain.

Table 2.

Parameters of the data sampling. The interval $[t_{min }/\tau _{\eta },\, t_{max }/\tau _{\eta }]$ gives the starting and ending times of the sampling, where $\tau _{\eta }$ is the Kolmogorov time scale. $N_{{case}}$ is the number of DNS realisations performed from the same reference flow with independently generated random perturbations. $N_{x}$ is the number of observation boxes used in the training dataset and $N_{t}$ is the length of the time sequence sampled in each observation box. $2L_{\mathbb{D}}$ is the size of the observation box and $L_{\varDelta }$ is the uncertainty-energy length scale. $L_{P_{\varDelta }}$ , $L_{\unicode{x1D64E}}$ , $L_{\boldsymbol{\varOmega }}$ and $L_{[\unicode{x1D64E},\boldsymbol{\varOmega }]}$ are the correlation lengths, defined from the integrals of the corresponding auto-correlation functions, of $P_{\varDelta }$ , $\Vert \unicode{x1D64E}\Vert _{F}$ , $\Vert \boldsymbol{\varOmega }\Vert _{F}$ and $\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }]\Vert _{F}$ , respectively. $\Delta t$ is the sampling time interval. $\mathbb{P}(\mathcal{A})$ and $\mathbb{P}(\mathcal{B})$ are the probabilities of the event sets $\mathcal{A}$ and $\mathcal{B}$ , respectively.

Table 2. Long description

A table with 4 rows and 17 columns comparing parameters of data sampling across different cases. The columns are labeled as follows: Case, t_min/η, t_max/η, N_case, N_x, N_1, 2L_D, L_Δ, L_P, L_S, L_Ω, L|S,Q|, Δt/η, P(A), P(B). The rows present data for cases F1_256, F1_512, F2_256, and F2_512. Each row lists specific values for each parameter, providing a detailed comparison of the sampling intervals and probabilities of event sets across different cases.

Figure 4.

Time evolution of p.d.f.s of $P_{\varDelta }$ for cases (a) F1 $_{128}$ , (b) F1 $_{512}$ and (c) F2 $_{128}$ and (d) F2 $_{512}$ . Similarly to Ge et al. (Reference Ge, Rolland and Vassilicos2023), p.d.f.s are normalised by their maximum $\mathcal{P}_{max }$ and are plotted versus $P_{\varDelta }/\sigma _{P_{\varDelta }}$ , where $\sigma _{P_{\varDelta }}$ is the standard deviation of $P_{\varDelta }$ with $\sigma ^{2}_{P_{\varDelta }}=\int _{P_{\Delta {min}}}^{P_{\Delta {max}}}(P_{\varDelta }-\langle P_{\varDelta }\rangle )^{2}\mathcal{P}(P_{\varDelta })\,\textrm {d}P_{\varDelta }$ .

Figure 5.

Time evolution of the average uncertainty production, normalised by its standard deviation. The top blue $x$ -axis corresponds to cases F1 $_{128}$ and F2 $_{128}$ , whereas the bottom red $x$ -axis corresponds to cases F1 $_{512}$ and F2 $_{512}$ .

4.3. Statistical analysis set-up

4.3.1. Data sampling

The data used to establish the analogue Markov chain are sampled from the DNS using a number of observation boxes. These boxes, denoted as $\mathbb{D}$ , are small boxes

(4.4) \begin{equation} \mathbb{D} = \big \{ \big ( x_1, x_2, x_3 \big ) \mid \big | x_i - x_i^{\mathbb{D}} \big | \lt L_{\mathbb{D}}, i = 1, 2, 3 \big \}, \end{equation}

where $\boldsymbol{x}^{\mathbb{D}}$ represents the geometrical centre of the sampling box, which is randomly determined within the simulation domain.

The size of the sampling box, $2L_{\mathbb{D}}$ , is chosen to be close to the average distance between two velocity stagnation points, which is equal to $\mathbb{P} (\boldsymbol{u}=0 )\lambda$ (Goto & Vassilicos Reference Goto and Vassilicos2009), where $\lambda$ is the Taylor microscale and $\mathbb{P} (\boldsymbol{u}=0 )$ is the probability density for the instantaneous fluid velocity to be zero in the reference flow. This characteristic length scale may be interpreted as being close to the average size of a turbulent eddy. In practice, the half-size of each observation box, denoted by $L_{\mathbb{D}}$ , is chosen in the form $L_{\mathbb{D}}=2^{n}\Delta x$ , where $\Delta x=2\pi /N$ is the DNS spatial resolution. The integer exponent $n$ is selected such that the full box size $2L_{\mathbb{D}}=2^{n+1}\Delta x$ is the smallest grid-compatible length scale larger than $\mathbb{P} (\boldsymbol{u}=0 )\lambda$ . This choice ensures that the observation boxes can exactly tile the computational domain.

To reduce spatial correlations between different observation boxes, the distance between their geometric centres is required to be larger than all relevant correlation lengths, including those of the candidate predictors (see § 4.3.3) and of the uncertainty-production term. Here, each correlation length is denoted as $L_{P_{\varDelta }}$ , $L_{\unicode{x1D64E}}$ , $L_{\boldsymbol{\varOmega }}$ and $L_{[\unicode{x1D64E},\boldsymbol{\varOmega }]}$ , defined as the integral of the corresponding auto-correlation function (Batchelor Reference Batchelor1953). Therefore, the distance between any two geometric centres is larger than $L_{{sep}}=\max \{L_{P_{\varDelta }}, L_{\unicode{x1D64E}}, L_{\boldsymbol{\varOmega }}, L_{[\unicode{x1D64E},\boldsymbol{\varOmega }]},L_{\varDelta }\}$ , where $L_{\varDelta }$ is the uncertainty length scale, defined as the integral length scale according to the uncertainty energy spectrum, see Ge et al. (Reference Ge, Rolland and Vassilicos2023) for details. Although $L_{{sep}}$ is chosen as the maximum among the relevant correlation lengths, it remains small compared with the size of the periodic domain, $\mathcal{L}=2\pi$ . From table 2, one obtains $L_{{sep}}/\mathcal{L}\in [0.02,0.2]$ for all the cases. Thus, even the largest separation length used to decorrelate the observation boxes is much smaller than the domain size. Under this constraint, different observation boxes do not overlap. To further reduce spatiotemporal correlations between samples collected from different boxes, the sampling interval in time is required to satisfy $\Delta t = L_{{sep}}/U_{max }$ , where $U_{max }$ is the maximum flow velocity.

The time sampling is restricted strictly to the exponential-growth regime of the uncertainty energy, namely the plateau region identified in figure 3. Together, this temporal restriction and the minimum sampling interval determine the number of temporal samples $N_t$ . The number of sampled observation boxes, denoted by $N_x$ , is then adjusted accordingly so that the total number of samples satisfies $N_xN_t \geq (\pi /L_{{sep}} )^3.$ All of these statistical parameters are listed in table 2.

Rather than using all the spatial points within each observation box to estimate the committor function, we use only the box-averaged observed events, namely $\langle \overline {\boldsymbol{Z}}_n \rangle _{\mathbb{D}}$ , where $\langle \boldsymbol{\cdot } \rangle _{\mathbb{D}}$ denotes the spatial average over the observation box and the overbar indicates normalisation by the standard deviation of the corresponding global field. As shown in § 4.2, this normalisation makes the statistics of these quantities approximately statistically steady, which is advantageous for the present data-driven construction. The normalisation here is intended to help appreciate how the input data deviate from their mean. The key parameters of the data sampling are gathered in table 2. Figure 6(a) illustrates the spatial distribution of observation boxes within the simulation domain in the case $N_{x}=28$ .

Figure 6.

(a) One example of spatial distribution of observation boxes ( $N_{x}=28$ ) within the simulation domain. Inside the boxes, the streamlines of the velocity of the reference flow are plotted. (b) Spatial distribution of all potential testing data boxes.

6.1. Prediction quality of the committor functions

A primary question of the statistical analysis is to assess which predictors, among various physical fields, possess the best predictive capability, or equivalently, which choice of predictor defines a coarse-grained phase space that best captures the dynamically relevant structure of the full Navier–Stokes phase space. For this matter, we follow an evaluation procedure used by Miloshevich et al. (Reference Miloshevich, Cozian, Abry, Borgnat and Bouchet2023, Reference Miloshevich, Lucente, Yiou and Bouchet2024) based on a comparison of committors trained on different predictors or combinations of predictors albeit with a different score. Note also that we are not concerned with the delay before the extreme event, we just monitor the improvement of the prediction as the training dataset size is increased. Indeed, since the statistical results depend on both the sample size and the sampling location, we train the committor functions on datasets of different sizes by varying the number of observation boxes $N_{x}$ included in the training set. For each case, the committor is trained using $\alpha N_{x}$ observation boxes, where $\alpha$ varies from $0$ to $1$ . As mentioned in § 4.3.4, each case is simulated $N_{{case}}$ times from the same reference field, but with different randomly generated perturbations. This yields $N_{{case}}$ Brier scores corresponding to committor functions trained from different realisations. The average values and standard deviations of these scores are represented in the figures as data points and shadow regions, respectively.

For both cases F1 and F2, it is observed that, all curves in figures 8 and 9 (except the curves with all four predictors in case F1 $_{128}$ ) reach a minimum that is approximately constant with $\alpha N_x$ as $\alpha$ is increased. Note also that the shadow regions (standard deviation) become progressively narrower with increasing $\alpha$ . The more or less sustained minimum suggests that the performance of the prediction given by the committor function is no longer sensitive neither to the amount of input sampling data nor their positions as $\alpha N_x$ increases. Meanwhile, the observed decrease of the error bars with increasing window size suggests that the influence of window position and inter-window correlation is progressively reduced. Equation (3.14) implies the existence of a non-zero minimum value of the Brier score and it is probably reached by the aforementioned trained committor functions. Among all the single predictors, for all cases, $E_{\varDelta }$ exhibits the best predictive capability with the lowest stable Brier score. This optimal Brier score is attributed to $E_{\varDelta }$ containing information from both the reference flow field and the perturbed flow field. However, in practice, one cannot simultaneously obtain information from both flow fields. Therefore, it is the predictive capability of physical quantities from only a single flow field that is most worthy of attention. In the low-Reynolds-number cases, we observe that the strain rate $\Vert \unicode{x1D64E} \Vert _{F}$ , although worse than uncertainty energy $E_{\varDelta }$ , returns the best predictions among all the predictors containing information from only one field in cases F1 and F2, while the vorticity $\Vert \boldsymbol{\varOmega } \Vert _{F}$ returns the worst predictions in both cases, with the highest Brier scores in both cases F1 and F2. However, for the higher-Reynolds-number cases, the differences among $\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert _{F}$ , $\Vert \boldsymbol{\varOmega } \Vert _{F}$ and $\Vert \unicode{x1D64E} \Vert _{F}$ become much less marked. In particular, in F2 $_{512}$ , as shown in figure 9(b), $\Vert \boldsymbol{\varOmega } \Vert _{F}$ gives the lowest Brier score among the three predictors. This trend is consistent with the Reynolds-number dependence of the local velocity-gradient structure and dynamics (Das & Girimaji Reference Das and Girimaji2019), and with the fact that the small-scale structure becomes more fully developed as the Reynolds number increases (Elsinga et al. Reference Elsinga, Ishihara, Goto, Boersma and Hunt2023). In that situation, the predictive information carried by $\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert _{F}$ , $\Vert \boldsymbol{\varOmega } \Vert _{F}$ and $\Vert \unicode{x1D64E} \Vert _{F}$ can become more similar. At the same time, the $(Q,R)$ -conditioned statistics show that strong uncertainty-production events extend over a larger fraction of the vortex-stretching quadrant at higher Reynolds numbers, as shown in figure 7, which makes the vorticity-based predictor comparatively more effective. In addition, the relative importance of the strain field depends not only on its magnitude but also on its alignment with the local structures associated with uncertainty-energy production (Ge et al. Reference Ge, Rolland and Vassilicos2023). At higher Reynolds numbers, where strong production events occupy a broader region of the $(Q,R)$ plane, such geometric alignment effects may become more relevant, further reducing the distinction among $\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }] \Vert _{F}$ , $\Vert \boldsymbol{\varOmega } \Vert _{F}$ and $\Vert \unicode{x1D64E} \Vert _{F}$ as predictors.

Figure 8.

The Brier score $\textit{BT}_{N}$ of the approximated committor function versus the number of observation boxes $\alpha N_{x}$ for (a,c,e) case F1 $_{128}$ and (b,d, f) case F1 $_{512}$ , obtained from (a,b) a single scalar field, (c,d) two scalar fields and (e, f) more than two scalar fields. For comparison, the Brier score of the approximated committor function obtained from the uncertainty energy is included in all plots. The data points and the shadowed zones around them show the average value and the standard deviation of the Brier scores obtained from different sampled data sets.

Figure 9.

The Brier score $\textit{BT}_{N}$ of the approximated committor function versus the number of observation boxes $\alpha N_{x}$ for (a,c,e) case F2 $_{128}$ and (b,d, f) case F2 $_{512}$ , obtained from (a,b) a single scalar field, (c,d) two scalar fields and (e, f) more than two scalar fields. For comparison, the Brier score of the approximated committor function obtained from the uncertainty energy is included in all plots. The data points and the shadowed zones around them show the average value and the standard deviation of the Brier scores obtained from different sampled data sets.

Compared with predictions based solely on the uncertainty energy $E_{\varDelta }$ , in the lower-Reynolds-number cases, a modest improvement is observed for the combinations $\{E_{\varDelta }, \Vert [\unicode{x1D64E},\boldsymbol{\varOmega }] \Vert _{F} \}$ and $\{E_{\varDelta }, \Vert \boldsymbol{\varOmega } \Vert _{F}\}$ , whereas $\{E_{\varDelta },\Vert \unicode{x1D64E} \Vert _{F}\}$ shows a clearly better predictive performance in both cases F1 $_{128}$ and F2 $_{128}$ . However, in the higher-Reynolds-number cases, the improvement becomes much less remarkable, as shown in figures 8(d) and 9(d). Although $\{E_{\varDelta }, \Vert \unicode{x1D64E} \Vert _{F}\}$ still gives the best performance among the three combinations, the difference is much smaller and all combinations involving $E_{\varDelta }$ even perform worse than the predictor based on $E_{\varDelta }$ alone. For predictors formed of combinations of vorticity $\Vert \boldsymbol{\varOmega } \Vert _{F}$ , vortex-stretching $\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }]\Vert _{F}$ and strain rate $\Vert \unicode{x1D64E} \Vert _{F}$ , the Brier scores are larger than the scores of the uncertainty energy $E_{\varDelta }$ , as shown in figures 8(e), 8( f), 9(e) and 9( f). For both lower-Reynolds-number cases, the lowest Brier score is obtained for the combination $\{ \Vert \boldsymbol{\varOmega } \Vert _{F}, \Vert \unicode{x1D64E} \Vert _{F}\}$ , as shown in figures 8(c) and 9(c), whereas the combination $\{ \Vert \boldsymbol{\varOmega } \Vert _{F}, \Vert [\unicode{x1D64E},\boldsymbol{\varOmega }] \Vert _{F}\}$ yields the highest Brier score, i.e. the worst predictive performance. Notably, as in the single-predictor case, the predictive performances of the combinations built from $\{\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }] \Vert _{F}, \Vert \boldsymbol{\varOmega } \Vert _{F}, \Vert \unicode{x1D64E} \Vert _{F}\}$ show almost no difference from one another in the higher-Reynolds-number cases, as shown in figures 8(d) and 9(d). Even the full three-component combination $\{\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }] \Vert _{F}, \Vert \boldsymbol{\varOmega } \Vert _{F}, \Vert \unicode{x1D64E} \Vert _{F}\}$ does not lead to any clear improvement, as shown in figures 8( f) and 9( f). A similar contrast is observed when $E_{\varDelta }$ is included: in the lower-Reynolds-number cases, combining $E_{\varDelta }$ with other predictors leads to an improvement, whereas in the higher-Reynolds-number cases, such combinations provide no clear benefit and may even degrade the predictive performance relative to using $E_{\varDelta }$ alone, as the committor function estimated using all candidate scalars $\{E_{\varDelta },\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert _{F}, \Vert \boldsymbol{\varOmega } \Vert _{F}, \Vert \unicode{x1D64E} \Vert _{F}\}$ achieves the lowest Brier score in the low-Reynolds-number cases, whereas in the higher-Reynolds-number cases, it gets higher Brier score than obtained with only $E_{\varDelta }$ .

Figure 10.

Committor functions estimated from different single-field predictors for cases (a) F1 $_{128}$ , (b) F1 $_{512}$ , (c) F2 $_{128}$ and (d) F2 $_{512}$ using the analogue Markov chain, namely $\{\overline {E_{\varDelta }},\overline {\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert }_{F},\overline {\Vert \unicode{x1D64E} \Vert }_{F}, \overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}\}$ . The lines show the averages over the $N_{{case}}$ realisations and the shaded regions indicate the corresponding standard deviations. The normalisation is defined as $\overline {F}\equiv F/\sigma _{F}$ . The markers represent the direct empirical probabilities computed from the training dataset.

6.2. Dependence of committor functions on predictor values

Figure 10 presents the probabilities for evolving into an extreme positive production event before a negative production event based on committor function obtained from single normalised field predictors. Note that this normalisation neither improves nor degrades the estimated committor. However, it helps compare how different committors depend on different predictor values. Indeed, a normalised variable close to zero indicates a state close to average, while an absolute value of one or several units indicates a state that is one or several standard deviations from average.

It can be seen that the committor functions obtained from the same predictor in different flow fields are qualitatively similar. For all single predictors in $\{\overline {E}_{\varDelta }, \overline {\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert }_{F}, \overline {\Vert \unicode{x1D64E} \Vert }_{F}, \overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}\}$ , the estimated probabilities increase with the input variable. Statistically, this implies that events characterised by high values of $\overline {E_{\varDelta }}$ , $\overline {\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }] \Vert }_{F}$ , $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ or $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ eventually evolve into extreme production events, but with different sensibility. In all our four cases, the committor function derived from $\overline {E_{\varDelta }}$ exhibits the steepest slope, followed by $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ , $\overline {\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert }_{F}$ and finally $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ . The rise in probability as $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ increases is consistent with the qualitative analysis based on (2.4). However, the increases of probability with $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ goes counter to our theoretical analysis in (2.4), which suggests that high vorticity should reduce the probability of extreme uncertainty production. The probabilities obtained directly from the training data are also shown in figure 10 as the markers. In all four cases, the committor curves obtained from the analogue Markov chain are broadly consistent with the directly sampled empirical probabilities, showing that the reduced Markov-chain model captures the main dependence of the future extreme-event probability on the predictor. This agreement is strongest for $\overline {E}_{\varDelta }$ , for which the markers closely follow the rapid increase and the subsequent plateau of the committor curve. For $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ , $\overline {\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }] \Vert }_{F}$ and $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ , the same overall monotonic trend is observed, although the empirical probabilities exhibit relatively larger deviations.

For case F1 $_{128}$ , the variation of the committor function trained using $\overline {E}_{\varDelta }$ can be divided into two parts. We find the committor curve starts from $q=0$ at $\overline {E}_{\varDelta }=0$ , i.e. at the point of minimum uncertainty energy. The committor increases as $\overline {E_{\varDelta }}$ grows, eventually converging to the highest probability platform of $q=0.9$ when $\overline {E_{\varDelta }}=4$ . Above this value of extremely high uncertainty energy, the committor is constant. This means that most high uncertainty production event are found at large values of positive standard deviation. The committor function trained using $\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }] \Vert _{F}$ monotonically increases from 0.1 to 0.6 as $\overline {\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert }_{F}$ increases from 0 to 5, while for the strain rate $\Vert \unicode{x1D64E} \Vert _{F}$ , the committor function also increases from 0 to 0.6 over the range $\overline {\Vert \unicode{x1D64E} \Vert }_{F}\in [0,5]$ . It should be noted, however, that the effective onset for the strain-rate predictor is at $\overline {\Vert \unicode{x1D64E} \Vert }_{F}\approx 1$ , below which the committor remains nearly zero. In contrast, the committor obtained from $\Vert \boldsymbol{\varOmega } \Vert _{F}$ is the flattest among the three, increasing only from approximately 0.1 to 0.4 as $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ varies from 0 to 5. The change of slope of the committor function indicates that the probability of evolving into extreme uncertainty production is most sensitive to the uncertainty energy and subsequently to the strain rate. For case F2 $_{128}$ , the slopes of the committor function are steeper than for case F1 $_{128}$ . For case F2 $_{128}$ , the committor function trained using $E_{\varDelta }$ increases from 0 at $\overline {E}_{\varDelta }=0$ to $q=0.9$ for $\overline {E_{\varDelta }}\gt 3$ . The committor function trained using $\Vert [\unicode{x1D64E} , \boldsymbol{\varOmega }]\Vert _{F}$ monotonically increases from 0 to 0.7 as $\overline {\Vert [\unicode{x1D64E} , \boldsymbol{\varOmega }] \Vert }_{F}$ varies from 0 to 4, while for the strain rate $\Vert \unicode{x1D64E} \Vert _{F}$ , the committor function increases from 0 at $\overline {\Vert \unicode{x1D64E} \Vert }_{F}=1$ to $q=0.85$ for $\overline {\Vert \unicode{x1D64E} \Vert }_{F}\gt 6$ . All in all, the committor functions in case F2 $_{128}$ appear more sensitive than in case F1 $_{128}$ .

For the higher-Reynolds-number cases, it can be observed that the slopes of the committor functions become systematically weaker in both F1 $_{512}$ and F2 $_{512}$ . Although the four predictors retain the same qualitative ordering as in the lower-Reynolds-number cases, namely a monotonic increase of the committor and the same hierarchy of plateau levels, the plateau probabilities are all reduced compared with their lower-Reynolds-number counterparts. This reduction is particularly obvious in case F2 $_{512}$ . Part of this trend can be attributed to the lower probability of extreme events, especially those in set $\mathcal{B}$ , as reported in table 2. However, this effect alone does not explain why the magnitude of the reduction differs across predictors and across cases. One possible reason is that, at higher Reynolds numbers, the small-scale velocity-gradient structure becomes more fully developed and more strongly coupled, so that the future evolution is less controlled by any single local velocity-gradient predictor taken in isolation. More importantly, uncertainty growth is influenced not only by the magnitude of the strain rate, but also by how the strain-rate field is aligned with the local uncertainty. At higher Reynolds numbers, such alignment effects may become increasingly important as the velocity-gradient dynamics becomes more strongly coupled. This may further weaken the predictive skill of predictors based on a single scalar measure, and thus contribute to the lower plateau values of the committor functions in F1 $_{512}$ and F2 $_{512}$ .

Figure 11.

Committor function for (a,c,e) case F1 $_{128}$ and (b,d, f) case F1 $_{512}$ , estimated with different predictors. (a,b) Normalised strain rate $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ and normalised vorticity $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ , respectively. The dashed line indicates $Q=0$ and separates the strain-dominated region ( $Q\lt 0$ ) from the rotation-dominated region ( $Q\gt 0$ ). (c,d) Normalised vortex stretching $\overline {\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert }_{F}$ and normalised strain rate $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ , respectively. (e, f) Normalised vorticity $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ and normalised vortex-straining $\overline {\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }]\Vert }_{F}$ , respectively. The normalisation is defined as $\overline {F}\equiv (F)/\sigma _{F}$ . We use all the testing boxes (see § 4.3.4) at initial time as input.

Figure 12.

Committor function for (a,c,e) case F2 $_{128}$ and (b,d, f) case F2 $_{512}$ estimated with different predictors. (a,b) Normalised strain rate $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ and normalised vorticity $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ , respectively. The dashed line indicates $Q=0$ and separates the strain-dominated region ( $Q\lt 0$ ) from the rotation-dominated region ( $Q\gt 0$ ). (c,d) Normalised vortex stretching $\overline {\Vert [\unicode{x1D64E}, \boldsymbol{\varOmega }]\Vert }_{F}$ and normalised strain rate $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ , respectively. (e, f) Normalised vorticity $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ and normalised vortex-straining $\overline {\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }]\Vert }_{F}$ , respectively. The normalisation is defined as $\overline {F}\equiv (F)/\sigma _{F}$ . We use all the testing boxes (see § 4.3.4) at initial time as input.

Figure 12. Long description

Panel A: A scatter plot showing the relationship between the normalized Frobenius norm of the strain rate and the normalized Frobenius norm of the vorticity. The horizontal axis represents the normalized Frobenius norm of the strain rate, and the vertical axis represents the normalized Frobenius norm of the vorticity. The data points are color-coded based on the value of q, ranging from 0 to 0.6. A dashed line indicates Q equals 0, separating the strain-dominated region from the rotation-dominated region. Panel B: A scatter plot showing the relationship between the normalized Frobenius norm of the strain rate and the normalized Frobenius norm of the vorticity for case F2. The horizontal axis represents the normalized Frobenius norm of the strain rate, and the vertical axis represents the normalized Frobenius norm of the vorticity. The data points are color-coded based on the value of q, ranging from 0 to 0.6. A dashed line indicates Q equals 0, separating the strain-dominated region from the rotation-dominated region. Panel C: A scatter plot showing the relationship between the normalized vortex stretching and the normalized Frobenius norm of the strain rate. The horizontal axis represents the normalized vortex stretching, and the vertical axis represents the normalized Frobenius norm of the strain rate. The data points are color-coded based on the value of q, ranging from 0 to 0.6. Panel D: A scatter plot showing the relationship between the normalized vortex stretching and the normalized Frobenius norm of the strain rate for case F2. The horizontal axis represents the normalized vortex stretching, and the vertical axis represents the normalized Frobenius norm of the strain rate. The data points are color-coded based on the value of q, ranging from 0 to 0.6. Panel E: A scatter plot showing the relationship between the normalized vorticity and the normalized vortex-straining. The horizontal axis represents the normalized vorticity, and the vertical axis represents the normalized vortex-straining. The data points are color-coded based on the value of q, ranging from 0 to 0.6. Panel F: A scatter plot showing the relationship between the normalized vorticity and the normalized vortex-straining for case F2. The horizontal axis represents the normalized vorticity, and the vertical axis represents the normalized vortex-straining. The data points are color-coded based on the value of q, ranging from 0 to 0.6.

Figures 11 and 12 show two-dimensional plots of committor functions trained using two normalised scalar variables as predictors. Each point in the plot corresponds to values of the used predictors and the colour of the point indicates the value of the committor. It can be seen that the distributions of $\{\overline {\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }]\Vert }_{F}, \overline {\Vert \boldsymbol{\varOmega }\Vert }_{F}\}$ and $\{\overline {\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }] \Vert }_{F}, \overline {\Vert \unicode{x1D64E}\Vert }_{F} \}$ clearly trend along the diagonal, suggesting a correlation between the fluctuations of $[\unicode{x1D64E},\boldsymbol{\varOmega }]$ with those of $\boldsymbol{\varOmega }$ and $\unicode{x1D64E}$ (fluctuations above and below their respective averages). This agrees with our aforementioned observation that the inclusion of $[\unicode{x1D64E},\boldsymbol{\varOmega }]$ in predictor combinations does not enhance prediction performance. Notably, in all four cases, the statistical outcomes derived from the committor function for the set $\{\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}, \overline {\Vert \unicode{x1D64E} \Vert }_{F} \}$ demonstrate that the probability of extreme events occurring in the future exhibits an increasing trend from the lower left to the upper right along the diagonal, suggesting that high $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}+\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ values are more likely to trigger future extreme uncertainty production. Furthermore, it can be seen that the trained committor functions are not symmetric with respect to the dashed line corresponding to $Q=0$ . This line separates the strain-dominated region, where $Q\lt 0$ , from the rotation-dominated region, where $Q\gt 0$ . For the same value of $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F} + \overline {\Vert \unicode{x1D64E} \Vert }_{F}$ , the probability of evolving towards an extreme uncertainty-production event is significantly higher on the strain-dominated side, i.e. in the region below the dashed line, as shown in figures 11(a), 11(b), 12(a) and 12(b). This indicates that appearance of extreme uncertainty production events are more likely when the strain rate is above average than when the vorticity is above average. Similar observations can be made for the other figures, where the probability increases with increasing $\overline {\Vert \unicode{x1D64E} \Vert }_{F}$ and decreasing $\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}$ for a given value of $\overline {\Vert [\unicode{x1D64E},\boldsymbol{\varOmega }]\Vert }_{F}$ . Such results indicate that extreme high strain rate is the most important factor for extreme positive uncertainty production.

To further clarify this point, we estimate the committor function at any state of turbulence existing in the phase space. Therefore, instead of displaying the committor for predictors averaged over testing boxes (as described in § 4.3.4), we compute this probability as a function of $\{\overline {\Vert \boldsymbol{\varOmega } \Vert }_{F}, \overline {\Vert \unicode{x1D64E} \Vert }_{F} \}$ at every grid point of the DNS (i.e. $N^3$ data points) at the starting time, for all four cases. For each of these input values, we obtain the probability by using the nearest neighbour method in § 3.3 and the committor functions. Figure 13 presents colour levels of our four interpolated committors. For all four cases, the committor functions reveal a complex region where both strain rate and vorticity are slightly larger than their spatial averaged value. In this region, the estimated committor function varies rapidly without a simple trend. Not only does the Navier–Stokes system have a strong deterministic sensitivity to initial conditions because of its chaoticity, the probability of evolution towards either an extreme uncertainty production event or a much more common uncertainty reducing event may also be very sensitive to the initial condition. Situations where both strain rate and vorticity are close to their spatial averaged value may display what has been termed probabilistic unpredictability in stochastic systems (Lucente et al. Reference Lucente, Herbert and Bouchet2022a ) or it may be that more predictor information is required, for example, in terms of the pressure field and two-point correlations rather than just one-point quantities. Note that the surface region of probabilistic unpredictability is qualitatively similar in all four cases, despite the differences in forcing and Reynolds number. For the region away from the complex region, where the strain rate or vorticity is close to the two tails of their p.d.f., it is observed that the contours of probability go through the complex region. In this region, more definitive probability predictions can be made. Similar to figures 11(a), 11(b), 12(a) and 12(b), the region with high velocity gradients, particularly on the side of high strain rates, is more likely to evolve into extreme uncertainty production events, while low probabilities are mainly concentrated on the high vorticity side with smaller strain rates.

Figure 13.

Probability (in colours) based on the estimated committor function for (a) F1 $_{128}$ , (b) F1 $_{512}$ , (c) F2 $_{128}$ and (d) F2 $_{512}$ , using the input $\{\Vert \boldsymbol{\varOmega } \Vert _{F}, \Vert \unicode{x1D64E} \Vert _{F}\}$ at every grid point of the DNS at its corresponding starting time of sampling. The vertical and horizontal dash-dotted lines denote the spatial average of $\Vert \unicode{x1D64E} \Vert _{F}$ and $\Vert \boldsymbol{\varOmega } \Vert _{F}$ , respectively.

In this section, we have investigated the future evolution of uncertainty production using the statistical committor function. Among all the candidate predictors suggested by the theoretical analysis, the strain rate is found to be the quantity most strongly related, in statistical terms, to future extreme positive uncertainty production in the lower-Reynolds-number cases, as indicated by its lowest Brier score. For the higher-Reynolds-number cases, however, this statistical advantage becomes much less obvious, as all three candidate predictors give very similar Brier scores. Furthermore, the non-uniform distribution of the probability of extreme events is notable in figure 13 and corresponds to the observation that extreme uncertainty production probabilities are localised in high convergence regions (§ 5). It is worth stressing that even probabilities may not be predictable in regions of the flow field where strain rate and vorticity are near their flow field average. In such regions, the committor function is very dependent on the input predictor values, making it potentially impossible to provide stable probabilistic forecasts.