Are extreme dissipation events predictable in turbulent fluid flows?
Patrick J. Blonigan Mohammad Farazmand Themistoklis P. Sapsis



Traces of exomoons in 

      computed flux and polarization phase curves 
 
                               of starlight reflected by exoplanets
    J. Berzosa MolinaNow at GMV AD, Calle de Isaac Newton, 11, 
        28760 Tres Cantos, Madrid, Spain 
        L. Rossi D. M. Stam
    Received 30 April 2018; accepted 23 July, 2018
================================================================================================================================================


We derive precursors of extreme dissipation events in a turbulent channel flow. Using a recently developed method that combines dynamics and statistics for the underlying attractor, we extract a characteristic state that precedes laminarization events that subsequently lead to extreme dissipation episodes. Our approach utilizes coarse statistical information for the turbulent attractor, in the form of second order statistics, to identify high-likelihood regions in the state space. We then search within this high probability manifold  for the state that leads to the most finite-time growth of the flow kinetic energy. This state has both high probability of occurrence and leads to extreme values of dissipation. We use the alignment between a given turbulent state and this critical state as a precursor for extreme events and demonstrate its favorable properties for prediction of extreme dissipation events. Finally, we analyze the physical relevance of the derived precursor and show its robust character for different Reynolds numbers. Overall, we find that our choice of precursor works well at the Reynolds number it is computed at and at higher Reynolds number flows with similar extreme events. 




§ INTRODUCTION

Turbulent fluid flows have been the most challenging paradigm of chaotic behavior with signatures of persistent and intermittent (i.e. over finite-times and at arbitrary time instants) instabilities leading to nonlinear energy transfers between scales. These nonlinear energy transfers are responsible for both the broad band character of the spectrum but also for the non-Gaussian statistics. More specifically, while non-zero third order statistics are primarily responsible for the persistent non-linear energy transfers (turbulent cascades) and the shape of the spectrum <cit.>, intermittent events, such as dissipation bursts, are primarily responsible for the heavy tail characteristics <cit.>. 

Here we are interested in the formulation of precursors for predicting these extreme events. These are  important in problems related to atmospheric and climate science, fluid-structure interactions, and turbulent fluid flow control, just to mention a few. We present our analysis on a standard configuration of a turbulent fluid flow, namely the channel flow, that exhibits extreme events in the form of large dissipation episodes occurring in random times <cit.>. These extreme events rise out of a high-dimensional turbulent attractor essentially without any clear warning. They have the form of a short-time excursion towards laminarization of the flow and a subsequent burst of turbulent kinetic energy which leads to a large dissipation episode pushing the flow away from the laminar state.

