Something I like to do at large research conferences is drop into sessions with intriguing titles about topics unfamiliar to me. At the last SIAM Conference on applications of dynamical systems, I dropped into “Set Oriented and Transfer Operator Methods for Turbulent Flows”. I learned about an appealing visualization called finite time Lyapunov exponent fields that I’ll present here.

Turbulence “is a flow regime in fluid dynamics characterized by chaotic changes in pressure and flow velocity (wikipedia: Turbulence).” Datasets of turbulent flows discussed by presenters in the session included ocean currents and particle tracers in a controlled laboratory fluid vat. The presenters and other researchers in the field approach the broad question “How do turbulent flows evolve over time?” and commonly use the approach of tracking coherent structures in the flow. Coherent structures are distinguished trajectories that exert a major influence on the overall flow. Examples include the water vortices classically illustrated by Leonardo da Vinci.

To my understanding, traditional approaches classify coherent structures with respect to the reference frame of a basis for the flow domain, yielding Eulerian coherent structures, but more recent approaches classify them with respect to the reference frame of the flow itself, yielding Lagrangian coherent structures. Distinguished sets in finite time Lyapunov exponent fields are examples of the latter.

Proper orthogonal decomposition, which refers to computing the SVD of a data matrix whose columns are vectorized observables from a grid superimposed on the flow domain, is a workhorse for detecting Eulerian coherent structures. The strongest modes of the SVD correspond to coherent structures and can appear as, say, vortices in the flow.

From dynamical systems theory, invariant manifolds emanate from the equilibria in flows generated by time-independent vector fields. These manifolds, and other Lagrangian coherent structures, partition the state space. They offer a powerful tool for understanding the entire flow since the asymptotic fate of any initial condition can be classified by simply identifying the state space partition in which it begins.

Many real-world flows, however, are generated by time-dependent vector fields and often lack equilibria, so the theory of time-independent invariant manifolds does not apply. There is a computational approach to approximate analogous separatrices for flows generated by time-dependent vector fields. As I will attempt to make more precise below, the finite time Lyapunov exponent quantifies at each point in the flow domain the degree to which nearby points diverge when advected by the flow.

Let $F:U\times[-c,c]\rightarrow M$, defined by $F(x_0;t)=x(t)$, be the flow that takes an initial condition $x_0\in U$ at time $t_0=0$ to position $x(t)$ at time $t\in[-c,c]$, where $U$ is an open set of manifold $M$ and $c$ is a real number. A trajectory of $F$ starting at $x_0$ with an initial perturbation $\Delta x_0$ evolves in the flow by $F(x_0+\Delta x_0;t)=x(t)+\Delta x(t)$. Expanding the lefthand side about $x_0$ yields a linear approximation of the evolving perturbation $\Delta x(t)$. Assuming $||\Delta x_0||=1$, the squared-magnitude of the evolving perturbation

\begin{equation}\label{eq:rayleigh} ||\Delta x(t)||_2^2\approx\Delta x_0^T\nabla F(x_0;t)^T\nabla F(x_0;t)\Delta x_0 \end{equation}

is bounded by the magnitude of the maximum eigenvalue of $\nabla F(x_0;t)^T \nabla F(x_0;t)$ since the righthand side of \eqref{eq:rayleigh} is a Rayleigh quotient of $\nabla F^T\nabla F$ (a matrix called the right Cauchy-Green deformation tensor). The scalar field

\begin{equation}\label{eq:ftle} \sigma(x;t)={1 \over |t|}\ln\sqrt{\lambda_{\max}(\nabla F(x;t)^T\nabla F(x;t))}, \end{equation}

where $\lambda_{\max}(\cdot)$ stands for maximum eigenvalue, is called the finite time Lyapunov exponent field of $F$ parameterized by $t$.

Ridges or elevated regions in the finite time Lyapunov exponent field $\sigma(x;t)$ act as barriers that partition the flow into sections that tend not to mix. Think of a saddle equilibrium in a two-dimensional flow generated by a time-independent vector field. Two initial conditions on either side of a branch of the saddle’s stable manifold tend toward opposite branches of the saddle’s unstable manifold, resulting in large finite time Lyapunov exponent values near the stable manifold. Especially in the context of flows generated by time-dependent vector fields, these ridges are examples of Lagrangian coherent structures.

Computing \eqref{eq:ftle} requires being able to approximate the flow $F$ and its first derivatives. This can be done directly from data. Numerical integration can be used to evolve the flow while interpolation is used to evaluate it at any point in the domain. Finite differences can be used to approximate $\nabla F$. As an example, I compute $\sigma(x;t)$ for Earth wind currents in this dataset that collects satellite measurements on a $628\times 1440$ latitude-longitude grid of wind vectors at 6-hr time resolution over 25 years.

Here’s a sample of the data for a month-long window (8/14-9/15/05) during which hurricane Katrina hit New Orleans. The hurricane appears as a rogue gust in the Gulf of New Mexico about halfway through the video.

I used forward Euler to do the numerical integration and cubic b-splines (via FITPACK) to do interpolation on a sphere with the Earth’s radius. Here are the resulting ridges in the finite time Lyapunov exponent field parameterized by $t$ equal to 2 weeks for a sliding window over the same month-long period.

Computing $\sigma(x;t)$ involves advecting a grid of trajectories by the flow. This is highly parallelizable. To capitalize on this, I wrote the numerical integration and differentiation routines in Fortran that I wrapped with Python (via f2py) and launched multiple processes (via Python’s multiprocessing module) to carry out parallel computations.