A variational approach to probing extreme events in turbulent dynamical systems

Variational formulation of extreme events

Consider the general evolution equations∂tu=N(u)(1A)K(u)=0(1B)u(⋅,t0)=u0(⋅)(1C)where u : Ω × ℝ⁺ → ℝ^d belongs to an appropriate function space X and completely determines the state of the system. The initial condition X ∋ u₀ : Ω → ℝ^d is specified at the time t₀, and Ω ⊂ℝ^d is an open bounded domain. The differential operators N and K can be potentially nonlinear. A wide range of physical models can be written as a set of partial differential equations (PDEs), as in Eqs. 1A to 1C. For instance, for incompressible fluid flows, Eq. 1A is the momentum equation and Eq. 1B is the incompressibility condition where K(u) = ∇ · u. For simplicity, we will denote a trajectory of Eqs. 1A to 1C by u(t).

Let I : X → ℝ denote an observable whose time series I(u(t)) along a typical trajectory u(t) exhibits intermittent bursts (see Fig. 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 Fig. 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 using a combination of observed data from the system and the governing equations of the system. We also require the instability regions to have a nonzero probability of occurrence under the natural flow of the dynamical system. This constraint is of particular importance because it excludes “exotic” states with extreme growth of I but with a negligible probability of being observed in practice [see the constraint C(u₀) = c₀ in Eq. 2B].

• Download high-res image
• Open in new tab
• Download Powerpoint

We formulate this task as a constrained optimization problem. Assume that there is a typical time scale τ ∈ ℝ⁺ over which the bursts in the observable I develop (see Fig. 1A). We seek initial conditions u₀ whose associated observable I(u(t)) attains a maximal growth within time τ. More precisely, we seek the solutions to the constrained optimization problemsupu0∈X(I(u(t0+τ))−I(u(t0)))(2A)where{u(t) satisfies (Eqs. 1A to 1C)C(u0)=c0(2B)where the optimization variable is the initial condition u(t₀) = u₀ of systems 1A to 1C. The set of critical states is required to satisfy the constraints in Eq. 2B to enforce two important properties. The first property ensures that u(t) obeys the governing Eqs. 1A to 1C as opposed to being an arbitrary one-parameter family of functions. The second property C(u₀) = c₀, where C : X → ℝ^k, is a codimension k constraint. This constraint is enforced to ensure the nonzero probability of occurrence, that is, 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 C(u₀) = c₀ is problem-dependent and will shortly be discussed in more detail. We point out that more general inequality constraints of the form c_min ≤ C(u₀) ≤ c_max may also be used. However, the treatment of these inequality constraints is not discussed in this paper.

We expect the set of solutions to problems 2A and 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 (that is, the instability regions of Fig. 1B) will result in a similar observable burst.

We emphasize that an optimization problem similar to Eqs. 2A and 2B has been pursued before in special contexts. The largest finite-time Lyapunov exponent can be formulated as Eqs. 2A and 2B, where the observable is the amplitude of infinitesimal perturbations after finite time. The maximizer is the corresponding finite-time Lyapunov vector (20, 21). In a similar context, Pringle and Kerswell (22) 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 Eqs. 2A and 2B, with the observable being the L² norm of the fluid velocity field. Ayala and Protas (23–25) 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 (Eqs. 2A and 2B) is an adjoint-based gradient iterative method (26–28). This method is computationally very expensive because, at each iteration, the gradient direction needs to be evaluated as the solution of an adjoint PDE. If the growth time scale τ is small compared to the typical time scales of the observable, then it is reasonable to replace the finite-time growth problem (Eqs. 2A and 2B) with its instantaneous counterpartsupu0∈Xddt|t=t0I(u(t))(3A)where{u(t) satisfies (Eqs. 1A to 1C)C(u0)=c0(3B)

Problems 3A and 3B seek initial states u₀ for which the instantaneous growth of the observable I along the corresponding solution u(t) is maximal.

We point out that the large instantaneous derivatives of I do not necessarily imply a subsequent burst in the observable because the growth may not always be sustained at later times along the trajectory u(t). As a result, the set of solutions to this instantaneous problem may differ significantly from the finite-time problem (Eqs. 2A and 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 Problems 3A and 3B. Taking the time derivative of the time series I(u(t)) yields (d/dt)I(u(t)) = dI(u; ∂_tu), where dI(u; v) := lim_ε→0 [I(u + εv) − I(u)]/ε denotes the Gâteaux differential of I at u evaluated along v. Using Eq. 1A to substitute for ∂_tu, we obtain the following optimization problem, which is equivalent to problems 3A and 3Bsupu∈XJ(u)(4A)subject to {K(u)=0C(u)=c0(4B)whereJ(u):=dI(u;N(u))(5)

Note that the first constraint in Eq. 3B is simplified because we have already used Eq. 1A, and it only remains to enforce Eq. 1B. For notational simplicity, we omit the subscript from u₀.

If J : X → ℝ is a continuous map and the subset S = {u ∈ X : K(u) = 0, C(u) = c₀} is compact in X, problems 4A and 4B have 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) ⊂ℝ is compact, which implies that J(S) is bounded and closed (29). 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 is a compact subset of S (30).

As we show in section S1, if X is a Hilbert space with the inner product 〈⋅,⋅〉_X and the operator K is linear, then every solution of the optimization problem (Eqs. 4A and 4B) satisfies the set of Euler-Lagrange equationsJ′(u)+K†(α)+∑i=1kβiCi′(u)=0(6A)K(u)=0(6B)C(u)=c0(6C)

Here, K^† is the adjoint of K, and J′(u) and C_i′(u) are the unique identifiers of the Gâteaux differentials dJ(u;⋅) and dCi(u;⋅), such that dJ(u;v) = 〈J′(u), v〉_X and dC_i(u;v) = 〈C_i′(u),v〉_X for all v. The existence and uniqueness of J′(u) and C_i′(u) are guaranteed by the Riesz representation theorem (31). Here, C_i are the components of the map C = (C₁, C₂, ⋯, C_k). The function α : Ω → ℝ and the vector β = (β₁, ⋯, β_k) ∈ ℝ^k are unknown Lagrange multipliers to be determined together with the optimal state u : Ω → ℝ^d.

Application to Navier-Stokes equations

We consider the Navier-Stokes equations∂tu=−u⋅∇u−∇p+νΔu+f,∇⋅u=0(7)where u : Ω × ℝ⁺ → ℝ^d is the fluid velocity field, p : Ω × ℝ⁺ → ℝ is the pressure field, and ν = Re⁻¹ is the nondimensional viscosity, which coincides with the reciprocal of the Reynolds number Re. Here, we consider two-dimensional flows (d = 2) over the domain Ω = [0, 2π] × [0, 2π] with periodic boundary conditions. The flow is driven by the monochromatic Kolmogorov forcing f(x) = sin(k_fy)e₁, where k_f = (0, k_f) is the forcing wave number and the vectors e_i denote the standard basis in ℝ^d. In the following, we assume that the velocity fields are square integrable for all times, that is, X = L² (Ω).

The kinetic energy E (per unit volume), the energy dissipation rate D, and the energy input rate I are defined, respectively, byE(u)=1|Ω|∫Ω|u|22dx, D(u)=ν|Ω|∫Ω|∇u|2dx, I(u)=1|Ω|∫Ωu⋅fdx(8)where |Ω| denotes the area of the domain, that is, |Ω| = (2π)². Along any trajectory u(t), these three quantities satisfy E.=I−D. We use the energy dissipation rate D to define the eddy turnover time, te=ν/E[D], where E denotes the expected value.

The Kolmogorov flow admits the laminar solution u=(Re/kf2) sin(kfy)e1. For the forcing wave number k_f = 1, the laminar solution is the global attractor of the system at any Reynolds number (32). If the forcing is applied at a higher wave number and the Reynolds number is sufficiently large, then 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 (18, 33, 34). This is manifested in Fig. 2A, showing the time series of the energy dissipation D at Reynolds number Re = 40 with k_f = 4.

• Download high-res image
• Open in new tab
• Download Powerpoint

A closer inspection reveals that each burst of the energy dissipation D is shortly preceded by a burst in the energy input I (see Fig. 2B). 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 section S2.1, the energy input is given by I(u) = −|a(k_f)| sin(φ(k_f)), where a(k) are the Fourier modes of the velocity field u and φ(k) are their corresponding phases such that a(k) = |a(k)| exp(iφ(k)). We refer to the Fourier mode a(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, that is, φ(k_f) → −π/2, and (ii) growth of the mean flow energy |a(k_f)|.

Examining the alignment between the forcing and the velocity field rules out mechanism (i) (cf. fig. S1). 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(k) through triad interactions of the modes whose wave numbers (k, p, q) satisfy k = p + q (see section S2.2 for further details). Because of the high number of active modes involved in the intricate network of triad interactions, it is unclear which triad (or triads) is (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. However, before obtaining the optimal solution, we need to specify the explicit form of the constraint C(u) = c₀.

Constraints

The constraint C(u) = c₀ is imposed to ensure that the optimal solutions are physically admissible, that is, that they are sufficiently close to the attractor and thus have a nonzero 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 (35). The maximizers of the functionals (Eq. 2A or Eq. 3A) that are far from this asymptotic attractor are physically irrelevant because 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 needs to be ensured otherwise. Here, we consider constraints of the formC(u):=1|Ω|∫Ω|A(u)|22dx(9)where A is a linear operator. Several physically important quantities can be written as function 9. For instance, if A is the identity operator, then C coincides with the kinetic energy, C(u) = E(u). If A is the gradient operator, A = ∇, then we have C(u) = D(u) × (Re/2), which can be used to constrain the energy dissipation rate. A more general class of these operators can be constructed as follows. Let u = ∑_iα_iv_i, where {v_i} is a principal component basis, that is, it diagonalizes the covariance operator of u. Define A as the diagonal linear operator such that A(v_i) = v_i/σ_i, where σi2 is the SD of α_i. Then, the constraint takes the form C(u)=(∑iαi2/σi2)/(2|Ω|). 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 (36). Note that the constraints of Eq. 9 are codimension one, C : X → ℝ, and hence a special case of the codimension k constraint in Eq. 2B.

Excluding the intermittent bursts, the energy dissipation of the Kolmogorov flow exhibits small oscillations around its mean value E[D] (see Fig. 2A). On the basis of this observation, we seek optimal solutions of Eqs. 4A and 4B, which are constrained to have the energy dissipation D = E[D]. This results into the constraint (Eq. 9) with A = ∇ and C(u) = c₀ = E[D] × (Re/2). We approximate the mean value E[D] from direct numerical simulations. At Re = 40, for instance, we have E[D] ≃ 0.117.

Probing the extreme energy transfers

The functional J (see Eq. 5) associated with the energy input I readsJ(u)=1|Ω|∫Ω[u⋅(u⋅∇f)+νu⋅(Δf)]dx(10)

The associated Euler-Lagrange equations (Eqs. 6A to 6C) read(∇f+∇fT)u+νΔf−∇α+βA†Au=0(11A)∇⋅u=0(11B)1|Ω|∫Ω|A(u)|22dx=c0(11C)where J′(u) = (∇f + ∇f^T) u + νΔf, K^†(α) = −∇α, and C′(u) = A^† Au (see section S2.3 for the derivations). We set A = ∇ to enforce a constant energy dissipation constraint. This implies that A^†Au = −Δu.

Using the symmetries of Eqs. 11A to 11C, we find that it admits the pair of exact solutions u±=±(2c0/kf) sin(kfy)e1, α±=±c0∫sin(2kfy)dy, and β±=±(νkf/2c0). More complex solutions, with unknown closed forms, may exist. We approximate these solutions using the Newton iterations described in section S3.

At each Re, we initiated several Newton iterations from random initial conditions. In addition to the pair of exact solutions (u_±, α_±, and β_±), the iterations yielded one nontrivial solution. Figure 3 shows the resulting three branches of solutions including the exact solution u₊ (solid black), the exact solution u₋ (dashed black), and the nontrivial solution (red circles). For small Reynolds numbers, our Newton searches only returned the exact solutions. At Re ≃ 3.1, a bifurcation takes place where a new nontrivial solution is born. This solution appears to be a global maximizer because no other solutions were found. Because the intermittent bursts are only observed for Re > 35, we focus the following analysis on this range of Reynolds numbers.

• Download high-res image
• Open in new tab
• Download Powerpoint

The nontrivial optimal solution converges to an asymptotic limit as the Re increases. This is discernible from the plateau of the red curve in Fig. 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 nonzero 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 4A shows the time series of I = −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 section S6).

• Download high-res image
• Open in new tab
• Download Powerpoint

This observation reveals that, before 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 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), and a(1, k_f) and their complex conjugates. Unfortunately, these severe truncations fail to illuminate because the dynamics of the truncated system severely departs from the original Navier-Stokes equations (37, 38).

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 λ(t) := |a(1, 0, t)|, along a solution u(t), signal the high probability of an upcoming burst in the energy dissipation.

To quantify this, we consider the conditional probabilityP(D1≤Dm(t)≤D2|λ1≤λ(t)≤λ2)(12)where Dm(t)=maxτ∈[t+ti,t+tf]D(u(τ)) 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₁, D₂], given that the present value of the indicator λ(t) belongs to [λ₁, λ₂]. 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 ≃ 2.2t_e and t_f = 2 ≃ 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 5A shows the conditional probability density corresponding to Eq. 12. We observe that relatively small values of λ correlate strongly with the high future values of the energy dissipation D. For instance, when λ < 0.4, the value of D_m is most likely larger than 0.2. Conversely, when λ is larger than 0.4, the future values of the future energy dissipation D_m are smaller than 0.2.

• Download high-res image
• Open in new tab
• Download Powerpoint

We seek an appropriate value λ₀ such that λ(t) < λ₀ 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 asPee(λ0)=P(Dm(t)>De|λ(t)=λ0)(13)which measures the likelihood that D_m(t) > D_e assuming that λ(t) = λ₀. Here, we set the threshold of the extreme event D_e as the mean value of the energy dissipation plus 2 SDs, De=E[D]+2E[D2]−E[D]2 ≃ 0.194. The extreme event probability P_ee can be computed from the probability density shown in Fig. 5A (see section S5 for details).

Figure 5B shows the probability of extreme events P_ee as a function of the parameter λ₀. If, at time t, the values of λ(t) are larger than 0.5, the probability of a future extreme event, that is, D_m(t) > D_e, is nearly zero. The probability of a future extreme event increases as λ(t) decreases. At λ(t) ≃ 0.4, the probability is 50%. If λ(t) < 0.3, then the likelihood of an upcoming extreme event is nearly 100%. The horizontal dashed line in Fig. 5A 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, divides 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; and (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 was 5.6%, and the number of correct rejections amounts to 93.3% of all predictions made. As we show in the Supplementary Materials, this amounts to a 95.6% success rate for the prediction of the extreme events (see eq. S24 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 turnover time t_e. In comparison, it takes approximately one eddy turnover time (on average) for the energy dissipation rate to grow from E[D]+E[D2]−E[D]2 to its extreme value De=E[D]+2E[D2]−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 ≃ 4.3t_e and t_f = 3 ≃ 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 eq. S24). 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.