A variational approach to probing extreme events in turbulent dynamical systems

TL;DR

This paper introduces a variational framework to identify triggers of extreme events in high-dimensional nonlinear dynamical systems, demonstrated on turbulent flows, enabling early prediction of energy dissipation bursts.

Contribution

The paper presents a novel variational approach to detect and predict extreme events in turbulent systems by identifying key energy transfer mechanisms as triggers.

Findings

Extreme energy dissipation events are driven by spontaneous energy transfer from large scales.

The variational method successfully identifies the responsible triad for extreme events.

A data-driven short-term predictor for energy bursts is developed and validated.

Abstract

Extreme events are ubiquitous in a wide range of dynamical systems, including turbulent fluid flows, nonlinear waves, large scale networks and biological systems. Here, we propose a variational framework for probing conditions that trigger intermittent extreme events in high-dimensional nonlinear dynamical systems. We seek the triggers as the probabilistically feasible solutions of an appropriately constrained optimization problem, where the function to be maximized is a system observable exhibiting intermittent extreme bursts. The constraints are imposed to ensure the physical admissibility of the optimal solutions, i.e., significant probability for their occurrence under the natural flow of the dynamical system. We apply the method to a body-forced incompressible Navier--Stokes equation, known as the Kolmogorov flow. We find that the intermittent bursts of the energy dissipation are independent of the external forcing and are instead caused by the spontaneous transfer of energy from large scales to the mean flow via nonlinear triad interactions. The global maximizer of the corresponding variational problem identifies the responsible triad, hence providing a precursor for the occurrence of extreme dissipation events. Specifically, monitoring the energy transfers within this triad, allows us to develop a data-driven short-term predictor for the intermittent bursts of energy dissipation. We assess the performance of this predictor through direct numerical simulations.

Tables

Table 1. Table S1: Simulation parameters including the Reynolds number R ​ e 𝑅 𝑒 Re , the resolution N × N 𝑁 𝑁 N\times N , the mean 𝔼 ​ [ D ] 𝔼 delimited-[] 𝐷 \mathbb{E}[D] and the standard deviation 𝔼 ​ [ D 2 ] − 𝔼 ​ [ D ] 2 𝔼 delimited-[] superscript 𝐷 2 𝔼 superscript delimited-[] 𝐷 2 \sqrt{\mathbb{E}[D^{2}]-\mathbb{E}[D]^{2}} of the energy dissipation rate D 𝐷 D . The eddy turn-over time t e subscript 𝑡 𝑒 t_{e} and the prediction time t i subscript 𝑡 𝑖 t_{i} are reported in terms of non-dimensional time units. The percentage of hits, correct rejections and the false positives and negatives of the extreme event predictions are also reported. The rate of successful predictions (RSP) and the rate of successful rejections (RSR) are computed from formula ( S37 ) and ( S38 ).

$Re$	40	60	80	100
$N$	128	256	256	256
$𝔼\left[D\right]$	0.1168	0.1159	0.1010	0.0903
$\sqrt{𝔼\left[D^2\right]-𝔼{\left[D\right]}^2}$	0.0384	0.0465	0.0369	0.0295
$t_e$	0.46	0.38	0.35	0.33
$t_i$	1.0	1.0	1.0	1.0
$t_f$	2.0	2.0	2.0	1.5
Hits	$5.60\%$	$17.7\%$	$15.3\%$	$11.3\%$
Correct Rejection	$93.3\%$	$77.8\%$	$78.5\%$	$81.7\%$
False Negatives	$0.26\%$	$2.3\%$	$3.5\%$	$4.3\%$
False Positives	$0.85\%$	$2.1\%$	$2.6\%$	$2.6\%$
RSP	95.6%	88.4%	81.2%	72.3%
RSR	99.1%	97.4%	96.8%	96.9%

Equations

$${\partial }_t\mathbf{u}=\mathcal{N}(\mathbf{u}),$$

$$\mathcal{K}(\mathbf{u})=0,$$

$$\mathbf{u}(\cdot ,t_0)={\mathbf{u}}_0(\cdot ),$$

$$\sup_{\bm{u}_{0}\in X}\big{(}I(\bm{u}(t_{0}+\tau))-I(\bm{u}(t_{0}))\big{)},$$