Many aspects of these intermittent bursts remain elusive primarily because of the intrinsic high dimensionality of the underlying turbulent attractor that limits the applicability of existing mathematical approaches <cit.>. In particular, extreme events due to their rare character cannot be `seen' effectively by energy-based methods, such as Proper Orthogonal Decomposition (POD). Even if one tries to consider conditional POD modes during extreme events these will not necessarily give the modes related to the triggering of the extreme events, as these do not necessarily obtain high energy, even during the extreme event. 

A different class of methods focus on the spectral analysis of the underlying Koopman operator <cit.>, and strive 
to extract unstable modes associated with certain observables of the system. Such unstable modes are typically estimated
by Dynamic Mode Decomposition (DMD) <cit.>.
This analysis, however, can only detect modes associated with long-term instabilities which
do not seem to explain short-term intermittent events 
observed in turbulent flows <cit.>.  Other variants, however, such as 
multi-resolution DMD <cit.> have been demonstrated to work well in systems 
with relatively low-dimensional attractors. 

Extreme events in complex dynamical systems have also been analyzed recently using Large Deviation Theory (LDT), e.g. in nonlinear water waves <cit.>. The basic idea is to search the phase space for initial conditions associated with a given magnitude of the objective function (observable of interest) and then from those pick the one with the highest probability of occurrence. However, these efforts have shown to work well in systems where the core of the attractor has Gaussian statistics. For different cases there is no rigorous foundation for LDT to operate. Even in the Gaussian case, the resulting optimization problem has very high dimensionality to be practically solvable for an application like the one considered here  – see <cit.> for a detailed discussion.

Here we apply a recently developed framework for the discovery of precursors to extreme events <cit.>. In contrast to LDT, in the present approach we consider a set of high probability initial conditions, based on a rough approximation of the attractor, and then search within this set for the initial state that has the highest growth for the quantity of interest. An adjoint solver is employed to compute the gradient of the objective function. Because the search is constrained within a low-dimensional (but high-probability) set, the feasibility set, the resulted optimization problem is computationally tractable. The solution provides with an initial state that has high probability to occur and leads to rapid growth of the objective function (in our case the kinetic energy). We discuss the physical relevance of the derived critical state in the context of the turbulent channel flow and use the alignment with this critical state as a scalar precursor for the prediction of future extreme events. We measure the effectiveness of the precursor through direct numerical experiments and examine its robustness over different Reynolds numbers. The success of our approach to an intermittently turbulent channel paves the way for studying transitional flows, such as bypass transition of boundary layers.

This paper is organized as follows. In Section <ref> we described the minimal channel flow used in this work
and discuss various aspects of the problem. The optimization problem for discovering the
precursor to extreme events is presented in Section <ref>. In Section <ref>
we present a statistical analysis quantifying the predictive power of the optimal precursor. 
Finally, our concluding remarks are presented in Section <ref>.



§ TURBULENT CHANNEL FLOW - PHENOMENOLOGY








 §.§ The minimal flow unit for channel flows


Turbulent channel flow has been a staple of numerical studies of turbulence for many years <cit.>. The chaotic nature of these simulations makes it difficult to analyze local spatiotemporal events and physical mechanisms in them, such as the formation and destruction of individual hairpin vortices in the near-wall region. To isolate these physical mechanisms and others, work has been done to find “minimal flow units” for various regions of the channel. Jimenez and Moin <cit.> found the minimal flow unit for near wall turbulence for several low Reynolds number flows by considering turbulent channel flow simulations with domains that were considerably smaller than conventional channel flow simulations. These smaller domains eliminate larger scale turbulent structures but accurately resolved the near-wall turbulent flow, matching turbulent flow statistics from experiments and prior numerical studies up to 40 wall units in the wall-normal direction. 

Various minimal flow units have been used in a range of different studies because of its isolation of a few physical mechanisms and its relatively low computational cost. Carlson and Lumley used the minimal flow unit to study flow control strategies for turbulent boundary layers <cit.>. Minimal flow units have also been used to study near-wall and log-layer turbulence <cit.>. These studies use forcing functions to achieve turbulent flows in half-channels, and to selectively damp larger scale flow structures. Recently, near-wall minimal flow units have been used to study and characterize the effects of wall-roughness on wall-bounded flows and to build models of wall roughness effects <cit.>. The near-wall minimal flow unit has also been used to demonstrate shadowing-based adjoint sensitivity analysis <cit.>. 

In addition, the near-wall minimal flow unit simulations routinely show highly intermittent behavior which is of interest for the purposes of our study. In certain low Reynolds number simulations, the flow on one wall exhibits turbulent behavior while the other remains laminar. The flow on both walls transitions at seemingly random intervals: the laminar wall would become turbulent, and then the turbulent wall would become laminar. Turbulent flow in the near-wall minimal flow unit is itself intermittent when it occurs. Time series of wall shear stress show that turbulence undergoes a cycle where it proliferates rapidly or “bursts”, then decays slowly. This observation led to numerous subsequent studies into the intermittent nature of near-wall turbulence using minimal flow units (see Refs <cit.>, for comprehensive reviews). 





 §.§ Flow Solver


In this study we consider near-wall minimal flow units similar to those considered in Ref. <cit.>. We use a Discontinuous-Galerkin spectral-element method (DGSEM) framework to simulate the minimal flow unit with the compressible Navier-Stokes equations <cit.>. The DGSEM framework has been successfully applied to a range of different flows including channels flows and the near-wall minimal flow unit  <cit.>. Also, it has an adjoint capability <cit.> that has been validated for the near-wall minimal flow unit in Ref. <cit.>. A detailed description of the discretization and implementation of this solver is available in Refs. <cit.> and <cit.>. 




























































We run a direct numerical simulation (DNS) of the channel flow with the compressible Navier-Stokes equations with a constant forcing in the axial direction to drive the channel. Since the Mach number is low, the effective governing equations reduce to the incompressible 
Navier–Stokes equations


    · u=0,


    ∂_t  u + u· u = f_0/ρ e_1 -1/ρ p +νΔ u,


    u( x,0)= u_0( x)


where u = (u,v,w) denotes the three-dimensional velocity field with streamwise component u, 
wall-normal component v and spanwise component w. The constant forcing in the
streamwise direction is denoted by f_0 and e_1=(1,0,0) is the unit vector in the streamwise direction. The boundary conditions are no-slip at the walls so that u(x,±δ,z,t)=0, 
and periodic in the spanwise and streamwise directions. Here, the channel height in the wall-normal direction is 2δ.



 §.§ Numerical experiment set-up


The case we consider has a domain size of πδ× 2δ× 0.34πδ in the streamwise, wall-normal, and spanwise directions, respectively.  The channel half-height δ was set to 1.0. The flow considered is at Reynolds number Re=2200 (corresponding to the friction Reynolds number Re_τ=110) unless stated otherwise. The Reynolds number Re is defined as
Re = ρ U δ/μ
where ρ is the fluid density, U is the centerline velocity of a laminar flow with the same mass flow rate (as in Ref. <cit.>), and μ is the dynamic viscosity. As in Ref. <cit.>, U was chosen so that the Mach number of the flow was under 0.3 and therefore the flow is nearly incompressible. The friction Reynolds number Re_τ is defined as 
Re_τ = ρ u_τδ/μ
where u_τ=√(τ_w/ρ) is the friction velocity and τ_w is the average shear stress at the wall. 
The channel constant forcing f_0 is set to balance the mean shear stress of the walls, so it is set by the choice of Reynolds number Re, and the friction Reynolds number Re_τ as follows

    f_0 = τ_w/δ = Re_τ^2/Reρ U


The domain is discretized with a 4 × 16 × 2 mesh with 8th order spatial elements, resulting in a 32 × 128 × 16 distribution of nodes (65536 total), similar to the mesh used in Ref. <cit.>. The choice of Re_τ=110.0 results in grid resolutions of Δ x^+ ≈ 11 and Δ z^+ ≈ 7 wall units per node, where
x^+ = u_τ x/ν , y^+ = u_τ y/ν, and z^+ = u_τ z/ν.

The wall-normal spacing corresponds to Δ y^+ ≈ 0.6 for the nodes closest to the walls, which ensures that the simulations are well resolved. We used space-time elements that were 4th order in time and the time slab (temporal element) was Δ t =0.05t_e, where
t_e = δ/U denotes the eddy turnover time, the time scale associated with the largest possible eddy in the channel.












The DGSEM flow solver used here has been validated for the minimal flow unit at Re=3000 <cit.> and we carried out a similar validation study for the Re=2200 case. Note that the statistics at Re=2200 were only computed over time intervals when both walls had turbulent flow, as was done in Ref. <cit.> for low Reynolds number cases. This was necessary because the Re=2200 flow exhibited intermittent behavior similar to that observed in previous studies of minimal flow units. 



 §.§ Extreme events




The intermittent behavior of the flow at Re=2200 can be seen in Fig. <ref> from the spikes in kinetic energy E(t) and dissipation Z(t). We define the kinetic energy,

    E(t) = ∭_Ωρ u · u  𝕀 x  𝕀 y  𝕀 z,

where Ω is the flow domain. Energy dissipation rate Z(t) is defined as

    Z(t) = ∭_Ωtr(τ∇ u )  𝕀 x  𝕀 y  𝕀 z,

where τ denotes the stress tensor, defined as 
τ = μ(∇ u + ∇ u^⊤)
for an incompressible flow.

Fig. <ref> shows that large increases in E(t) are followed by spikes in Z(t) and a subsequent decrease in E(t). These large spikes in kinetic energy occur when there is laminar flow near one wall, and nearly laminar flow near the other wall. Laminar flows correspond to higher kinetic energy E(t) because the channel is driven by a constant body force in the axial direction. This body force acts as a fixed axial pressure gradient. The body force and wall shear stress balance one another when the flow is in an equilibrium state where we have

    μ∂u/∂ n  = δ∂p/∂ x,

where n is the wall-normal direction. For a given centerline velocity, a laminar flow will exert less shear stress μ∂u/∂ n on the walls than a turbulent flow, so for a given wall shear stress μ∂u/∂ n the laminar flow will have a larger centerline velocity than a turbulent flow. 














Therefore, the large spikes in E(t) are the result of a flow laminarization event. Of course, the flow never completely reaches the laminar
state, though it gets very close it.
Figs. <ref> and <ref> show an example of a typical laminarization event where
the flow undergoes the following stages.


  * The flow on the bottom wall becomes laminar (Fig. <ref>). 

  * The flow on the top wall becomes nearly laminar (Fig. <ref>).

  * As the entire channel becomes nearly laminar, the axial velocity increases. 

  * The higher velocities make the effective Reynolds number of the flow larger. This increases the likelihood of turbulent burst occuring since the flow is less stable to perturbations at a larger Reynolds number. 

  * A turbulent burst occurs on the top wall, which causes Z(t) to increase rapidly (Fig. <ref>).

  * The bottom wall transitions to turbulence (Fig. <ref>).

  * The turbulent flow on both walls causes E(t) to decrease as it returns to the equilibrium mean turbulent flow profile (Fig. <ref>). 


These flow laminarization events, and the resulting bursts of energy 
dissipation rate, are the extreme events we will consider in this paper. 




§ OPTIMAL STATES FOR EXTREME EVENTS

We now describe the constrained optimization
problem whose solutions determine the most likely 
triggers of extreme events. The method is presented in detail in 
Ref. <cit.> and is reviewed here for 
completeness. We describe the optimization problem in the context of the channel flow, 
outline the numerical method for obtaining its solutions and present our numerical 
results.  



 §.§ High-Likelihood triggers of extreme events


We seek initial states u_0 that after a given integration time τ
realize the largest possible energy growth. More precisely, 
we seek initial states u_0 such that E( u(τ)) - E( u_0)
is maximized.
This is a PDE-constrained optimization problem,
since the velocity field u(t) is required to satisfy the channel flow (<ref>).

In addition to this PDE constraint, we also enforce a feasibility constraint, by requiring the initial state u_0
to belong to the system attractor. This second constraint is essential in order to guarantee that 
the optimal solution is probabilistically relevant. As in many dissipative PDEs, the channel flow has an attractor towards which the solutions converge asymptotically after an initial transient. We are interested in the self-sustained and recurrent extreme events on this attractor as opposed
to transient extreme events off the attractor that may occur but are not recurrent. In order to prevent the optimizer from considering such transient 
events, we enforce a feasibility constraint which is further elucidated in section <ref>.

With this prelude, the optimization problem can be formulated as
 

    max_ u_0∈𝒰[E( u(τ))-E( u_0) ],


    u(τ)   u(0)= u_0,


    u_0∈𝒜⊂𝒰,


where τ>0 is a prescribed time, related to the growth time scale of a typical extreme event. The constraint (<ref>) implies that u(τ)
is a solution of the channel flow. Constraint (<ref>) ensures that the 
optimizer belongs to the attractor 𝒜 and is therefore probabilistically relevant. 
In the next section, we describe the method used to approximate the attractor 𝒜. 




 §.§ Feasibility constraint and proper orthogonal decomposition


The attractors of dissipative dynamical systems are often very complex sets.
Their estimation has been the subject of many studies (see e.g.  <cit.>). Here, we approximate the attractors
via Proper Orthogonal Decomposition (POD) of longterm simulations of the channel flow
(this method is also known as the Principal Component Analysis). 
The POD approximation assumes that the attractor has a Gaussian
distribution with mean u( x) and covariance matrix C( x, x')
where 


    u( x) = lim_T→∞1/T∫_0^T u( x,t)𝕀 t,


    C( x, x')=lim_T→∞1/T∫_0^T 
    ( u( x,t)- u( x))⊗( u( x',t)- u( x'))𝕀 t.


Let the vector fields v_i:Ω→ℝ^3 denote the eigenfunctions of the
covariance tensor, so that

    ∫_Ω C( x, x') v_i( x')𝕀 x' = λ_i  v_i( x),   i∈ℕ

where λ_i∈ℝ are the corresponding eigenvalues. The eigenvectors are
ordered such that λ_1≥λ_2≥⋯. 
Since the 
covariance tensor is symmetric and positive definite, the eigenvalues are real-valued and non-negative, and furthermore the eigenfunctions
are orthogonal with respect to the L^2 inner product, 
i.e., ⟨ v_i, v_j⟩_L^2(Ω)=δ_ij.
We refer to the eigenfunctions v_i as the POD modes.

In the POD approximation, any state on the attractor is approximated as 

    u( x,t) =  u( x) +∑_i=1^nξ_i(t) v_i( x),

which is a finite-dimensional truncation to the first n POD modes. Each component of the
vector ξ = (ξ_1,⋯,ξ_n)∈ℝ^n is given by 
ξ_i(t) =⟨ u(t)- u, v_i⟩_L^2(Ω).

With this POD approximation of the attractor, the optimization problem (<ref>)
can be written more explicitly as
 

    max_ξ∈ℝ^n[E( u(τ))-E( u_0) ],


    u( x,τ)   u( x,0)= u_0( x),


    u_0( x)=  u( x) +∑_i=1^nξ_i v_i( x),


    ∑_i=1^nξ_i^2/λ_i≤ r_0^2


where r_0∈ℝ is a prescribed parameter that is set equal to 1.0 in this study. Note that the form of the constraint essentially restricts our optimization within states that have the highest 
probability, given second-order statistics for the attractor. Constraint (<ref>) enforces that 
the mean-zero initial condition u_0- u belongs to the subspace spanned by the first n POD modes. 
Constraint (<ref>), which describes an ellipsoid in the subspace 
{ v_1,⋯, v_n}, ensures that the initial conditions are not two far from the 
mean flow u (see Fig. <ref>).

Although the initial condition u_0 is 
constrained to the subspace spanned by the first n POD modes, the final state u(τ)
may not belong to this subspace. This is because the POD decomposition is only an approximation 
of the attractor, which represents initial states for our analysis, and not the exact invariant attractor. However, we take n large enough
so that only an insignificant fraction of the energy content of the states on the attractor are neglected. 
More precisely, in the following, we set n=50 so that the truncation of the turbulent states to 
these POD modes contain at least 90% of the kinetic energy, as shown in Fig. <ref>. 


As discussed in section <ref>, previous studies of the minimal flow unit have shown that intermittent bursts of the flow
originate from the near-wall activities. In light of these observations,
we make an important modification to the computation of the POD modes by 
multiplying the zero-mean velocity fields by a weight function that emphasizes the 
near-wall contribution of the flow. 
More precisely, we compute the weighted velocity fields

    u_ϵ,h( x,t)=( u( x,t)- u( x))q_ϵ,h(y),
 
by multiply the original velocity fields u (after removing the mean u) with the weight function 







    q_ϵ,h(y)= 1/2{
    2+tanh[1/ϵ(y-(δ-h)) ]  
    -tanh[1/ϵ(y+(δ-h)) ]
     }.

For ϵ≪ h ≪δ, the weight function q_ϵ,h vanishes near the 
center of the channel and approaches unity near the walls at y=± h (see figure <ref>).
The parameter h is the width of the near-wall region that we 
would like to emphasize and the small parameter ϵ determines how quickly the function 
q_ϵ,h decays to zero far from the walls.

 

The POD modes are computed through the equations (<ref>) and (<ref>)
except that in computing the covariance (<ref>) instead of the terms u( x,t)- u( x), we use the mean-zero weighted velocity u_ϵ,h( x,t).
In the following, we set ϵ = 0.1 to keep the decay smooth. We selected h = 25δ/110 to emphasize the 
near-wall flow up to ignore any flow above y^+≈50, where the minimal flow unit does 
not capture all length scales. 
This modified POD not only emphasizes the near-wall contributions, but also speeds up the convergence of the 
numerical optimization of problem (<ref>).

The POD modes were computed using 1000 snapshots taken from flows computed from 10 different randomly chosen, initial conditions. The snapshots were taken at intervals of 50 eddy turnover times to minimize any correlation in time. Fig. <ref> shows the energy and cumulative energy of the POD modes. Almost 80% of the energy is captured by the first six modes, indicating that near-wall dynamics emphasized by our weight function are low dimensional. 



Fig. <ref> shows a few of these modes. The shape of the most energetic modes is dominated by long axial streaks which are know to be a main feature of the near-wall region <cit.>. The more energetic modes such as modes 1 and 3 contain wider streaks, while less energetic modes such as mode 11 contain thinner, less coherent streaks that meander slightly. 

Finally, we point out that an energy maximizing optimization similar to (<ref>) was used by Pringle and Kerswell <cit.> in 
the context of transition to turbulence in the pipe flow. The optimal perturbations to the laminar state that
lead to turbulence were found as solutions to an optimization problem. 
In that context, the perturbations to a known state (i.e., the laminar state) 
are sought and therefore the knowledge of the turbulent attractor 
does not enter the optimization problem. In contrast, our study demands that the initial states belong to 
the turbulent attractor as enforced through the constraint (<ref>) and implemented numerically in this section.
Similarly, Farano et al. <cit.> investigated the observed bursts in a turbulent flow similar to our channel flow. However, they only
require that the energy of the optimal state is prescribed, i.e., E( u_0)=E_0 for a prescribed energy level E_0.
This does not necessarily imply that the optimal state belongs to the attractor and therefore unphysical 
optimizers are not ruled out.



 §.§ Numerical implementation


The optimization problem, equation (<ref>), was solved using the Python package scipy <cit.>. Specifically, the “optimize” package was used and the optimization was carried out using Sequential Least SQuares Programming (SLSQP). The first constraint, (<ref>), is implicitly satisfied by the flow solver, which is called from Python using the “multiprocessing” module. The second and third constraints (<ref>)–(<ref>) and their gradients are implemented directly in Python. The convergence tolerance for the objective function was set to 10^-7. Otherwise the default convergence criterion were used. An adjoint solver was used to minimize the cost of computing the gradient of the objective function (<ref>). The DGSEM solver has a dual consistent, discrete adjoint formulation, details of which are discussed in Refs. <cit.>. This adjoint solver allows us to compute the gradient at a computational cost similar to solving the primal, which is much cheaper than using finite differences to compute the gradient with respect to all n=50 POD modes. 



 §.§ Optimization results


The optimization was run from five different initial guesses u_0 and with three different integration times of τ=52.5t_e, τ=105t_e, and τ=210t_e. The integration times correspond to roughly 1/8, 1/4, and 1/2 of the time scale of laminarization event similar to the one shown in Fig. <ref>. All five optimizations with τ=52.5t_e, and the optimization with τ=105t_e computed very similar optimal solutions u_0^*. The optimization with an integration time of τ=210t_e failed to converge because the adjoint grew very large in magnitude, and the gradient caused the optimizer to consider a non-physical solution in the ensuing line search. 




Fig. <ref> shows the convergence history of a typical optimization for τ=52.5t_e. The corresponding POD mode weights for the initial guess and optimal solution are shown in  Fig. <ref>. The two largest POD modes in the optimal solution are modes 6 and 26, shown in Fig. <ref>. Together, these modes create the axial velocity deficit in the near wall region shown in Fig. <ref>. The maximum velocity deficit occurs roughly y^+=18 units from each wall, and decays to zero roughly y^+ ≈ 55 units from the wall. 







The region in which the velocity deficit occurs is known to have a major impact on near-wall turbulence. In <cit.> it was observed that damping axial velocity streaks or quasi-streamwise vortex structures between y^+ ≈ 20 to y^+ ≈ 60 led to laminarization of the near wall region. It appears that our optimal solution modifies the lower portion of this range. The relative uniformity of the axial velocity deficit results in the absence of axial velocity streaks and quasi-streamwise vortex structures. This is in contrast with the flow snapshots in Fig. <ref>, where the presence of low- and high- axial velocity streaks can be inferred from streaks in the wall shear stress, and quasi-axial structures are revealed by the q-criterion iso-contours. Without the axial velocity streaks and quasi-streamwise vorticity, the “streak-cycle” mechanism for the regeneration of near-wall turbulent fluctuations is broken and the flow laminarizes. 



Therefore, although the optimal solution u_0^* is turbulent, the absence of axial velocity streaks and quasi-streamwise vortices make it a precursor for a flow laminarization. By tracking how close a given state is to u_0^*, we can determine if a laminarization event is likely to occur or not.



§ PREDICTING EXTREME EVENTS

Figure <ref> shows a close-up of the evolution of 
energy E and dissipation Z together with the indicator 

    λ =⟨ u -  u , u_0^* -  u⟩/ u -  u_2  u_0^* -  u_2,

at Re=2200.
Large values of energy and dissipation are proceeded with relatively large values of the indicator. 
This behavior turns out to be generic and not specific to this time window. As a result, 
one can use the indicator λ to predict the upcoming extreme events in the channel flow.
In order to quantify these predictions, we first review some statistical tools in Section <ref>. 
Subsequently, in Section <ref>, we apply these tools to longterm simulations of the 
channel flow and report the results. 




 §.§ Statistical preliminaries

In this section, we show that the trigger state obtained previously can be used for the prediction of 
the extreme events in the channel flow. In order to make quantitative statements, we 
use joint and conditional statistics between the trigger mode and the energetic observables of the
turbulent flow, namely kinetic energy and energy dissipation rate.

For a given random variable X_t, we would like to find an indicator Y_t (another random variable)
whose values signal an upcoming extreme event of X_t.
We identify extreme events of X_t as any instant where X_t>x_e for a 
prescribed extreme value threshold x_e.
First we define the maximum value of the random 
variable X_t over a future time interval [t+t_0,t+t_0+Δ t]
for some t_0,Δ t≥ 0,

    X̃_t(t_0,Δ t) = max_s∈ [t+t_0,t+t_0+Δ t] X_s

The maximum X̃_t(t_0,Δ t) is a new random variable which depends on the
parameters t_0 and Δ t. At any time t, 
X̃_t(t_0,Δ t) measures the maximum value that X_t will take over 
the future time interval [t+t_0,t+t_0+Δ t].
For notational simplicity, we omit the parameters
t_0 and Δ t and simply write X̃_t. 

The joint probability distribution of the pair (X̃_t,Y_t) is defined by

    F_X̃_t,Y_t(x,y)   =ℙ(X̃_t≤ x,Y_t≤ y)
       =∫_-∞^x∫_-∞^yp_X̃_t,Y_t(x,y)𝕀 x𝕀 y,

where p_X̃_t,Y_t is the probability density associated with the probability distribution F_X̃_t,Y_t. The probability density p_X̃_t,Y_t(x,y)
measures the probability that at time t we observe Y_t=y and X̃_t=x.

The conditional probability of X̃_t=x given that Y_t=y is defined through the Bayes' formula

    p_X̃_t|Y_t = p_X̃_t,Y_t/p_Y_t,

where p_Y_t is the probability density associated with the random variable Y_t.
We use the conditional PDF p_X̃_t|Y_t to quantify the extent to which the
behavior of Y_t is indicative of the extreme events of X_t over the 
future time interval [t+t_0,t+t_0+Δ t].
More precisely, given an extreme event threshold x_e, we define the probability of upcoming extreme events by

    P_ee(y) = ∫_x_e^∞ p_X̃_t|Y_t(x,y)𝕀 x.

This quantity measures the probability of an extreme event over the future time interval 
[t+t_0,t+t_0+Δ t] given the current value of the indicator Y_t=y.

For a reliable indicator, P_ee should be monotonic so that the
probability of upcoming extreme events increases with y. 
More precisely, P_ee should be nearly zero for small values of y
and increase monotonically towards 1 as y increases. 
We predict an upcoming extreme only if P_ee>0.5. 
This defines an extreme event threshold y_e for the indicator Y_t
where P_ee(y_e)=0.5. 
If Y_t>y_e an upcoming extreme event is predicted and conversely if Y_t<y_e
it is predicted that no extreme events will occur over the future time interval [t+t_0,t+t_0+Δ t].

This classification leads to four possible prediction outcomes in terms of the indicator value Y_t and
the future observable values X̃_t:

    X̃_t<x_e      Y_t<y_e,
           X̃_t>x_e     Y_t>y_e,
           X̃_t>x_e     Y_t<y_e,
           X̃_t<x_e     Y_t>y_e.

Therefore, the skill of an indicator for predicting upcoming extreme events can be quantifies as follows


    = /+
        =∫_x_e^∞∫_y_e^∞ p_X̃_t|Y_t(x,y)𝕀 y𝕀 x/∫_x_e^∞∫_-∞^∞ p_X̃_t|Y_t(x,y)𝕀 y𝕀 x


    = /+
        =∫_-∞^x_e∫_-∞^y_e p_X̃_t|Y_t(x,y)𝕀 y𝕀 x/∫_-∞^x_e∫_-∞^∞ p_X̃_t|Y_t(x,y)𝕀 y𝕀 x



A skillful indicator is one that returns relatively small percentage of false negatives (respectively, false positives) 
compared to the number of correct predictions (respectively, correct rejections).
In the following, we use the statistical quantities introduced above to quantify the predictive skill of the 
indicator (<ref>).



 §.§ Prediction results

We first present the joint and conditional statistics for the energy E, 
dissipation rate Z and the indicator λ.
In this first step, we do not include any time shifts, thus setting
t_0=Δ t=0 so that X̃_t=X_t in equation (<ref>). 
Figure <ref> shows the joint and conditional PDFs of the
indicator versus energy E and the energy dissipation Z
at Re=2200. This figure is generated from an ensemble of longterm simulations with data recorded
every one eddy turnover time collecting a total of 106,063 data points.

The shape of the conditional PDFs p_E|λ and p_Z|λ
shows that the extreme values of the indicator λ 
correlate strongly with relatively large energy and large dissipation episodes (see Fig. <ref>(a,b)).
We also note that that the correlation is much stronger for dissipation and
that λ tends to increase monotonically with the dissipation Z. 



The next step is to introduce a time lag in order to investigate whether
the extreme indicator values proceed the extreme episodes of energy and dissipation. 
This is clearly the case in figure <ref>, however,
we show that this is true most of the time during the longterm simulations. 
Figure <ref> shows the conditional PDFs p_Ẽ|λ
and p_Z̃|λ where the future maxima Ẽ and Z̃ 
are computed with t_0=t_e and Δ t = 10t_e (see equation (<ref>)).
The extreme value threshold E_e (respectively, Z_e) are 
set as the mean of energy (respectively, dissipation) plus one standard deviation.
Figure <ref> also shows the 
corresponding P_ee computed from equation (<ref>). 
The extreme event threshold according to the indicator λ is 
the point at which P_ee(λ_e)=0.5.





As the prediction time t_0 increases, we expect the prediction skill of the indicator 
to deteriorate. This is shown in figure <ref> where the 
conditional PDF p_Z̃|λ are shown for increasing prediction times t_0. 
For prediction times as large as t_0=10t_e the prediction skill of the indicator 
is still reasonably satisfactory. However, as the prediction time increases to 
close t_0=25t_e the predictor returns significant amount of false positives and 
false negatives, thus losing its predictive value.
The rates of successful predictions and successful rejections, defined in equation (<ref>),
are reported in Table <ref> for a few prediction time horizons. 



Similar results are observed at higher Reynolds numbers. Figure <ref>
 shows the joint and conditional PDFs of the indicator, energy and dissipation at Re=3000.
These PDFs resemble those of figure <ref> for the lower Reynold number Re=2200. 
This demonstrates clearly the robustness of the derived indicator. 
Table <ref> also contains the rates of success in predicting extreme (and non-extreme) events
at Re=3000.
We point out that the channel flow at some intermediate Reynolds numbers 
between Re=2200 and Re=3000
did not exhibit extreme events in the time horizons we simulated. 





§ CONCLUSIONS


We have demonstrated an original approach for the derivation of precursors of extreme events in a challenging problem involving a turbulent channel flow. The extreme events in this case have the form of random near-laminarization episodes that lead to bursts of the kinetic energy and the energy dissipation rate. We formulate a constrained optimization problem that searches for initial states with the most intense growth of kinetic energy, within a constrained set in the core of the underlying turbulent attractor. By searching over a high probability set, we achieve a numerically tractable optimization problem, while at the same time we exclude exotic states that may correspond to intense growth of energy but have very low probability to occur if we are close to the attractor. 

The derived precursor is demonstrated to successfully capture extreme dissipation episodes several eddy turnover times before the event. We have discussed its physical relevance and have demonstrated its robustness as the Reynolds number of the flow changes. Because the developed scheme utilizes full simulations and not linearized approximations, it has the potential to be extended to larger, more complicated flows although divergence of the adjoint might be a limiting factor for obtaining the solution of the associated optimization problem. The success of the presented approach to an intermittently turbulent
channel flow implies the potential of the method for studying transitional flows, such as bypass transition of boundary layers. Our future endeavors include the utilization of these precursors for the control and suppression of extreme events in these systems.  



  §.§.§ Acknowledgments

PJB was supported by an appointment to the NASA Postdoctoral Program at the NASA Ames Research Center, administered by Universities Space Research Association under contract with NASA. This research was sponsored by NASA's Transformational Tools and Technologies (TTT) Project of the Transformative Aeronautics Concepts Program under the Aeronautics Research Mission Directorate. TPS and MF have been supported through the ARO MURI W911NF-17-1-0306 and the ONR grant N00014-15-1-2381.
Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.




10

babaee17
H. Babaee, M. Farazmand, G. Haller, and T. P. Sapsis.
Reduced-order description of transient instabilities and computation
  of finite-time Lyapunov exponents.
Chaos, 27(6):063103, 2017.

Blonigan_2017
P. J. Blonigan.
Adjoint sensitivity analysis of chaotic dynamical systems with
  non-intrusive least squares shadowing.
Journal of Computational Physics, 348:803–826, 2017.

Carlson_Lumley_1996
H. A. Carlson and J. L. Lumley.
Active control in the turbulent wall layer of a minimal flow unit.
Journal of Fluid Mechanics, 329:341–371, 1996.

Ceze_Diosady_Murman_2016
M. A. Ceze, L. T. Diosady, and S. M. Murman.
Development of a high-order space-time matrix-free adjoint solver.
AIAA 2016-0833, January 2016.

chen12
K. K. Chen, J. H. Tu, and C. W. Rowley.
Variants of dynamic mode decomposition: boundary condition,
  Koopman, and Fourier analyses.
Journal of Nonlinear Science, 22(6):887–915, 2012.

Chung_2015
D. Chung, L. Chan, M. MacDonald, N. Hutchins, and A. Ooi.
A fast direct numerical simulation method for characterising
  hydraulic roughness.
Journal of Fluid Mechanics, 773:418–431, 2015.

Carton_Diosady_2018
C. Carton de Wiart, L. T. Diosady, A. Garai, N. K. Burgess, P. J. Blonigan,
  D. Ekelschot, and S. M. Murman.
Design of a modular monolithic implicit solver for multi-physics
  applications.
In AIAA Aerospace Sciences Meeting, AIAA SciTech Forum,
  number 2018-1400. AIAA, 2018.

Carton_Murman_2017
C. Carton de Wiart and S. M. Murman.
Assessment of wall-modeled les strategies within a
  discontinuous-galerkin spectral-element framework corentin carton de wiart,
  nasa ames research center.
In 55th AIAA Aerospace Sciences Meeting and Exhibit, 2017.

dematteis2017
G. Dematteis, T. Grafke, and E. Vanden-Eijnden.
Rogue Waves and Large Deviations in Deep Sea.
Proceedings of the National Academy of Sciences,
  (10.1073/pnas.1710670115), 2018.

Diosady_Murman_2013
L. T. Diosady and S. M. Murman.
Design of a Variational Multiscale Method for Turbulent Compressible
  Flows.
In 21th AIAA Computaional Fluid Dynamics Conference, 2013.

Diosady_Murman_2014
L. T. Diosady and S. M. Murman.
Dns of flows over periodic hills using a discontinuous galerkin
  spectral element method.
AIAA paper, 2014.
accepted for publication.

Diosady_Murman_2015
L. T. Diosady and S. M. Murman.
Higher-order methods for compressible turbulent flows using entropy
  variables.
In 53rd AIAA Aerospace Sciences Meeting and Exhibit, 2015.

farano2017
M. Farano, S. Cherubini, J.-C. Robinet, and P. De Palma.
Optimal bursts in turbulent channel flow.
J. Fluid Mech., 817:35–60, 2017.

faraz_adjoint
M. Farazmand.
An adjoint-based approach for finding invariant solutions of
  Navier-Stokes equations.
J. Fluid Mech., 795:278–312, 2016.

Farazmand2016
M. Farazmand and T. Sapsis.
Dynamical indicators for the prediction of bursting phenomena in
  high-dimensional systems.
Physical Review E, 032212:1–31, 2016.

Farazmand2017a
M. Farazmand and T. P. Sapsis.
A variational approach to probing extreme events in turbulent
  dynamical systems.
Science Advances, 3(9):e1701533, 2017.

Farazmand2018a
M. Farazmand and T. P. Sapsis.
Extreme Events: Mechanisms and Prediction.
arXiv preprint arXiv:1803.06277, 2018.

Flores_Jimenez_2010
O. Florez and J. Jimenez.
Hierarchy of minimal flow units in the logarithmic layer.
Physics of Fluids, 22(7):071704, 2010.

Garai_Diosady_Murman_Madavan_2017
A. Garai, L. T. Diosady, S. M. Murman, and N. K. Madavan.
High-pressure turbine cascade using a spectral-element
  discontinuous-Galerkin method.
In Proceedings of ASME Turbo Expo 2017, no. GT2017-64697,
  Charlotte, USA, June 2017.

Holmes_et_al96
P. Holmes, J. Lumley, and G. Berkooz.
Turbulence, Coherent Structures, Dynamical Systems and
  Symmetry.
Cambridge University Press, 1996.

Jimenez_2012
J. Jimenez.
Cascades in wall-bounded turbulence.
Annual Review of Fluid Mechanics, 44:27–45, 2012.

jimenez2018
J. Jiménez.
Coherent structures in wall-bounded turbulence.
J. Fluid Mech., 842:P1, 2018.

Jimenez_Moin_1991
J. Jimenez and P. Moin.
The minimal flow unit in near-wall turbulence.
Journal of Fluid Mechanics, 225:213–240, 1991.

Jimenez_Pinelli_1999
J. Jimenez and A. Pinelli.
The autonomous cycle of near-wall turbulence.
Journal of Fluid Mechanics, 389:335–359, 1999.

Jimenez_Simens_2001
J. Jimenez and M. P. Simens.
Low-dimensional dynamics of a turbulent wall flow.
Journal of Fluid Mechanics, 435:81–91, 2001.

Scipy
E. Jones, T. Oliphant, P. Peterson, et al.
SciPy: Open source scientific tools for Python, 2001–.
[Online; accessed 7/11/18].

Kim_Moin_Moser_1986
J. Kim, P. Moin, and R. Moser.
Turbulence statistics in fully developed channel flow at low
  Reynolds number.
Journal of Fluid Mechanics, 177:133–166, 1986.

Kutz2015
J. N. Kutz, Xi. Fu, and S. L. Brunton.
Multi-Resolution Dynamic Mode Decomposition.
15(2):713–735, 2015.

MacDonald_2017
M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi, and R. Garcia-Mayoral.
The minimal-span channel for rough-wall turbulent flows.
Journal of Fluid Mechanics, 816:5–42, 2017.

mezic2013
I. Mezić.
Analysis of fluid flows via spectral properties of the Koopman
  operator.
Annual Review of Fluid Mechanics, 45:357–378, 2013.

mohamad2016b
M. A. Mohamad, W. Cousins, and T. P. Sapsis.
A probabilistic decomposition-synthesis method for the
  quantification of rare events due to internal instabilities.
Journal of Computational Physics, 322:288–308, 2016.

mohamad2015
M. A. Mohamad and T. P. Sapsis.
Probabilistic description of extreme events in intermittently
  unstable dynamical systems excited by correlated stochastic processes.
SIAM/ASA Journal on Uncertainty Quantification, 3(1):709–736,
  2015.

pringle2010
C. C. T. Pringle and R. R. Kerswell.
Using nonlinear transient growth to construct the minimal seed for
  shear flow turbulence.
Phys. Rev. Lett., 105(15):154502, 2010.

rowley09
C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson.
Spectral analysis of nonlinear flows.
J. Fluid Mech., 641:115–127, 2009.

sapsis2018
T. P. Sapsis.
New perspectives for the prediction and statistical quantification
  of extreme events in high-dimensional dynamical systems.
Phil. Trans. R. Soc. Lond. A, 376:20170133, 2018.

sapsis_majda_mqg
T. P. Sapsis and A. J. Majda.
A statistically accurate modified quasilinear Gaussian closure for
  uncertainty quantification in turbulent dynamical systems.
Physica D, 252:34–45, 2013.

sapsis_majda_tur
T. P. Sapsis and A. J. Majda.
Statistically accurate low-order models for uncertainty
  quantification in turbulent dynamical systems.
Proceedings of the National Academy of Sciences of the United
  States of America, 110(34), 2013.

schmid10
P. J. Schmid.
Dynamic mode decomposition of numerical and experimental data.
J. Fluid Mech., 656:5–28, 2010.

williams2015
M. O. Williams, I. G. Kevrekidis, and C. W. Rowley.
A data–driven approximation of the Koopman operator: Extending
  dynamic mode decomposition.
Journal of Nonlinear Science, 25(6):1307–1346, 2015.




