class: center, middle .title[Capturing Diffusive and Ballistic Transport from Lagrangian Observations] .author[Jonathan Lilly
1
] .coauthor[Sofia Olhede
2
, Adam Sykulski
1,2
, Jeffrey Early
1
] .institution[
1
NorthWest Research Associates,
2
University College London] .date[July 27, 2016] .center[[{www.jmlilly.net}](http://www.jmlilly.net)] .footnote[Created with [{Remark.js}](http://remarkjs.com/) using [{Markdown}](https://daringfireball.net/projects/markdown/) + [{MathJax}](https://www.mathjax.org/)] --- class: left ## Motivation Transport in turbulence may be broadly categorized into two types: 1. Diffusive transport, in which we can represent particle dispersion as a diffusive process; and 1. Non-diffusive transport where we cannot. One notable type of non-diffusive transport is the long-distance transport of material by the propagation of long-lived eddies. This kind of transport has the sense of material being thrown from one place to another, and has been called “ballistic” transport. --- class: center ## The Global Surface Drifter Dataset
Note: non-uniform sampling distribution, temporal variation, superposition of spatial scales --- class: center ## Mean Surface Current Speed
Formed by binning in latitude and longitude, then averaging. Easy to compute maps of low-order statistics: mean, variance, etc. Does not capture full richness of dataset. What else can be done? -- Strategy: create mathematical/statistical methods to isolate particular processes: (i) turbulent background (ii) eddies. --- class: center, middle #Part I: A Stochastic Model for Turbulent Dispersion --- class: center ## Particle Trajectories in 2D Turbulence
Shading is speed. Only non-eddy trapped trajectories are shown. --- class: left ## Modeling Dispersion in 2D Turbulence
We have identified a suitable, and very simple, stochastic model for representing Lagrangian trajectories in geostrophic turbulence. .cite[Lilly, J. M., A. M. Sykulski, J. J. Early, and S. C. Olhede (2016).
Fractional Brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion.
Submitted. Available at
{arXiv}.
] --- class: left ## A Hierarchy of Stochastic Models With advance apologies for notation: the following SDEs should be written as the corresponding stochastic integral equations. Note that `\(z(t)=u(t)+iv(t)\)` represents a complex-valued
velocity
. Brownian motion: `\[\frac{dz}{d t} = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{\omega^2}\]` Fractional Brownian motion (fBm): `\[\frac{d^\alpha z}{d t^\alpha} = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{\omega^{2\alpha}}\]` Damped Fractional Brownian motion a.k.a. the Matérn process: `\[\left[\frac{d}{d t} + \lambda \right]^\alpha z = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{(\omega^2+\lambda^2)^\alpha}\]` In Lilly et al. (2016), we showed that the relatively little-known Matérn process is in fact equivalent to adding a damping to fBm. --- class: left ## Another Hierarchy of Stochastic Models Damped Brownian motion
a.k.a.
Ornstein-Uhlenbeck (OU) process: `\[\frac{dz}{d t} + \lambda z = A\,\frac{d w}{d t} \label{maternsde}\quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{\omega^2+\lambda^2}\]` Damped Brownian motion + spin
a.k.a.
complex OU (cOU) process: `\[\frac{dz}{d t} - i\Omega z+ \lambda z = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{(\omega-\Omega)^2+\lambda^2}\]` Damped fBm + spin
a.k.a.
Matérn + spin
a.k.a.
fractional cOU: `\[\left[\frac{d}{d t} - i\Omega+ \lambda \right]^\alpha z = A\,\frac{d w}{d t} \quad\longrightarrow\quad S_{zz}(\omega) = \frac{A^2}{[(\omega-\Omega)^2+\lambda^2]^\alpha}\]` Matérn + spin generalizes cOU to fractional orders, thus unifying OU+cOU+fBm into a single larger family. --- class: left ##Diffusivity is the Spectrum's Value at Zero The velocity autocovariance and spectrum are defined as `\[R_{zz} (\tau)\equiv\left\langle z(t+\tau)\,z^*(t)\right\rangle = \frac{1}{2\pi}\int_{-\infty}^\infty e^{i \omega \tau} S_{zz}(\omega) \, d\omega\]` The corresponding diffusivity is given by `\[\kappa \equiv \lim_{t\longrightarrow\infty} \frac{1}{4} \frac{d}{d t} \left\langle x^2(t) + y^2(t)\right\rangle =\frac{1}{4} \int_{-\infty}^\infty R_{zz}(\tau)\, d \tau =\frac{1}{4}S_{zz}(0)\]` with the last equality following from the Fourier transform `\(\int_{-\infty}^\infty e^{-i \omega \tau} R_{zz}(\tau) \, d\tau= S_{zz}(\omega)\)` evaluated at zero frequency. Matérn: `\( \,\,\,\,S_{zz}(\omega) = \frac{A^2}{(\omega^2+\lambda^2)^\alpha} \quad\longrightarrow\quad \kappa = \frac{1}{4}\frac{A^2}{\lambda^{2\alpha}}\)` fBm: `\( \quad\quad S_{zz}(\omega)= \frac{A^2}{\omega^{2\alpha}} \quad\quad\,\,\longrightarrow\quad \kappa = \infty\)` Thus the Matérn, unlike fBm, can model diffusive processes. --- class: center ##Comparison with Other Models
--- class: center ##Application to the Global Drifter Dataset
Using a parametric fit to the Matérn, only on the anti-inertial side, the nondimensionalized damping parameter. Shorter ‘memory’ in more energetic regions. --- class: center ##Application to the Global Drifter Dataset
Using a parametric fit to the Matérn, only on the anti-inertial side, the first global map of the slope parameter `\(\alpha\)`, with `\(|\omega|^{-2\alpha}\)`. Slopes vary from `\(|\omega|^{-1}\)` to `\(|\omega|^{-3}\)`, with `\(|\omega|^{-2}\)` over major currents. Does not fit the conventional wisdom of only `\(|\omega|^{-2n}\)` slopes. --- class: left ### Modeling Trajectories in Ocean Turbulence -- We have developed a method for stochastic modeling of trajectories in ocean turbulence, and inferring parameters from large datasets. -- * Create an appropriate stochastic model for particle trajectories. .cite[Sykulski, Olhede, Lilly, and Danioux (2015)] -- * A key ingredient is a *damped* version of fractional Brownian motion. .cite[Lilly, Sykulski, Early, and Olhede (2016), submitted] -- * Parameter estimation is best done in the frequency domain. .cite[Whittle (1953)] -- * Parameter estimation must be adjusted to handle *complex-valued* time series. .cite[Sykulski, Olhede, Lilly, and Early (2016a), submitted] -- * Parameter estimation must be corrected for bias due to small sample size effects. .cite[Sykulski, Olhede, Lilly, and Early (2016b), in preparation] --- class: left ##Summary, Part I The Matérn process is the simplest random process that can simultaneously reproduce (i) the velocity variance (ii) the diffusivity and (iii) the spectral slope or degree of small-scale roughness. This is a promising stochastic model for Lagrangian trajectories, and diffusive processes in general. It is useful for (i) parameteric spectral modeling, (ii) as model for synthetic trajectories, (iii) as a null hypothesis in eddy studies. --- class: center, middle # Part II: Extracting Coherent Eddy Motions --- class: left ## Goal Create a method capable of automatically extracting velocity signals of particles trapped in coherent eddies, a.k.a. “loopers”. Use these to estimate physical properties such as radius, vorticity, Rossby number, etc., and eventually populations and fluxes. -- ## Why this is difficult The oscillations of a particle trapped in an eddy are highly nonstationary due to eddy evolution, non-Lagrangian behavior, etc. In particular, frequency is unknown, and may occupy a broad band. Thus, neither Fourier nor baseband demodulation methods are suitable. -- ## Method Our approach is to use nonstationary signal analysis techniques, using the wavelet transform as a “lens”, with the analytic signal theory providing a foundation. --- class: center ### A Numerical Simulation of an Unstable Current
A highly idealized version of the Gulf Stream. Note formation of vortices or “coherent eddies”. --- class: center ### Identifying Vortices from Particle Trajectories
Vortices (loops) are identified using only particle trajectories (dots). --- class: left ### How Does this Work? First, we conceptualize the orbit of a particle trapped in an eddy as a
time-varying ellipse.
`\[z_o(t)=x_o(t)+\mathrm{i}y_o(t)=\mathrm{e}^{\mathrm{i}\theta(t)}\left[a(t)\cos\phi(t) + \mathrm{i} b(t)\sin \phi(t)\right]\]`
We construct a type of “filter” that can isolate such signals. This lets us approximately split a time series `\(z(t)=z_o(t)+z_\epsilon(t)\)` into an oscillatory `\(z_o(t)\)` and stochastic `\(z_\epsilon(t)\)` portion. --- class: center ### Extraction of Modulated Elliptical Signals
Vortices can be identified, studied, and distinguished from waves. .cite[Lilly, Scott, and Olhede (2011), *GRL*] --- class: center ## Radius-Velocity Distribution of Signals
Coherent eddies are primarily high-frequency and are circularly polarized. Jet meander is lower-frequency and linearly polarized. --- class: left ##Recovering Modulated Oscillations The wavelet transform is a tool for
recovering
or
estimating
the properties of modulated oscillation immersed in background noise. The idea is to project the time series onto an oscillatory test signal, or
wavelet
, to find a “best fit” frequencies at each moment. The fits are then chained together into a continuous curve called a
ridge.
For a vector-valued signal `\(\mathbf{x}_o(t)=\begin{bmatrix}x_1(t) & x_2(t) & \cdots & x_N(t)\end{bmatrix}^T\)` in noise `\(\mathbf{x}_\epsilon(t)\)`, `\(\mathbf{x}(t)=\mathbf{x}_o(t)+\mathbf{x}_\epsilon(t)\)`, define the wavelet transform `\[\mathbf{w}_\mathbf{x}(t ,s) \equiv \int_{-\infty}^{\infty} \frac{1}{s} \psi^*\left(\frac{\tau-t}{s}\right)\,\mathbf{x}(\tau)\,\mathrm{d} \tau\]` (also a vector) and then find the
wavelet ridges
`\(s(t)\)` from `\[\frac{\partial}{\partial s}\, \left\|\mathbf{w}_\mathbf{x}(t ,s)\right\| = 0,\quad\quad \frac{\partial^2}{\partial s^ 2}\, \left\|\mathbf{w}_\mathbf{x}(t ,s)\right\| < 0.\]` The oscillation is estimated simply by `\(\widehat{\mathbf{x}_o}(t)\equiv\Re\left\{\mathbf{w}_\mathbf{x}(t,s(t))\right\}.\)` .cite[See Delprat et al. (1992), Lilly and Olhede (2009, 2010, 2012a). ] --- class: center ##Example of Multivariate Ridge Analysis
The
ridge
is the curve made by tracing out the maximum modulus of the wavelet transform. --- class: center ##Example of Multivariate Ridge Analysis
The transform value along this curve is an estimate of the oscillatory signal component, visualized as ellipses. Ridge analysis separates modulated elliptical signal from low-frequency meandering and higher-frequency variability. Use `\(\mathtt{wavetrans}\)` and `\(\mathtt{ridgewalk}\)` from jLab at [{www.jmlilly.net}](http://www.jmlilly.net). --- class: center ###Deep Connection to Vortex Dynamics
Blue = time-varying integrated vorticity, red = one-point estimate. From a simulation by Jeffrey Early. --- class: left ###Vorticity from a single Lagrangian trajectory An analytic expression for vortex circulation estimated from a single Lagrangian trajectory: `\[\Gamma(t)= \pi \Im \left\{ \mathbf{x}_+^H(t) \frac{\mathrm{d}}{\mathrm{d} t} \mathbf{x}_+ (t) \right\}\]` where `\(\mathbf{x}_+(t)\)` is a complex-valued version of `\(\mathbf{x}(t)\)`, the
analytic signal
`\[\mathbf{x}_+(t)\equiv\mathbf{x}(t)+\mathrm{i}\mathcal{H}\mathbf{x}(t) \equiv \mathbf{x}(t)+\mathrm{i}\frac{1}{\pi}\,\int_{-\infty}^\infty\frac{\mathbf{x}(u)}{t-u}\,\mathrm{d} u\]` with `\(\mathcal{H}\)` being the Hilbert transform. For a particle on a slowly varying elliptical material curve, this expression is
exactly correct
. Uniqueness: If you wish to infer the eddy properties using a
linear filter
, and if you wish your method to obtain the correct results for a strictly periodic eddy flow (harmonic correspondence), there is
no alternative
to using the analytic signal. Rewording Vakman (1996). --- class: left ##Sufficient Condition for Recoverability Q. Under what conditions are the time-varying properties of an ellipse exactly recoverable from the trajectory of one particle? -- A. When the ellipse geometry evolves more slowly than circulation of particles around the periphery. --
Ellipse parameters for top half are exactly recoverable. --- class: center ## Application to the Global Drifter Dataset
Apply eddy detection algorithm to 20K time series. Use a range of frequencies compared to the Coriolis frequency `\(f\)`. --- class: center ### Eddy Detection in the Global Drifter Dataset
Red = cyclonic rotation, blue = anticyclonic. --- class: center ### Statistically Significant Eddies
After comparison with a null hypothesis of red noise. --- class: left ### Method for Extracting Eddy Currrents -- We have developed a general method for extracting and analyzing *time-varying* and quasi-periodic signals, such as eddy currents. * The notion of the *analytic signal* lets an *instantaneous* amplitude and frequency be associated with a given time series. .cite[Gabor (1946), Vakman and Vainshtein (1977), Picinbono (1997)] -- * Similarly, a bivariate (x,y) signal defines the geometrical properties of a time-varying ellipse. .cite[Lilly and Olhede (2010a)] -- * This may be extended to 3D (e.g. seismographs) or N-D signals. .cite[Lilly (2011), Lilly and Olhede (2012a)] -- * Modulated oscillations can be extracted using *wavelet ridge analysis*, a local best fit onto an oscillatory test function. .cite[Delprat et al. (1992), Mallat (1999), Lilly and Olhede (2010b)] -- * Modulated multivariate oscillations can be treated similarly .cite[Lilly and Olhede (2012a)] -- * Errors are proportional to modulation strength, and are minimized by a suitable choice of wavelet. .cite[Lilly and Olhede (2010b), Lilly and Olhede (2012a)] -- * Best choice of wavelet is found by considering a superfamily encompassing all other continuous analytic wavelets. .cite[Lilly and Olhede (2009, 2012b)] --- class: left ## Summary, Part II We have shown that coherent eddy velocity signals can be extracted from the turbulent background, and that polarization and frequency can distinguish eddies from waves. By isolating eddies, we will be able to examine their statistical distributions, propagation characteristics, and contribution to non-diffusive or ballistic transport. -- Taken together, these two methods will let us isolate and separately quantify two very different contributors to turbulent transport. --- class: center, middle # Thank you! .center[Visit [{www.jmlilly.net}](http://www.jmlilly.net) for papers, this talk, and a Matlab toolbox of all numerical code.] .center[.footnote[P.S. Like the way this presentation looks? Check out [{Liminal}](https://github.com/jonathanlilly/liminal).]]