Abstract
Extreme events are ubiquitous in a wide range of dynamical systems, including turbulent fluid flows, nonlinear waves, large-scale networks, and biological systems. 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, that is, 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 SDs away from its mean, resulting in heavy tails of the corresponding probability distribution. Important examples include climate phenomena (1, 2), rogue waves in oceanic and optical systems (3–5), and large deviations in turbulent flows (6–8). Because these extreme phenomena typically have marked consequences, their quantification and prediction are 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 the work of Touchette (9) for a review]. Although 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 (10, 11). For systems that are not perturbations from these 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 (12–17).
Here, we propose a variational framework to probe the underlying conditions that lead to extreme events in these 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 lifetime of the extreme events are short compared to the typical dynamical time scales of the system, then 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 (18, 19). 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. However, the optimal solution to our variational method 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 use 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
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(u0) = c0 in Eq. 2B].
(A) A depiction of intermittent bursts of an observable. The highlighted regions mark an approximation of the growth phase of the extreme events. (B) In the state space, extreme events are viewed as fast excursions away from the background attractor (blue ball) owing to small regions of stability.
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 u0 whose associated observable I(u(t)) attains a maximal growth within time τ. More precisely, we seek the solutions to the constrained optimization problem
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 L2 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 counterpart
Problems 3A and 3B seek initial states u0 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 3B
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 u0.
If J : X → ℝ is a continuous map and the subset S = {u ∈ X : K(u) = 0, C(u) = c0} 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 equations
Here, K† is the adjoint of K, and J′(u) and Ci′(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 dCi(u;v) = 〈Ci′(u),v〉X for all v. The existence and uniqueness of J′(u) and Ci′(u) are guaranteed by the Riesz representation theorem (31). Here, Ci are the components of the map C = (C1, C2, ⋯, Ck). The function α : Ω → ℝ and the vector β = (β1, ⋯, β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
The kinetic energy E (per unit volume), the energy dissipation rate D, and the energy input rate I are defined, respectively, by
The Kolmogorov flow admits the laminar solution
(A) The time series of energy dissipation rate D at Reynolds number R = 40. (B) A close-up of the energy input I (solid red curves) and the energy dissipation D (dashed black curves) at Re = 40. The bursts in the energy dissipation are slightly preceded with a burst in the energy input. A similar behavior is observed for all bursts and at higher Reynolds numbers. (C) The vorticity field ∇ × u(x, t) = ω(x, t)e3 at time t = 433 over the domain x ∈ [0, 2π] × [0, 2π].
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(kf)| sin(φ(kf)), 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(kf) 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, φ(kf) → −π/2, and (ii) growth of the mean flow energy |a(kf)|.
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) = c0.
Constraints
The constraint C(u) = c0 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 form
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) = c0 = 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 reads
The associated Euler-Lagrange equations (Eqs. 6A to 6C) read
Using the symmetries of Eqs. 11A to 11C, we find that it admits the pair of exact solutions
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.
The solid (dashed) black line corresponds to the exact solution u– (u+). The red solid line (circles) corresponds to the global maximizer. The outset shows the scalar vorticity ∇ × u(x, y) = ω(x, y)e3 and the Fourier spectrum |a(kx, ky)| of the global maximizer at select Reynolds numbers.
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, kf) and the wave numbers (1, 0) and (1, kf) together with their complex conjugate pairs. Incidentally, these wave numbers form a triad, (0, kf) + (1, 0) = (1, kf). 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, kf), a(1, 0), and a(1, kf). We find that the energy transfers within this triad underpin the intermittent bursts of the mean flow energy |a(0, kf)| and, hence, the energy input rate I. Figure 4A shows the time series of I = −Im[a(0, kf)] 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).
(A) Time series of the energy input I and the modulus of the Fourier mode a(1, 0) at Re = 40. The eddy turnover time at this Reynolds number is te = 0.46. (B) The joint probability density of the energy input I and the real and imaginary parts of the mode a(1, 0), approximated from 100,000 samples. The density decreases from dark green to light blue. The cone-shaped density indicates the strong correlation between the large values of the energy input rate I and small values of |a(1, 0)|. The axisymmetric nature of the probability density is a consequence of the translation invariance of the Kolmogorov flow.
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, kf) through the triad interaction of the modes a(0, kf), a(1, 0), and a(1, kf). 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, kf). 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, kf), a(1, 0), and a(1, kf) 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 probability
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 Dm is most likely larger than 0.2. Conversely, when λ is larger than 0.4, the future values of the future energy dissipation Dm are smaller than 0.2.
(A) The probability density associated with the conditional probability (Eq. 12). The vertical dashed line marks the extreme event threshold
We seek an appropriate value λ0 such that λ(t) < λ0 predicts an extreme burst of energy dissipation over the future time interval [t + ti, t + tf]. Denote the extreme event threshold by De, such that D > De constitutes an extreme burst of energy dissipation. We define the probability of an upcoming extreme event Pee as
Figure 5B shows the probability of extreme events Pee as a function of the parameter λ0. If, at time t, the values of λ(t) are larger than 0.5, the probability of a future extreme event, that is, Dm(t) > De, 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 Pee < 0.5 to the higher likelihood Pee > 0.5. This line, together with the vertical line Dm = De, divides the conditional probability density into four regions: (I) correct rejections [Pee < 0.5 and Dm(t) < De]: correct prediction of no upcoming extremes; (II) false positives [Pee > 0.5 but De(t) < De]: the indicator predicts an upcoming extreme event but no extreme event actually takes place; (III) hits [Pee > 0.5 and Dm(t) > De]: correct prediction of an upcoming extreme event; and (IV) false negatives [Pee < 0.5 but Dm(t) > De]: 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 ti = 1 is approximately twice the eddy turnover time te. In comparison, it takes approximately one eddy turnover time (on average) for the energy dissipation rate to grow from
The prediction time ti can always be increased at the cost of increasing false positives and/or false negatives. For instance, with the choice ti = 2 ≃ 4.3te and tf = 3 ≃ ti + 2.2te, prediction of the extreme events Dm > De 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 ti = 2 still yields satisfactory predictions. Upon increasing ti 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 time scale of the extreme events is short compared to the typical time scales 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 can 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 were solved numerically with a standard pseudospectral code with N × N Fourier modes and 2/3 dealiasing. For Re = 60, 80, and 100, we used N = 256 to fully resolve the velocity fields. However, at Re = 40, this resolution was unnecessarily high, and hence, we used N = 128. The temporal integration of the Navier-Stokes equations was carried out with a fourth-order Runge-Kutta scheme.
SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/9/e1701533/DC1
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 versus 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.
movie S1. The prediction of an extreme event in the Kolmogorov flow.
This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.
REFERENCES AND NOTES
- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).