$$\text{mbox}where\{\begin{array}{l}\mathbf{u}(t)\text{mbox}satisfies\text{eqref}eq:pde,\\ \mathbf{C}({\mathbf{u}}_0)={\mathbf{c}}_0,\end{array}$$

$$\sup_{\bm{u}_{0}\in X}\frac{\mathrm{d}}{\mathrm{d}t}\Big{|}_{t=t_{0}}I\big{(}\bm{u}(t)\big{)},$$

$$\text{mbox}where\{\begin{array}{l}\mathbf{u}(t)\text{mbox}satisfies\text{eqref}eq:pde,\\ \mathbf{C}({\mathbf{u}}_0)={\mathbf{c}}_0.\end{array}$$

$$\underset{\mathbf{u}\in X}{\sup }J(\mathbf{u}),$$

$$\text{mbox}subjectto\{\begin{array}{l}\mathcal{K}(\mathbf{u})=0,\\ \mathbf{C}(\mathbf{u})={\mathbf{c}}_0,\end{array}$$

$$J(\mathbf{u}):=\mathrm{d}I(\mathbf{u};\mathcal{N}(\mathbf{u})).$$

$$J^{\prime }(\mathbf{u})+{\mathcal{K}}^{†}(\alpha )+\sum _{i=1}^k{\beta }_i{C}_i^{\prime }(\mathbf{u})=0,$$

$$\mathcal{K}(\mathbf{u})=0,$$

$$\mathbf{C}(\mathbf{u})={\mathbf{c}}_0.$$

$${\partial }_t\mathbf{u}=-\mathbf{u}\cdot \nabla \mathbf{u}-\nabla p+\nu \Delta \mathbf{u}+\mathbf{f},\nabla \cdot \mathbf{u}=0,$$

$$E(\mathbf{u})=\frac{1}{|\Omega |}{\int }_{\Omega }\frac{|\mathbf{u}{|}^2}{2}\mathrm{d}\mathbf{x},D(\mathbf{u})=\frac{\nu }{|\Omega |}{\int }_{\Omega }|\nabla \mathbf{u}{|}^2\mathrm{d}\mathbf{x},I(\mathbf{u})=\frac{1}{|\Omega |}{\int }_{\Omega }\mathbf{u}\cdot \mathbf{f}\mathrm{d}\mathbf{x},$$

$$C(\mathbf{u}):=\frac{1}{|\Omega |}{\int }_{\Omega }\frac{|A(\mathbf{u}){|}^2}{2}\mathrm{d}\mathbf{x},$$

$$J(\bm{u})=\frac{1}{|\Omega|}\int_{\Omega}\big{[}\bm{u}\cdot(\bm{u}\cdot\boldsymbol{\nabla}\bm{f})+\nu\,\bm{u}\cdot\left(\Delta\bm{f}\right)\big{]}\mathrm{d}\bm{x}.$$

$$\left(\nabla \mathbf{f}+\nabla {\mathbf{f}}^{\top }\right)\mathbf{u}+\nu \Delta \mathbf{f}-\nabla \alpha +\beta A^{†}A\mathbf{u}=0,$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{1}{|\Omega |}{\int }_{\Omega }\frac{|A(\mathbf{u}){|}^2}{2}\mathrm{d}\mathbf{x}=c_0,$$

$$\displaystyle\mathcal{P}\left(D_{1}\leq D_{m}(t)\leq D_{2}\big{|}\,\lambda_{1}\leq\lambda(t)\leq\lambda_{2}\right),$$

$$P_{ee}\left(\lambda_{0}\right)=\mathcal{P}\left(D_{m}(t)>D_{e}\big{|}\,\lambda(t)=\lambda_{0}\right),$$

$${\mathcal{L}}_c(\mathbf{u},\alpha ,\mathbf{\beta }):=J(\mathbf{u})+⟨\mathcal{K}(\mathbf{u}),\alpha {⟩}_X+\mathbf{\beta }\cdot (\mathbf{C}(\mathbf{u})-{\mathbf{c}}_0),$$

$$\frac{\delta {\mathcal{L}}_c}{\delta \mathbf{u}}(\mathbf{v})$$

$$=\mathrm{d}J(\mathbf{u};\mathbf{v})+⟨\mathcal{K}(\mathbf{v}),\alpha ⟩+\sum _{i=1}^k{\beta }_i\mathrm{d}{C}_i(\mathbf{u};\mathbf{v})$$

$$=⟨J^{\prime }(\mathbf{u}),\mathbf{v}⟩+⟨{\mathcal{K}}^{†}(\alpha ),\mathbf{v}⟩+\sum _{i=1}^k{\beta }_i⟨C_i^{\prime }(\mathbf{u}),\mathbf{v}⟩$$

$$\displaystyle=\langle J^{\prime}(\bm{u})+\mathcal{K}^{\dagger}(\alpha)+{\color[rgb]{0,0,0}\sum_{i=1}^{k}\beta_{i}C^{\prime}_{i}(\bm{u})},\bm{v}\rangle.$$

$$J^{\prime}(\bm{u})+\mathcal{K}^{\dagger}(\alpha)+{\color[rgb]{0,0,0}\sum_{i=1}^{k}\beta_{i}C^{\prime}_{i}(\bm{u})}=\bm{0}.$$

$$\frac{\delta {\mathcal{L}}_c}{\delta \alpha }(\tilde{\alpha })=⟨\mathcal{K}(\mathbf{u}),\tilde{\alpha }{⟩}_X,\frac{\delta {\mathcal{L}}_c}{\delta \mathbf{\beta }}=\mathbf{C}(\mathbf{u})-{\mathbf{c}}_0.$$

$${\partial }_t\mathbf{u}=-\mathbf{u}\cdot \nabla \mathbf{u}-\nabla p+\nu \Delta \mathbf{u}+\mathbf{f},$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{u}(\mathbf{x},t)=\sum _{\mathbf{k}\in {\mathbb{Z}}^2}\frac{a(\mathbf{k},t)}{k}\left(\begin{array}{c}{k}_2\\ -k_1\end{array}\right)e^{\hat{i}\mathbf{k}\cdot \mathbf{x}},$$

$$I(\mathbf{u}(t))=-\text{mbox}Im[a({\mathbf{k}}_f,t)]=-r({\mathbf{k}}_f,t)\sin \left(\varphi ({\mathbf{k}}_f,t)\right),$$

$$\mathbf{u}(\mathbf{x},t)=\overline{\mathbf{u}}(\mathbf{x},t)+\sum _{j=1}^{\infty }{\alpha }_j(t){\mathbf{v}}_j(\mathbf{x}),$$

$${\partial }_t{\hat{u}}_i(\mathbf{k})=-\hat{i}{P}_{ij}(\mathbf{k})\sum _{\mathbf{p}+\mathbf{q}=\mathbf{k}}{q}_m{\hat{u}}_m(\mathbf{p}){\hat{u}}_j(\mathbf{q})-\nu k^2{\hat{u}}_i(\mathbf{k})+{\hat{f}}_i(\mathbf{k}),$$

$${\partial }_t{\hat{u}}_1(\mathbf{k})=$$

$${\partial }_t{\hat{u}}_2(\mathbf{k})=$$

Full text

A variational approach to probing extreme events in turbulent dynamical systems

Mohammad Farazmand,∗ Themistoklis P. Sapsis*∗*

Department of Mechanical Engineering, Massachusetts Institute of Technology

77 Massachusetts Avenue, Cambridge MA 02139-4307

∗To whom correspondence should be addressed; E-mail: mfaraz@mit.edu, sapsis@mit.edu

Extreme events are ubiquitous in a wide range of dynamical systems, including turbulent fluid flows, nonlinear waves, large scale networks and biological systems. Here, we propose a variational framework for probing conditions that trigger intermittent extreme events in high-dimensional nonlinear dynamical systems. We seek the triggers as the probabilistically feasible solutions of an appropriately constrained optimization problem, where the function to be maximized is a system observable exhibiting intermittent extreme bursts. The constraints are imposed to ensure the physical admissibility of the optimal solutions, i.e., significant probability for their occurrence under the natural flow of the dynamical system. We apply the method to a body-forced incompressible Navier–Stokes equation, known as the Kolmogorov flow. We find that the intermittent bursts of the energy dissipation are independent of the external forcing and are instead caused by the spontaneous transfer of energy from large scales to the mean flow via nonlinear triad interactions. The global maximizer of the corresponding variational problem identifies the responsible triad, hence providing a precursor for the occurrence of extreme dissipation events. Specifically, monitoring the energy transfers within this triad, allows us to develop a data-driven short-term predictor for the intermittent bursts of energy dissipation. We assess the performance of this predictor through direct numerical simulations.

Introduction

A plethora of dynamical systems exhibit intermittent behavior manifested through sporadic bursts in the time series of their observables. These extreme events produce values of the observable that are several standard deviations away from its mean, resulting in heavy tails of the corresponding probability distribution. Important examples include climate phenomena  (?, ?), rogue waves in oceanic and optical systems  (?, ?, ?) and large deviations in turbulent flows  (?, ?, ?). Since such extreme phenomena typically have dramatic consequences, their quantification and prediction is of great interest.

Significant progress has been made in the computation of extreme statistics, both through direct numerical simulations and through indirect methods such as the theory of large deviations (see  (?) for a review). While these methods estimate the probability distribution of the extreme events, they do not inform us about the underlying mechanisms that lead to the extremes nor are they capable of predicting individual extreme events.

For systems operating near an equilibrium or systems that are nearly integrable, perturbative methods have been successful in identifying the resonant interactions that cause the extreme events  (?, ?). For systems which are not perturbations from such trivial limits, a general framework for probing the transition mechanism to extreme states is missing. These systems, such as turbulent fluid flows and water waves, are also typically high-dimensional and nonlinear where the nonlinearities create a complex network of interdependent interactions among many degrees of freedom  (?, ?, ?, ?, ?, ?).

Here, we propose a variational framework to probe the underlying conditions that lead to extreme events in such high-dimensional complex systems. More specifically, we derive precursors of extreme events as the solutions of a finite-time constrained optimization problem. The functional to be maximized is the observable whose time series exhibit the extreme events. The constraints are designed to ensure that the triggers belong to the system attractor and therefore reflect physically relevant phenomena. If the life-time of the extreme events are short compared to the typical dynamical time scales of the system, the finite-time optimization problem can be replaced with its instantaneous counterpart. For the instantaneous problem, we derive the Euler-Lagrange equations that can be solved numerically using Newton-type iterations.

We apply the variational framework to the Kolmogorov flow, a two-dimensional Navier–Stokes equation driven by a monochromatic body forcing. At sufficiently high Reynolds numbers, this flow is known to exhibit intermittent bursts of energy dissipation rate  (?, ?). We first show that these extreme bursts are due to the internal transfers of energy through nonlinearities, as opposed to phase locking with the external forcing. Because of the high number of involved degrees of freedom and their complex interactions, deciphering the modes responsible for the extreme energy dissipation rate is not straightforward. The optimal solution to our variational method, however, isolates the triad interaction responsible for this transfer of energy. Monitoring this triad along trajectories of the Kolmogorov flow, we find that, on the onset of the extreme bursts, the energy is transferred spontaneously from a large-scale Fourier mode to the mean flow, leading to growth of the energy input rate and consequently the energy dissipation rate.

We then utilize the derived large-scale mode as a predictor for the extreme events. Specifically, by tracking the energy of this mode we develop a data-driven short-term predictor of intermittency in the Kolmogorov flow. We assess the effectiveness of the prediction scheme on extensive direct numerical simulations by explicitly quantifying its success rate, as well as the false positive and false negative rates.

Results

Variational formulation of extreme events

Consider the general evolution equations,

[TABLE]

where $\bm{u}:\Omega\times\mathbb{R}^{+}\to\mathbb{R}^{d}$ belongs to an appropriate function space $X$ and completely determines the state of the system. The initial condition $X\ni\bm{u}_{0}:\Omega\to\mathbb{R}^{d}$ is specified at the time $t_{0}$ and $\Omega\subset\mathbb{R}^{d}$ is an open bounded domain. The differential operators $\mathcal{N}$ and $\mathcal{K}$ can be potentially nonlinear. A wide range of physical models can be written as a set of partial differential equations (PDEs) as in (1c). For instance, for incompressible fluid flows,  (1a) is the momentum equation and (1b) is the incompressibility condition where $\mathcal{K}(\bm{u})=\boldsymbol{\nabla}\cdot\bm{u}$. For simplicity, we will denote a trajectory of (1c) by $\bm{u}(t)$.

Let $I:X\rightarrow\mathbb{R}$ denote an observable whose time series $I(\bm{u}(t))$ along a typical trajectory $\bm{u}(t)$ exhibits intermittent bursts (see figure 1a). Drawing upon the near-integrable case, we view the system as consisting of a background chaotic attractor which has small regions of instability (see figure 1b). Once a trajectory reaches an instability region, it is momentarily repelled away from the background attractor, resulting in a burst in the time series of the observable. Our goal here is to probe the instability region(s) by utilizing a combination of observed data from the system and the governing equations of the system. We also require the instability regions to have non-zero probability of occurrence under the natural flow of the dynamical system. This constraint is of particular importance since it excludes ‘exotic’ states with extreme growth of $I$ but with negligible probability of being observed in practice (see the constraint $\bm{C}(\bm{u}_{0})=\bm{c}_{0}$ in equation (2b) below).

We formulate this task as a constrained optimization problem. Assume that there is a typical timescale $\tau\in\mathbb{R}^{+}$ over which the bursts in the observable $I$ develop (see figure 1a). We seek initial conditions $\bm{u}_{0}$ whose associated observable $I(\bm{u}(t))$ attains a maximal growth within time $\tau$. More precisely, we seek the solutions to the constrained optimization problem,

[TABLE]

where the optimization variable is the initial condition $\bm{u}(t_{0})=\bm{u}_{0}$ of system (1c). The set of critical states are required to satisfy the constraints in (2b) in order to enforce two important properties. The first property ensures that $\bm{u}(t)$ obeys the governing equation (1c) as opposed to being an arbitrary one-parameter family of functions. The second property $\bm{C}(\bm{u}_{0})=\mathbf{c}_{0}$, where $\bm{C}:X\to\mathbb{R}^{k}$, is a codimension-$k$ constraint. This constraint is enforced to ensure the non-zero probability of occurrence, i.e. states that are sufficiently close to the chaotic background attractor. The set of probabilistically feasible states can be generally described by exploiting basic physical properties of the chaotic attractor such as average energy along different components of the state space or the second-order statistics. The precise form of the constraint $\bm{C}(\bm{u}_{0})=\bm{c}_{0}$ is problem dependent and will shortly be discussed in more detail. We point out that more general inequality constraints of the form $\bm{c}_{\mathrm{min}}\leq\bm{C}(\bm{u}_{0})\leq\bm{c}_{\mathrm{max}}$ may also be employed. The treatment of such inequality constraints, however, is not discussed in this paper.

We expect the set of solutions to problem (2b) to unravel the mechanisms underpinning the intermittent bursts of the observable. Although it is unlikely that a generic trajectory of the system passes exactly through one of the maximizers, by continuity, any trajectory passing through a sufficiently small open neighborhood of the maximizer (i.e., the instability regions of figure 1b) will result in a similar observable burst.

We emphasize that an optimization problem similar to (2b) has been pursued before in special contexts. The largest finite-time Lyapunov exponent can be formulated as (2b) where the observable is the amplitude of infinitesimal perturbations after finite time. The maximizer is the corresponding finite-time Lyapunov vector  (?, ?). In a similar context, Pringle and Kerswell  (?) seek optimal finite-amplitude perturbations that trigger transition to turbulence in the pipe flow. They formulate the unknown optimal perturbation as the solution of a constrained optimization problem similar to (2b), with the observable being the $L^{2}$ norm of the fluid velocity field. Ayala and Protas  (?, ?, ?) consider the finite-time singularity formation for Navier–Stokes equations. They also use a variational method to seek the initial conditions that could lead to finite-time singularities. In these studies, the emphasis is given to the analysis of the most ‘unstable’ states but the physical properties of the attractor are not taken into account.

The standard approach for solving the PDE-constrained optimization problem (2b) is an adjoint-based gradient iterative method  (?, ?, ?). This method is computationally very expensive since, at each iteration, the gradient direction needs to be evaluated as the solution of an adjoint PDE. If the growth timescale $\tau$ is small compared to the typical time scales of the observable, it is reasonable to replace the finite-time growth problem (2b) with its instantaneous counterpart,

[TABLE]

Problem (3b) seeks initial states $\bm{u}_{0}$ for which the instantaneous growth of the observable $I$ along the corresponding solution $\bm{u}(t)$ is maximal.

We point out that the large instantaneous derivatives of $I$ do not necessarily imply a subsequent burst in the observable as the growth may not always be sustained at later times along the trajectory $\bm{u}(t)$. As a result, the set of solutions to this instantaneous problem may differ significantly from the finite-time problem (2b). Nonetheless, the solutions to the instantaneous problem can still be insightful. In addition, as we show below, these solutions can be obtained at a relatively low computational cost.

Optimal solutions

First, we derive an equivalent form of problem (3b). Taking the time derivative of the time series $I(\bm{u}(t))$ yields $(\mathrm{d}/\mathrm{d}t)I(\bm{u}(t))=\mathrm{d}I(\bm{u};\partial_{t}\bm{u})$, where $\mathrm{d}I(\bm{u};\bm{v}):=\lim_{\varepsilon\to 0}\left[I(\bm{u}+\varepsilon\bm{v})-I(\bm{u})\right]/\varepsilon$ denotes the Gâteaux differential of $I$ at $\bm{u}$ evaluated along $\bm{v}$. Using (1a) to substitute for $\partial_{t}\bm{u}$, we obtain the following optimization problem which is equivalent to problem (3b):

[TABLE]

where

[TABLE]

Note that the first constraint in (3b) simplified since we have already used (1a) and it only remains to enforce (1b). For notational simplicity, we omit the subscript from $\bm{u}_{0}$.

If $J:X\to\mathbb{R}$ is a continuous map and the subset $S=\{\bm{u}\in X:\mathcal{K}(\bm{u})=0,\;\bm{C}(\bm{u})=\bm{c}_{0}\}$ is compact in $X$, problem (4b) has at least one solution. This follows from the fact that the image of a compact set under a continuous transformation is compact. Therefore $J(S)\subset\mathbb{R}$ is compact which implies $J(S)$ is bounded and closed  (?). Therefore $J$ is bounded and attains its maximum (and minimum) on $S$. The uniqueness of the maximizer is not generally guaranteed. However, the set of maximizers (and minimizers) of $J$ are compact subsets of $S$  (?).

As we show in the Supplementary Material (section S1), if $X$ is a Hilbert space with the inner product $\langle\cdot,\cdot\rangle_{X}$ and the operator $\mathcal{K}$ is linear, every solution of the optimization problem (4b) satisfies the set of Euler–Lagrange equations,

[TABLE]

Here $\mathcal{K}^{\dagger}$ is the adjoint of $\mathcal{K}$ and $J^{\prime}(\bm{u})$ and $C_{i}^{\prime}(\bm{u})$ are the unique identifier of the Gâteaux differentials $\mathrm{d}J(\bm{u};\cdot)$ and $\mathrm{d}C_{i}(\bm{u};\cdot)$, such that $\mathrm{d}J(\bm{u};\bm{v})=\langle J^{\prime}(\bm{u}),\bm{v}\rangle_{X}$ and $\mathrm{d}C_{i}(\bm{u};\bm{v})=\langle C_{i}^{\prime}(\bm{u}),\bm{v}\rangle_{X}$ for all $\bm{v}$. The existence and uniqueness of $J^{\prime}(\bm{u})$ and $C_{i}^{\prime}(\bm{u})$ are guaranteed by the Riesz representation theorem  (?). Here, $C_{i}$ are the components of the map $\bm{C}=(C_{1},C_{2},\cdots,C_{k})$. The function $\alpha:\Omega\to\mathbb{R}$ and the vector $\bm{\beta}=(\beta_{1},\cdots,\beta_{k})\in\mathbb{R}^{k}$ are unknown Lagrange multipliers to be determined together with the optimal state $\bm{u}:\Omega\to\mathbb{R}^{d}$.

Application to Navier–Stokes equations

We consider the Navier–Stokes equations

[TABLE]

where $\bm{u}:\Omega\times\mathbb{R}^{+}\to\mathbb{R}^{d}$ is the fluid velocity field, $p:\Omega\times\mathbb{R}^{+}\to\mathbb{R}$ is the pressure field and $\nu=Re^{-1}$ is the non-dimensional viscosity which coincides with the reciprocal of the Reynolds number $Re$. Here, we consider two-dimensional flows ($d=2$) over the domain $\Omega=[0,2\pi]\times[0,2\pi]$ with periodic boundary conditions. The flow is driven by the monochromatic Kolmogorov forcing $\bm{f}(\bm{x})=\sin(k_{f}y)\bm{e}_{1}$ where $\bm{k}_{f}=(0,k_{f})$ is the forcing wave number and the vectors $\bm{e}_{i}$ denote the standard basis in $\mathbb{R}^{d}$. In the following, we assume that the velocity fields are square integrable for all times, i.e., $X=L^{2}(\Omega)$.

The kinetic energy $E$ (per unit volume), the energy dissipation rate $D$ and the energy input rate $I$ are defined, respectively, by

[TABLE]

where $|\Omega|$ denotes the area of the domain, i.e., $|\Omega|=(2\pi)^{2}$. Along any trajectory $\bm{u}(t)$ these three quantities satisfy $\dot{E}=I-D$. We use the energy dissipation rate $D$ to define the eddy turn-over time, $t_{e}=\sqrt{\nu/\mathbb{E}[D]}$, where $\mathbb{E}$ denotes the expected value.

The Kolmogorov flow admits the laminar solution $\bm{u}=(Re/k_{f}^{2})\sin(k_{f}y)\bm{e}_{1}$. For the forcing wave number $k_{f}=1$, the laminar solution is the global attractor of the system at any Reynolds number  (?). If the forcing is applied at a higher wave number and the Reynolds number is sufficiently large, the laminar solution becomes unstable. In particular, numerical evidence suggests that for $k_{f}=4$ and sufficiently large Reynolds numbers, the Kolmogorov flow is chaotic and exhibits intermittent bursts of energy dissipation  (?, ?, ?). This is manifested in figure S10(a), showing the time series of the energy dissipation $D$ at Reynolds number $Re=40$ with $k_{f}=4$.

A closer inspection reveals that each burst of the energy dissipation $D$ is shortly preceded by a burst in the energy input $I$ (see figure S10(b)). Therefore, we expect the mechanism behind the bursts in the energy input to be also responsible for the bursts in the energy dissipation. As we show in Supplementary Material (section S2.1), the energy input is given by $I(\bm{u})=-|a(\bm{k}_{f})|\sin(\phi(\bm{k}_{f}))$ where $a(\bm{k})$ are the Fourier modes of the velocity field $\bm{u}$ and $\phi(\bm{k})$ are their corresponding phases such that $a(\bm{k})=|a(\bm{k})|\exp(i\phi(\bm{k}))$. We refer to the Fourier mode $a(\bm{k}_{f})$ as the mean flow. The energy input $I$ can grow through two mechanisms: (i) Alignment between the phase of the mean flow and the external forcing, i.e., $\phi(\bm{k}_{f})\to-\pi/2$ and (ii) Growth of the mean flow energy $|a(\bm{k}_{f})|$.

Examining the alignment between the forcing and the velocity field rules out mechanism (i) (cf. figure S6 of the Supplementary Material). The remaining mechanism (ii) is possible through the nonlinear term in the Navier–Stokes equation. This nonlinearity redistributes the system energy among various Fourier modes $a(\bm{k})$ through triad interactions of the modes whose wave numbers $(\bm{k},\bm{p},\bm{q})$ satisfy $\bm{k}=\bm{p}+\bm{q}$ (see Supplementary Material, section S2.2, for further details). Due to the high number of active modes involved in the intricate network of triad interactions, it is unclear which triad (or triads) are responsible for the nonlinear transfer of energy to the mean flow during the extreme events. As we show below, our variational approach identifies the modes involved in this transfer. Before obtaining the optimal solution, however, we need to specify the explicit form of the constraint $C(\bm{u})=c_{0}$.

Constraints

The constraint $\bm{C}(\bm{u})=\bm{c}_{0}$ is imposed in order to ensure that the optimal solutions are physically admissible, i.e. that they are sufficiently close to the attractor and thus have non-zero probability of occurrence. For instance, the solutions to a wide range of dissipative PDEs are known to converge asymptotically to a finite-dimensional subset of the state space  (?). The maximizers of the functionals (2a) or (3a) that are far from this asymptotic attractor are physically irrelevant as they correspond to a transient phase that cannot be sustained.

In most applications, including the Kolmogorov flow, the system attractor is not explicitly known. Therefore, the physical relevance of the optimal solutions need to be ensured otherwise. Here, we consider constraints of the form

[TABLE]

where $A$ is a linear operator. Several physically important quantities can be written as the function (9). For instance, if $A$ is the identity operator, $C$ coincides with the kinetic energy, $C(\bm{u})=E(\bm{u})$. If $A$ is the gradient operator, $A=\boldsymbol{\nabla}$, we have $C(\bm{u})=D(\bm{u})\times(Re/2)$ which can be used to constrain the energy dissipation rate. A more general class of such operators can be constructed as follows. Let $\bm{u}=\sum_{i}\alpha_{i}\bm{v}_{i}$ where $\{\bm{v}_{i}\}$ is a principal component basis, i.e. it diagonalizes the covariance operator of $\bm{u}$. Define $A$ as the diagonal linear operator such that $A(\bm{v}_{i})=\bm{v}_{i}/\sigma_{i}$, where $\sigma_{i}^{2}$ is the standard deviation of $\alpha_{i}$. Then the constraint takes the form $C(\bm{u})=\left(\sum_{i}\alpha_{i}^{2}/\sigma_{i}^{2}\right)/(2|\Omega|)$. This corresponds to an ellipsoid, describing points that have equal probability of occurrence when we approximate the statistics of the background attractor by a Gaussian measure  (?). Note that the constraints of the form (9) are codimension-one, $C:X\to\mathbb{R}$, and hence a special case of the codimension-$k$ constraint in equation (2b).

Excluding the intermittent bursts, the energy dissipation of the Kolmogorov flow exhibits small oscillations around its mean value $\mathbb{E}[D]$ (see figure S10a). Based on this observation, we seek optimal solutions of (4b) which are constrained to have the energy dissipation $D=\mathbb{E}[D]$. This result into the constraint (9) with $A=\boldsymbol{\nabla}$ and $C(\bm{u})=c_{0}=\mathbb{E}[D]\times(Re/2)$. We approximate the mean value $\mathbb{E}[D]$ from direct numerical simulations. At $Re=40$, for instance, we have $\mathbb{E}[D]\simeq 0.117$.

Probing the extreme energy transfers

The functional $J$ (see Eq. (5)) associated with the energy input $I$ reads

[TABLE]

The associated Euler–Lagrange equations (6c) read

[TABLE]

where $J^{\prime}(\bm{u})=\left(\boldsymbol{\nabla}\bm{f}+\boldsymbol{\nabla}\bm{f}^{\top}\right)\bm{u}+\nu\Delta\bm{f}$, $\mathcal{K}^{\dagger}(\alpha)=-\boldsymbol{\nabla}\alpha$ and $C^{\prime}(\bm{u})=A^{\dagger}A\bm{u}$ (see Supplementary Material, section S2.3, for the derivations). We set $A=\boldsymbol{\nabla}$ in order to enforce a constant energy dissipation constraint. This implies $A^{\dagger}A\bm{u}=-\Delta\bm{u}$.

Using the symmetries of Eq. (S30c), we find that it admits the pair of exact solutions $\bm{u}_{\pm}=\pm(2\sqrt{c_{0}}/k_{f})\sin(k_{f}y)\bm{e}_{1}$, $\alpha_{\pm}=\pm\sqrt{c_{0}}\int\sin(2k_{f}y)\mathrm{d}y$, $\beta_{\pm}=\pm(\nu k_{f}/2\sqrt{c_{0}})$. More complex solutions, with unknown closed forms, may exist. We approximate these solutions using the Newton iterations described in Supplementary Material, section S3.

At each $Re$, we initiated several Newton iterations from random initial conditions. In addition to the pair of exact solutions $(\bm{u}_{\pm},\alpha_{\pm},\beta_{\pm})$, the iterations yielded one non-trivial solution. Figure 3 shows the resulting three branches of solutions including the exact solution $\bm{u}_{+}$ (solid black), the exact solution $\bm{u}_{-}$ (dashed black) and the non-trivial solution (red circles). For small Reynolds numbers, our Newton searches only returned the exact solutions. At $Re\simeq 3.1$, a bifurcation takes place where a new non-trivial solution is born. This solution appears to be a global maximizer as no other solutions were found. Since the intermittent bursts are only observed for $Re>35$, we focus the following analysis on this range of Reynolds numbers.

The non-trivial optimal solution converges to an asymptotic limit as the $Re$ increases. This is discernible from the plateau of the red curve in figure 3 and the select solutions shown in its outset. The three most dominant Fourier modes present in this asymptotic solution are the forcing wave number $(0,k_{f})$ and the wave numbers $(1,0)$ and $(1,k_{f})$ together with their complex conjugate pairs. Incidentally, these wave numbers form a triad, $(0,k_{f})+(1,0)=(1,k_{f})$. The dominant mode of the optimal solution corresponds to the wave number $(1,0)$ whose modulus $|a(1,0)|$ is one order of magnitude larger than the other non-zero modes.

Next, we turn to the direct numerical simulations of the Kolmogorov flow and monitor the three Fourier modes $a(0,k_{f})$, $a(1,0)$ and $a(1,k_{f})$. We find that the energy transfers within this triad underpin the intermittent bursts of the mean flow energy $|a(0,k_{f})|$, and hence the energy input rate $I$. Figure 4(a) shows time series of $I=-\mbox{Im}[a(0,k_{f})]$ and the Fourier mode $|a(1,0)|$ along a typical trajectory of the Kolmogorov flow at $Re=40$. The bursts of the energy input rate $I$ are nearly concurrent with extreme dips in the modulus $|a(1,0)|$. A similar concurrent behavior was observed for other trajectories and at higher Reynolds numbers (see Supplementary Material, section S6).

This observation reveals that, prior to a burst, the mode $a(1,0)$ transfers a significant portion of its energy budget to the mean flow $a(0,k_{f})$ through the triad interaction of the modes $a(0,k_{f})$, $a(1,0)$ and $a(1,k_{f})$. This leads to the increase in the energy of the mean flow, and therefore the energy input rate $I$, which in turn leads to the growth of the energy dissipation rate $D$.

One can go one step further and inquire about the reason for the release of energy from mode $a(1,0)$ to the mean flow $a(0,k_{f})$. The answer to this question, involving the relative phases of the the modes and their interactions with other triads, is beyond the scope of the present work and will be addressed elsewhere. It is tempting to study these interactions by truncating the Kolmogorov flow to the modes $a(0,k_{f})$, $a(1,0)$, $a(1,k_{f})$ and their complex conjugates. Unfortunately, such severe truncations fail to be illuminating since the dynamics of the truncated system severely departs from the original Navier–Stokes equations  (?, ?).

Prediction of extreme events

Given the above observation that the optimal solution primarily consists of mode $(1,0)$, we choose this mode to formulate our data-driven prediction scheme. The decrease in the energy of the mode $(1,0)$ precedes the increase in the energy of the mean flow. This enables the data-driven short-term prediction of extreme bursts of the energy dissipation by observing the modulus $|a(1,0,t)|$. More specifically, relatively small values of $\lambda(t):=|a(1,0,t)|$, along a solution $\bm{u}(t)$, signal the high probability of an upcoming burst in the energy dissipation.

To quantify this, we consider the conditional probability

[TABLE]

where $D_{m}(t)=\max_{\tau\in[t+t_{i},t+t_{f}]}D(\bm{u}(\tau))$ is the maximum of the energy dissipation over the future time interval $[t+t_{i},t+t_{f}]$. This conditional probability measures the likelihood of the future maximum value of the energy dissipation belonging to the interval $[D_{1},D_{2}]$, given that the present value of the indicator $\lambda(t)$ belongs to $[\lambda_{1},\lambda_{2}]$. The constant parameters $t_{f}>t_{i}>0$ determine the future time interval $[t+t_{i},t+t_{f}]$. In the following, we set $t_{i}=1\simeq 2.2t_{e}$ and $t_{f}=2\simeq t_{i}+2.2t_{e}$. The length of the time window $t_{f}-t_{i}$ is long enough to ensure that the extreme event (if it exists) is contained in the time interval $[t+t_{i},t+t_{f}]$. The choice of the prediction time $t_{i}$ will be discussed shortly. The reported results are robust to small perturbations to all parameters.

Figure 5(a) shows the conditional probability density corresponding to (12). We observe that relatively small values of $\lambda$ correlate strongly with the high future values of the energy dissipation $D$. For instance, when $\lambda<0.4$, the value of $D_{m}$ is most likely larger than $0.2$. Conversely, when $\lambda$ is larger than $0.4$, the future values of the future energy dissipation $D_{m}$ are smaller than $0.2$.

We seek an appropriate value $\lambda_{0}$ such that $\lambda(t)<\lambda_{0}$ predicts an extreme burst of energy dissipation over the future time interval $[t+t_{i},t+t_{f}]$. Denote the extreme event threshold by $D_{e}$, such that $D>D_{e}$ constitutes an extreme burst of energy dissipation. We define the probability of an upcoming extreme event $P_{ee}$ as

[TABLE]

which measures the likelihood that $D_{m}(t)>D_{e}$ assuming that $\lambda(t)=\lambda_{0}$. Here, we set the threshold of the extreme event $D_{e}$ as the mean value of the energy dissipation plus two standard deviations, $D_{e}=\mathbb{E}[D]+2\sqrt{\mathbb{E}[D^{2}]-\mathbb{E}[D]^{2}}\simeq 0.194$. The extreme event probability $P_{ee}$ can be computed from the probability density shown in figure 5(a) (see Supplementary Material, section S5, for the details).

Figure 5(b) shows the probability of extreme events $P_{ee}$ as a function of the parameter $\lambda_{0}$. If at time $t$, the values of $\lambda(t)$ is larger than $0.5$ the probability of a future extreme event, i.e. $D_{m}(t)>D_{e}$, is nearly zero. The probability of a future extreme event increases as $\lambda(t)$ decreases. At $\lambda(t)\simeq 0.4$, the probability is $50\%$. If $\lambda(t)<0.3$, the likelihood of an upcoming extreme event is nearly $100\%$. The horizontal dashed line in figure 5(a) marks the transition line from the low likelihood of an upcoming extreme event $P_{ee}<0.5$ to the higher likelihood $P_{ee}>0.5$. This line together with the vertical line $D_{m}=D_{e}$ divide the conditional probability density into four regions: (I) Correct rejections ($P_{ee}<0.5$ and $D_{m}(t)<D_{e}$): Correct prediction of no upcoming extremes. (II) False positives ($P_{ee}>0.5$ but $D_{e}(t)<D_{e}$): The indicator predicts an upcoming extreme event but no extreme event actually takes place. (III) Hits ($P_{ee}>0.5$ and $D_{m}(t)>D_{e}$): Correct prediction of an upcoming extreme event. (IV) False negatives ($P_{ee}<0.5$ but $D_{m}(t)>D_{e}$): An extreme event takes place but the indicator fails to predict it.

A reliable indicator of upcoming extreme events must maximize the number of correct rejections (quadrant I) and hits (quadrant III), while at the same time having minimal false positives (quadrant II) and false negatives (quadrant IV). From nearly $100,000$ predictions made only $0.26\%$ false negatives and $0.85\%$ false positives were recorded. The number of hits were $5.6\%$ and the number of correct rejections amount to $93.3\%$ of all predictions made. As we show in the Supplementary Material, this amounts to $95.6\%$ success rate for the prediction of the extreme events (see Equation (S37) and Table S1). Note that the high percentage of correct rejections compared to the hits is a mere consequence of the fact that the extreme events are rare.

An additional desirable property of an indicator is its ability to predict the upcoming extremes well in advance of the events taking place. The chosen prediction time $t_{i}=1$ is approximately twice the eddy turn-over time $t_{e}$. In comparison, it takes approximately one eddy turn-over time (on average) for the energy dissipation rate to grow from $\mathbb{E}[D]+\sqrt{\mathbb{E}[D^{2}]-\mathbb{E}[D]^{2}}$ to its extreme value $D_{e}=\mathbb{E}[D]+2\sqrt{\mathbb{E}[D^{2}]-\mathbb{E}[D]^{2}}$.

The prediction time $t_{i}$ can always be increased at the cost of increasing false positives and/or false negatives. For instance, with the choice $t_{i}=2\simeq 4.3t_{e}$ and $t_{f}=3\simeq t_{i}+2.2t_{e}$, prediction of the extreme events $D_{m}>D_{e}$ returns $1.2\%$ false negatives and $0.6\%$ false positives. The number of hits decreases slightly to $5.3\%$, as does the number of correct rejections $92.9\%$, which amounts to a success rate of $82\%$ in the extreme event prediction (see Equation (S37)). Therefore, the prediction time $t_{i}=2$ still yields satisfactory predictions. Upon increasing $t_{i}$ further, eventually the number of hits becomes comparable to the number of false negatives at which point the predictions are unreliable.

Discussion

A method for the computation of precursors of extreme events in complex turbulent systems is introduced here. The new approach combines basic physical properties of the chaotic attractor (such as energy distribution along different directions of phase space) obtained from data, with stability properties induced by the governing equations. The method is formulated as a constrained optimization problem which can be solved explicitly if the timescale of the extreme events is short compared to the typical timescales of the system. To demonstrate the approach we consider a stringent test case, the Kolmogorov flow, which has a turbulent attractor with positive Lyapunov exponents and intermittent extreme bursts of energy dissipation. We are able to correctly identify the triad of modes associated with the extreme events. Moreover, the derived precursors allow for the formulation of an accurate short-term prediction scheme for the intermittent bursts. The results demonstrate the robustness and applicability of the approach on systems with high-dimensional chaotic attractors.

Materials and Methods

The Navier–Stokes equations and the corresponding Euler–Lagrange equations are solved numerically with a standard pseudo-spectral code with $N\times N$ Fourier modes and $2/3$ dealiasing. For $Re=60,80$ and $100$, we use $N=256$ to fully resolve the velocity fields. At $Re=40$, however, this resolution is unnecessarily high and hence we use $N=128$. The temporal integration of the Navier–Stokes equations are carried out with a forth-order Runge–Kutta scheme.

Supplementary Material

section S1. Derivation of the Euler–Lagrange equation

section S2. The Navier–Stokes equation

section S3. Newton iterations

section S4. Sensitivity to parameters

section S5. Computing the probability of extreme events

section S6. Supporting computational results

fig. S1. Evolution of the energy input vs. mean flow

fig. S2. Triad interactions

fig. S3. Sensitivity of the optimal solutions

fig. S4. Joint PDFs for higher Reynolds numbers

fig. S5. Prediction of intermittent bursts at higher Reynolds numbers

table S1: Simulation parameters

Acknowledgments: This work has been supported through the ONR grant N00014-15-1-2381, the AFOSR grant FA9550-16- 1-0231, and the ARO grant 66710-EG-YIP.

Data and materials availability: Certain codes and data used in this manuscript are publicly available on GitHub at https://github.com/mfarazmand/VariExtEvent.

Supplementary Material

S1 Derivation of the Euler–Lagrange equation

In this section, we detail the derivation of the Euler–Lagrange equations. We first form the constrained Lagrangian functional,

[TABLE]

where the function $\alpha:\Omega\to\mathbb{R}$ and the vector $\bm{\beta}=(\beta_{1},\cdots,\beta_{k})\in\mathbb{R}^{k}$ are the Lagrange multipliers. Taking the first variation of the constrained Lagrangian with respect to $\bm{u}$, we obtain

[TABLE]

Since the first variation $\delta\mathcal{L}_{c}/\delta\bm{u}$ must vanish for all $\bm{v}$, we obtain

[TABLE]

Similarly, the first variations of the Lagrangian $\mathcal{L}_{c}$ with respect to the Lagrange multipliers $\alpha$ and $\bm{\beta}$ read

[TABLE]

Since they must vanish for all $\tilde{\alpha}$, we obtain the constraints $\mathcal{K}(\bm{u})=0$ and $\bm{C}(\bm{u})=\bm{c}_{0}$.

S2 The Navier–Stokes equation

S2.1 Preliminaries

Recall the Navier–Stokes equations

[TABLE]

with the Kolmogorov forcing $\bm{f}(\bm{x})=\sin(k_{f}y)\bm{e}_{1}$ for some forcing wave number $\bm{k}_{f}=(0,k_{f})$. In two dimensions, a divergence free velocity field $\bm{u}:\Omega\to\mathbb{R}^{2}$ admits the following Fourier series expansion,

[TABLE]

where $\bm{k}=(k_{1},k_{2})$, $k=|\bm{k}|$ and $\hat{i}=\sqrt{-1}$ (see Ref.  (?)). Since the velocity field is real-valued, we have $a(-\bm{k})=-\overline{a(\bm{k})}$.

For the Kolmogorov forcing, the energy input rate satisfies

[TABLE]

where Im denotes the imaginary part and $a(\bm{k},t)=r(\bm{k},t)\exp(\hat{i}\phi(\bm{k},t))$ is the Fourier coefficient with phase $\phi(\bm{k},t)\in(-\pi,\pi]$ and amplitude $r(\bm{k},t)\in\mathbb{R}^{+}$. For simplicity, we may omit the dependence of these variables on time $t$. For reasons that will become clear in the next section, we refer to the Fourier mode $a(\bm{k}_{f},t)$ as the mean flow.

Examining equation (S20), the energy input $I$ may grow through two mechanisms:

1. (1)

The phase $\phi(\bm{k}_{f})$ approaching $-\pi/2$, 2. (2)

The amplitude $r(\bm{k}_{f})$ growing.

Noting that the phase of the external forcing is also $-\pi/2$, scenario (1) corresponds to an alignment between the phases of the external force and the mean flow $a(\bm{k}_{f})$. It is therefore tempting to attribute the intermittent bursts of the energy input $I$ to the intermittent alignments between the forcing $\bm{f}$ and the velocity field $\bm{u}$. This postulate, however, does not stand further scrutiny. Figure S6 shows the phase $\phi(\bm{k}_{f},t)$ of the mean flow along a typical Kolmogorov trajectory $\bm{u}(t)$. This phase oscillates around $-\pi/2$ for all times. Note that $-\pi/2$ corresponds to perfect alignment between the mean flow and external forcing. Figure S6 also shows the evolution of the energy input $I$ along the same trajectory. No positive correlation exists between intermittent growth of the mean flow energy $I$ and the phase of the mean flow being $-\pi/2$. In fact, the phase $\phi(\bm{k}_{f})$ seems to deviate from $-\pi/2$ during the bursts. Contrast this with the strong correlation between the growth of the energy input rate and the amplitude $r(\bm{k}_{f})$ of the mean flow.

This observation shows that the intermittent energy input bursts are triggered through mechanism (2), that is the growth of the amplitude $r(\bm{k}_{f})$. A similar observation is made at higher Reynolds numbers (not shown here). This growth of the mean flow amplitude, in turn, is possible through the internal transfer of energy via nonlinear terms as discussed below.

S2.2 Nonlinear triad interactions

The velocity field $\bm{u}(\bm{x},t)$ can be written in the general Fourier-type expansion

[TABLE]

where $\overline{\bm{u}}$ is the statistical mean and $\{\bm{v}_{j}\}$ is a set of prescribed functions that form a complete basis for the function space $X$ ($=L^{2}(\Omega)$). Under certain assumptions which are met by the Navier–Stokes equation (S18b), the energy is injected into the mean flow $\overline{\bm{u}}$ by the external forcing $\bm{f}$ (see Refs.  (?, ?)). The nonlinear term $\bm{u}\cdot\boldsymbol{\nabla}\bm{u}$, coupling the mean flow and the modes $\bm{v}_{j}$, redistributes the injected energy to all modes $\bm{v}_{j}$. This nonlinear term conserves the total energy of the system. At the same time, each mode dissipates energy due to the viscous term $\nu\Delta\bm{u}$ (see figure S7, for an illustration). A convenient choice of the basis $\{\bm{v}_{j}\}$ is problem dependent. Here, we choose the conventional Fourier basis as described in equation (S19). In case $\bm{f}$ is the Kolmogorov forcing, the symmetries of the system dictate $\overline{\bm{u}}(\bm{x},t)=\alpha_{0}(t)\bm{f}(\bm{x})=\alpha_{0}(t)\sin(k_{f}y)\bm{e}_{1}$ (see, e.g.,  (?, ?)).

In order to make the above statements more explicit, we write the Navier–Stokes equation in the Fourier space. Following  (?), we have

[TABLE]

where the hat signs denote the Fourier transform, $P_{ij}(\bm{k})=\delta_{ij}-k_{i}k_{j}/k^{2}$ is the Leray projection onto the space of diverge-free vector fields and the convention of summation over repeated indices is used. Equation (S22) can be written more explicitly as

[TABLE]

Recall the Fourier expansion (S19) which implies $\hat{u}_{1}(\bm{k})=k_{2}a(\bm{k})/k$ and $\hat{u}_{2}(\bm{k})=-k_{1}a(\bm{k})/k$. Upon substitution in equation (S23b) and noting that

[TABLE]

and

[TABLE]

we obtain

[TABLE]

We rewrite the above equation more compactly,

[TABLE]

where $\bm{k}\cdot\bm{q}=k_{1}q_{1}+k_{2}q_{2}$ and $\mu(\bm{p},\bm{q}):=p_{1}q_{2}-p_{2}q_{1}$ is the two-form measuring the surface area of the parallelogram with sides $\bm{p}$ and $\bm{q}$. Writing the modes in terms of their amplitudes and phases, $a(\bm{k})=r(\bm{k})\exp[\hat{i}\phi(\bm{k})]$, and using equation (S25), we obtain

[TABLE]

Note that $a(-\bm{k})=-\overline{a(\bm{k})}$ implies $\phi(-\bm{k})=\pi-\phi(\bm{k})$.

We now focus on the amplitude of the mean flow $r(\bm{k}_{f})$ (and its corresponding conjugate at $\bm{k}=-\bm{k}_{f}$). The negative definite term $-\nu k_{f}^{2}r(\bm{k}_{f})$ representing the dissipation acts to decrease the mean flow amplitude. This decay is counteracted by the external forcing $\frac{1}{2}\cos\left[\frac{\pi}{2}+\phi(\bm{k}_{f})\right]$. Recall from figure S6 that the phase $\phi(\bm{k}_{f})$ oscillates around $-\pi/2$ for all times, implying $\cos\left[\frac{\pi}{2}+\phi(\bm{k}_{f})\right]>0$. The complications arise from the summation term in (S26a) which couples the mean flow to all other modes that form the wave vector triads, $\bm{p}+\bm{q}=\bm{k}_{f}$. The contribution from these other modes depends on the amplitudes, $r(\bm{p})$ and $r(\bm{q})$, and the relative phases $\phi(\bm{k}_{f})-\phi(\bm{p})-\phi(\bm{q})$. Even the modes that do not form a triad with $\bm{k}_{f}$, affect the mean flow amplitude indirectly through their coupling to the modes that do form a triad with $\bm{k}_{f}$ (see the schematic figure S7).

S2.3 Derivation of Euler-Lagrange equation for Navier–Stokes

We first derive the functional $J$ corresponding to the Navier–Stokes equation and the energy input rate $I$. For the function space $X$ we set $X=L^{2}(\Omega)$ assuming that the state $\bm{u}$ belongs to the space of square integrable vector fields. By definition, we have $J(\bm{u})=dI(\bm{u};\mathcal{N}(\bm{u}))$ which implies

[TABLE]

where we used integration by parts. The term involving the pressure $p$ vanishes since the forcing is divergence free, $\boldsymbol{\nabla}\cdot\bm{f}=0$. Since $\|\bm{f}\|_{2}$ is constant, we can safely omit the second term and let

[TABLE]

Next we compute the Gâteaux differential of $J$. By definition, we have

[TABLE]

On the other hand, by Riesz representation theorem, we have $\mathrm{d}J(\bm{u};\bm{v})=\langle J^{\prime}(\bm{u}),\bm{v}\rangle_{L^{2}}$ which implies

[TABLE]

Similarly, the Gâteaux differential of the constraint $C(\bm{u})=\langle A\bm{u},A\bm{u}\rangle_{L^{2}}/(2|\Omega|)$ is given by

[TABLE]

implying $C^{\prime}(\bm{u})=A^{\dagger}A\bm{u}/|\Omega|$. Finally, the adjoint of the divergence operator, $\mathcal{K}=\boldsymbol{\nabla}\cdot$, with respect to the $L^{2}$ inner product is the gradient operator, $\mathcal{K}^{\dagger}=-\boldsymbol{\nabla}$. Substituting the above in the Euler–Lagrange equation (S16) and (S17), we obtain

[TABLE]

A few remarks about equations (S30c) are in order: (i) The PDE (S30a) is inhomogeneous due to the term $\nu\Delta\bm{f}=-\nu k_{f}^{2}\sin(k_{f}y)\bm{e}_{1}$. (ii) The equations are nonlinear in the constraint (S30c). (iii) With the Kolmogorov forcing $\bm{f}=\sin(k_{f}y)\bm{e}_{1}$, the translations $\bm{u}(x,y)\mapsto\bm{u}(x+\ell,y)$, with $\ell\in\mathbb{R}$, are a symmetry transformation of equations (S30c). That is, if $\bm{u}(x,y)$ solves (S30c), so does $\tilde{\bm{u}}(x,y)=\bm{u}(x+\ell,y)$ for all $\ell\in\mathbb{R}$.

S3 Newton iterations

In this section, we outline the Newton iterations for solving the system (S30c). Define

[TABLE]

The zeros of $\mathcal{F}$ coincide with the solutions of (S30c). We find these zeros numerically using damped Newton iterations

[TABLE]

At each iteration, the Newton direction $(\tilde{\bm{u}},\tilde{\alpha},\tilde{\beta})$ is obtained as the solution of the linear equation

[TABLE]

where $\mathcal{L}(\bm{u},\alpha,\beta;\cdot,\cdot,\cdot)$ is the Gateaux differential of $\mathcal{F}$ at $(\bm{u},\alpha,\beta)$ and is given explicitly as

[TABLE]

The solution of the linear PDE (S33) is approximated by the generalized minimal residual (GMRES) algorithm  (?). At each iteration, the step size $\epsilon\in(0,1]$ is adjusted to achieve maximal decrease in the error $\|\mathcal{F}(\bm{u}_{n+1},\alpha_{n+1},\beta_{n+1})\|_{L^{2}}$  (?). The standard Newton iterations correspond to $\epsilon=1$.

S4 Sensitivity to parameters

Recall that the constraint $\int_{\Omega}|\boldsymbol{\nabla}\bm{u}|^{2}\mathrm{d}\bm{x}/(2|\Omega|)=c_{0}$ enforces a constant energy dissipation rate. This constraint is motivated by the fact that, away from extreme bursts, the energy dissipation rate $D$ exhibits small oscillations around its mean value. Nonetheless, $D$ is not exactly constant, prompting the question whether the optimal solution is robust with respect to small perturbations to the constant $c_{0}$.

To examine this robustness, we have computed the optimal solution for a wide range of parameters $c_{0}$. We find that the optimal solution is in fact robust even with respect to relatively large variations in the parameter $c_{0}$. Figure S8, for instance, shows the optimal solution for three different values of $c_{0}$ at $Re=40$ and $100$ (the results are similar for $Re=60$ and $80$).

The insensitivity of the optimal solution with respect to the constant $c_{0}$ also implies that the equality constraint $\int_{\Omega}|\boldsymbol{\nabla}\bm{u}|^{2}\mathrm{d}\bm{x}/(2|\Omega|)=c_{0}$ can be replaced with an inequality constraint of the form $c_{1}\leq\int_{\Omega}|\boldsymbol{\nabla}\bm{u}|^{2}\mathrm{d}\bm{x}/(2|\Omega|)\leq c_{2}$. For a wide range of values for $c_{2}>c_{0}>c_{1}>0$, the optimal solutions corresponding to the two constraints will be similar.

S5 Computing the probability of extreme events

We approximate the conditional PDFs using the following steps. For any two observables $\lambda$ and $\gamma$, we assume that their joint probability density function $p_{\gamma,\lambda}$ exists such that

[TABLE]

Similarly, we also assume that the observable $\lambda$ has a probability density $p_{\lambda}$. Once the PDF $p_{\lambda}$ and the joint PDF $p_{\gamma,\lambda}$ are approximated using direct numerical simulations, the conditional PDF $p_{\gamma|\lambda}$ can be evaluated by the Bayesian formula,

[TABLE]

Computation of the extreme event probability $P_{ee}$ from the conditional probability is straightforward. Let $\gamma_{e}$ denote the threshold such that $\gamma>\gamma_{e}$ denotes an extreme event. Then by definition, we have

[TABLE]

where $\gamma^{\prime}$ is a dummy integration variable. In the present paper, the variable $\lambda$ is the indicator $|a(1,0)|$ and the variable $\gamma$ is the future maximum of the energy dissipation rate, $\gamma(t)=D_{m}(t)=\max_{\tau\in[t+t_{i},t+t_{f}]}D(\bm{u}(\tau))$. At each Reynolds number, the joint probability $p_{\gamma,\lambda}$ is approximated from the $100,000$ computed data points on a $20\times 30$ grid over the $(\gamma,\lambda)$ plane.

S6 Supporting computational results

In this section, we present the numerical results for Reynolds numbers $Re=40$, $60$, $80$ and $100$. The relevant parameters and variables are summarized in Table S1. At each Reynolds number, the statistics are computed from long trajectory data of length $10,000$ time units. The states (i.e. the velocity fields $\bm{u}$) are saved along these trajectories at every $0.1$ time units, amounting to a combined $100,000$ distinct states at each Reynolds number. Before recording any data, we evolved random initial conditions for $500$ time units to ensure the decay of transients.

The Navier–Stokes equations are solved numerically with a standard pseudo-spectral code with $N\times N$ Fourier modes and 2/3 dealiasing  (?) and a forth-order Runge–Kutta scheme for the temporal evolution. For $Re=60$, $80$ and $100$, we use $256\times 256$ Fourier modes to fully resolve the velocity fields. At $Re=40$, however, this resolution is unnecessarily high and hence we use $128\times 128$ modes.

Figure S9 shows the joint PDFs of the mode $a(1,0)$ versus the energy input $I$. At all Reynolds numbers the joint PDFs have a cone shape reflecting the fact that small values of $|a(1,0)|$ correspond to large values of the energy input rate.

As in $Re=40$, we use the evolution of $|a(1,0)|$ to predict an upcoming burst of the energy dissipation $D$. Figure S10 shows the computational results at higher Reynolds numbers. For $Re=60$, $80$ and $100$, we set the threshold $D_{e}$ for the extreme energy dissipation to be the mean plus one standard deviation of the energy dissipation. The measured mean $\mathbb{E}[D]$ and standard deviation $\sqrt{\mathbb{E}[D^{2}]-\mathbb{E}[D]^{2}}$ are reported in Table S1. The corresponding extreme dissipation thresholds $D_{e}$ are marked by vertical red dashed lines in the middle panel of figure S10. The horizontal dashed line marks the critical $\lambda_{0}$ for which $P_{ee}=0.5$, that is $50\%$ probability of an upcoming extreme event.

We recall from the main body of the paper that the four quadrants in the conditional PDFs (middle column of figure S10) correspond to:

1. (I)

Correct rejection ($P_{ee}<0.5$ and $D_{m}(t)<D_{e}$): Correct prediction of no upcoming extremes. 2. (II)

False positives ($P_{ee}>0.5$ but $D_{e}(t)<D_{e}$): The indicator predicts an upcoming extreme event but no extreme event actually takes place. 3. (III)

Hit ($P_{ee}>0.5$ and $D_{m}(t)>D_{e}$): Correct prediction of an upcoming extreme. 4. (IV)

False negatives ($P_{ee}<0.5$ but $D_{m}(t)>D_{e}$): An extreme event takes place but the indicator fails to predict it.

Table S1 also shows the results of the extreme event prediction. In order to quantify the success of these predictions, we define

[TABLE]

which measures the ratio of the number of extreme events that were successfully predicted to the total number of extreme events. Similarly, the quantity,

[TABLE]

measures the ratio of the number of non-extreme events that were correctly rejected to the total number of non-extreme events.