class: center, middle .title[Progress in Stochastic Modeling of Lagrangian Trajectories]           .author[Jonathan Lilly¹] <!--.center[[www.jmlilly.net](http://www.jmlilly.net)] (CICESE)-->   .coauthor[Adam Sykulski², Sofia Olhede³, Jeffrey Early⁴]   .institution[¹Theiss Research, ²Lancaster University, ⁴École Polytechnique Fédérale de Lausanne, ⁴NorthWest Research Associates<br/> ]     .date[June 21, 2019]   .footnote[Created with [{Remark.js}](http://remarkjs.com/) using [{Markdown}](https://daringfireball.net/projects/markdown/) + [{MathJax}](https://www.mathjax.org/)] --- class: left ## An Aside on Kinematic Properties In a 2D flow in which velocity depends linearly on the spatial coordinates, `$\mathbf{u}=\mathbf{U} \mathbf{x}$`, there is a geometric intepreation of the vorticity `$\zeta$`, divergence `$\delta$`, and normal and shear strains `$\nu$` and `$\sigma$`. Draw any ellipse on the flow. The following is true. “For any ellipse advected by an arbitrary linear two-dimensional flow, the rates of change of the ellipse parameters are uniquely determined by the four parameters of the velocity gradient matrix, and vice versa. This result, termed *ellipse/flow equivalence*, provides a stronger version of the well-known result that a linear velocity field maps an ellipse into another ellipse. ” .cite[(Lilly, 2018)] More generally, if you Taylor-expand any 2D flow about a point where the linear dependence on velocity does not vanish, then this should approximately hold for some suitably small ellipse. --- class: center ## Ellipse/Flow Equivalence `\begin{equation}\label{basisdef} \mathbf{I} \equiv \begin{bmatrix} 1 & 0 \\ 0 &1\end{bmatrix},\quad \mathbf{J} \equiv \begin{bmatrix} 0 & -1 \\ 1 & 0\end{bmatrix},\quad \mathbf{K} \equiv \begin{bmatrix} 1 & 0 \\ 0 &-1\end{bmatrix},\quad \mathbf{L} \equiv \begin{bmatrix} 0 & 1 \\ 1 & 0\end{bmatrix} \end{equation}` `\begin{equation}\label{Udecomp2} \mathbf{U}(t) = \frac{1}{2}\left(\delta \mathbf{I} + \zeta \mathbf{J}+\nu\mathbf{K} +\sigma\mathbf{L}\right) \end{equation}` <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/ijkl-vectors.jpg"> `\begin{multline} \mathbf{U}=\mathbf{R}(\theta)\left\{ \frac{d \ln \rho}{d t}\mathbf{I} +\left(\frac{d\theta}{d t}+\frac{\eta^2+1}{2\eta}\frac{d\phi}{d t}\right)\mathbf{J}\right.\\\left. + \frac{1}{2}\frac{d\ln \eta}{d t}\mathbf{K}-\frac{\eta^2-1}{2\eta}\frac{d\phi}{d t}\mathbf{L}\right\}\mathbf{R}^T(\theta) \end{multline}` `$x+iy=\mathrm{e}^{\mathrm{i}\theta(t)}\left[a(t)\cos \phi(t)+\mathrm{i}b(t)\sin \phi(t)\right]$`, `$\rho\equiv \sqrt{ab}$`, `$\eta=a/b$` --- class: left ## Overview We have found a simple three-parameter stochastic model that adequately describes the essential second-order features of Lagrangian velocities. One is the absolute diffusivity, one is a time scale related to Lagrangian integral time scale, and the third is the spectral slope. These three parameters provide a convenient way to summarize the bulk properties of a Lagrangian time series. It's essential to have all three. We don't have any theory for why the Lagrangian velocities should match this form. More generally, there has been almost no work done on Lagrangian spectral slopes—which makes this a promising area for future work. --- class: center ##Ocean Turbulence <video preload="auto" width="60%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/turbulencemovie.mp4" type="video/mp4" /></video> An idealized model of oceanic turbulence by Jeffrey Early. Shading is speed with particle trajectories shown in color. --- class: center ##Modeling Lagrangian Trajectories <video style="margin-top:-0.5em;" preload="auto" width="100%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/dispersionmovie.mp4" type="video/mp4" /></video> <!-- One of these is 2D turbulence, one is a 3-parameter stochastic model.--> A simple stochastic model for trajectories in oceanic turbulence is the *Matérn process* or <i>damped</i> Fractional Brownian motion. This three-parameter model lets us match (i) the kinetic energy, (ii) the dispersion, and (iii) the degree of small-scale roughness. --- class: left ##The Matérn Process A relatively little-known random process called the Matérn process .cite[(Matérn, 1960)] has a spectrum given by, with `$z(t)=u(t)+iv(t)$`, `\[S_{zz}(\omega) = \frac{A^2}{(\omega^2+\lambda^2)^\alpha}\]` where `$A$` sets the energy level, `$\alpha$` sets the high-frequency slope, and `$\lambda$` determines a low-frequency transition to a constant spectral value. Recently, we investigated this process in detail, and showed that it is an excellent model for trajectories in ocean turbulence. We speculate it will prove suitable for many other processes as well. The Matérn process generalizes of *fractional Brownian motion* .cite[(Mandelbrot and Van Ness, 1968)] to include a damping term, `$\lambda$`, that is the essential ingredient for modeling dispersion. .cite[Lilly, Sykulski, Early, and Olhede (2017). Fractional Brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion. <a href="./papers/lilly17-npg.pdf"><nobr>[3.2 Mb pdf]</nobr> </a>] --- class: left ## An Hierarchy of Stochastic Models We have the following correspondences between stochastic differential equations (SDEs) and velocity spectra `$S_{zz}(\omega)$`. 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}\]` (These should be written as the corresponding *integral* equations.) --- class: left ## Another Hierarchy of Stochastic Models Damped Brownian motion <i>a.k.a.</i> Ornstein-Uhlenbeck (OU) process: (Same as a damped-slab mixed layer model with white noise forcing.) `\[\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 <i>a.k.a.</i> complex OU (cOU) process: (Same spin parameter used by Veneziani et al. to model eddies.) `\[\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 <i>a.k.a.</i> Matérn + spin <i>a.k.a.</i> 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}\]` <!--Note spin is a deterministic effect!--> <!--`\[z(t) = e^{i\Omega t} \tilde z(t)\longrightarrow R_{z z}= e^{i\Omega \tau} R_{\tilde z\tilde z}(\tau) \longrightarrow S_{zz}(\omega)=S_{\tilde z\tilde z}(\omega-\Omega) \]`--> Matérn + spin generalizes cOU to fractional orders, thus unifying OU+cOU+fBm into a single larger family. --- class: left ##Meanings of the Parameters The Matérn process is controlled by three parameters, `\[S_{zz}(\omega) = \frac{A^2}{(\omega^2+\lambda^2)^\alpha}.\]` Fractional Brownian Motion is contolled by two, `\[S_{zz}(\omega) = \frac{A^2}{\omega^{2\alpha}}.\]` What are the meanings of these parameters? --- class: left ##Slope and Roughness <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/fbmroughness.png"> The `$\alpha$` parameter sets the slope in fractional Brownian motion, or at high frequencies in the Matérn process, as `$S(\omega)\sim \omega^{-2\alpha}$`. It also sets the degree of small-scale roughness in the trajectories. In the above, we go from slopes of `$\omega^{-3}$` at left to `$\omega^{-1}$` at right. Steeper spectral slope = more smooth. Spectral slope can be formally linked to a mathematical measure of roughness called the *fractal dimension*, though in my opinion, fractal dimension seems to have little practical value. --- class: center ##Slope and Self-Similarity <img style="width:90%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/fbmselfsimilarity.png"> The `$\alpha$` parameter also sets *aspect ratio* of statistical self-similarity. --- 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 absolute 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 `$\int_{-\infty}^\infty e^{-i \omega \tau} R_{zz}(\tau) \, d\tau= S_{zz}(\omega)$`. 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$` Infinity means the diffusivity increases without bound. Thus the Matérn, unlike fBm, can model diffusive processes. --- class: left ##Another Take on Diffusivity Most simply, the absolute diffusivity can be understood as a fundamental second-order statistical property of the Lagrangian velocity. The diffusivity is the twin of the variance. `\[\sigma^2 =\frac{1}{2\pi} \int_{-\infty}^\infty S_{zz}(\tau)\, d \tau =R_{zz}(0)\]` `\[\kappa =\frac{1}{4} \int_{-\infty}^\infty R_{zz}(\tau)\, d \tau =\frac{1}{4}S_{zz}(0)\]` Thus diffusivity is a basic second-order property of any time series. The zero frequency value of the spectrum of a process describes how the temporal integral of the process spreads out over time. You don't have to think that diffusion is a good physical model in order to usefully quantify a process using `$\kappa$` as one of the parameters. --- class: left ##The Damping Time Scale The `$\lambda$` parameter sets the frequency, or timescale `$2\pi/\lambda$`, at which the Matérn process transitions from behaving like fractional Brownian motion to behaving like white noise: `\[S_{zz}(\omega) = \frac{A^2}{(\omega^2+\lambda^2)^\alpha}.\]` `\[S_{zz}(\omega) \approx \frac{A^2}{\omega^{2\alpha}},\quad\quad\quad \omega\gg \lambda\]` `\[S_{zz}(\omega) \approx \frac{A^2}{\lambda^{2\alpha}},\quad\quad\quad \omega\ll \lambda\]` If you zoom in to the velocity process, you see fBm. If you zoom out, you see white noise (and thus finite diffusivity). Note that `$2\pi/\lambda$` is not the same as the Lagrangian integral timescale. --- class: center ##Comparison of Models <img style="width:65%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/randomwalk_dispersion.png"> Top left: numerical simulations. Top right: Matérn process. Bottom left: white noise. Bottom right: red noise (Brownian). --- class: left ##Argument for Generality We conjecture that the Matérn may be of general usefulness: “Systems are often characterized by a pressure to grow—represented by a forcing—together with some drag or resistance on that growth, represented by a damping. After a sufficiently long time, the forcing and the damping equilibrate and one reaches a bounded state. This leads to the speculation that many time series that are well described as fBm over relatively short timescales may be better matched by the Matérn process over longer timescales. More generally, the Matérn process adds a third parameter (damping) to the two parameters (amplitude together with spectral slope or the Hurst parameter) of fBm, thus permitting a wider range of spectral forms to be accommodated. It is therefore reasonable to think that the Matérn process could be of broad interest in many areas in which fBm has already proven itself useful.” .cite[Lilly, Sykulski, Early, and Olhede (2017)] --- class: left ##A Conjecture on the Spectral Slope Is there a general relationship between wavenumber slope `$k^{-\beta}$` and frequency slope `$\omega^{-2\alpha}$` ? Here is a naïve try: <img style="width:90%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/slopeconjecture.png"> --- class: center, ## The Global Surface Drifter Dataset <video preload="auto" width="100%" height="auto" data-setup="{}" autoplay loop controls><source src="../videos/driftermovie.mp4" type="video/mp4" /></video> Now we will apply the Matérn process to the global drifter dataset. --- class: center, ## Spectra of Drifter Velocities <img style="width:45%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/highreslatitudespectra.png"> We see a low-frequency maximium and tidal maxima—all vertical lines—and a maximum along a curving line, the inertial frequency. <!--Inertial frequency becomes the vertical line `$\omega/f_o=1$`. Tides become curved. Note roughly universal appearance in new coordinates.--> --- class: center, ##A Conceptual Model for the Spectrum <img style="width:65%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/schematictransparency2.png"> \begin{equation} S(\omega) = \overset{\mathrm{background}}{\overbrace{\frac{A_b^2}{\left[\omega^2+\lambda_b^2\right]^\alpha}}} +\overset{\mathrm{inertial\,\, oscillations}}{\overbrace{\frac{A_o^2}{\left(\omega-f_o\right)^2+\lambda_o^2}}}+ \overset{\mathrm{semidiurnal\,\,tide}}{\overbrace{\frac{A_s^2}{\left(\omega-f_s\right)^2+\lambda_s^2}}} \end{equation} This is just several versions of the Matérn spectrum added together. --- class: left ## Details of the Stochastic Modeling We have developed a method for stochastic modeling of trajectories in ocean turbulence, and inferring parameters from large datasets. This is an example of *parametric spectral analysis*: rather than trying to estimate the spectrum directly, one begins with a *conceptual model*, controlled by some parameters, and then finds the choice of parameters which provide the best match to the data. The full method consists of several ingredients: * Create an appropriate stochastic model for particle trajectories. .cite[Sykulski, Olhede, Lilly, and Danioux (2015)] * Employ the Matérn process as a building block. .cite[Lilly, Sykulski, Early, and Olhede (2017)] * 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 (2017)] * Parameter estimation must be corrected for subtle bias due to small sample size effects. .cite[Sykulski, Olhede, Guillaumin, Lilly, and Early (2019)] --- class: center ## Spectra in Latitude Bands <img style="width:75%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/drifterspectralinestransparency.png"> Averages in `$\pm 5 ^\circ $` bins. Frequency is nondimensionalized as `$\omega/f_o$`. Inertial peak broadens equatorward. --- class: center ## Comparison with Spectral Fit <img style="width:75%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/drifterspectralinestransparencywithfullfit.png"> In general, spectral model compares very well with estimated spectra. Inertial and semidiurnal peaks are well captured. --- class: center ## Fit to Background vs. Cyclonic Side <img style="width:75%;margin-top:-0.7em;margin-bottom:-0.8em" src="../figures/drifterspectralinestransparencycyclonicwithfit.png"> “Background” fit to anticyclonic (inertial side), matches well the spectrum on the cyclonic (anti-inertial) side, apart from tides. --- class: center ##Application to the Global Drifter Dataset <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/globalantiinertiallambda.png"> The nondimensionalized damping parameter from a parametric fit to the Matérn, only on the anti-inertial side. Shorter ‘memory’ in more energetic regions. We are not aware of any theory for this. --- class: center ##Application to the Global Drifter Dataset <img style="width:100%;margin-top:-0.2em;margin-bottom:-0.2em" src="../figures/globalantiinertialalpha.png"> 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, as in e.g. .cite[Berloff and McWilliams (2002b)] <!--Lagrangian integral timescale: `\[T \equiv \frac{1}{2}\frac{\int_{-\infty}^\infty R_{zz}(\tau) \,d \tau}{R_{zz}(0)} = \frac{1}{2}\frac{S_{zz}(0)}{\frac{1}{2\pi} \int_{-\infty}^\infty S_{zz}(\omega) \,d \omega} =\frac{2\kappa}{\sigma^2}\]`--> <!--Spectral slope / fractal dimension: `\[R_{zz}(\tau) \sim |\tau|^{2\alpha-1}, \quad\tau \approx 0 \quad \longleftrightarrow\quad S_{zz}(\omega) \sim |\omega|^{-2\alpha},\quad \omega\gg 0\]`--> <!--which holds for $1/2<\alpha<3/2$; see \citet{goff88-jgr},--> <!-- $\sigma^2 =c_\alpha A^2/\lambda^{2\alpha-1}$ and $\kappa = \frac{1}{4}A^2/(\Omega^2+\lambda^2)^\alpha$, where $c_\alpha\equiv\Gamma\left(\alpha-\frac{1}{2}\right)/\Gamma(\alpha)/(2\sqrt{\pi})$ with $\Gamma(\cdot)$ being the gamma function, see \citet{lilly16-itit}, while the Lagrangian integral timescale is just~$T=2\kappa/\sigma^2$. --> --- class: left ##Summary Understanding the form of Lagrangian velocity spectrum, in particular the spectral slope, is a promising avenue for future work. 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. It could be useful for summarizing velocity statistics through parameteric spectral modeling, as a model for synthetic trajectories, and as a null hypothesis in e.g. eddy studies. The Matérn process is conjectured to be of broad applicability. The paper describing the Matérn process is intended to serve as an accessible introduction to stochastic processes in general. .cite[Lilly, Sykulski, Early, and Olhede (2017). <a href="https://doi.org/10.5194/npg-24-481-2017">{Fractional Brownian motion, the Matérn process, and stochastic modeling of turbulent dispersion}</a>, <i>Nonlinear Processes in Geophysics</i>.] Code can be found in [{jLab}](http://www.jmlilly.net). --- class: center, middle # Thank you! .center[This talk is available at].center[[http://www.jmlilly.net/talks/lilly19b-lapcod.html](http://www.jmlilly.net/talks/lilly19b-lapcod.html)] .center[.footnote[P.S. Like the way this presentation looks? Check out [{Liminal}](https://github.com/jonathanlilly/liminal).]]