class: center, middle .title[New Directions in Oceanographic Time Series Analysis] .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[March 22, 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: 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. Clearly does not capture full richness of dataset. What else can be done? --- class: left ## Properties of the Surface Drifter Dataset -- 1. Non-homogeneously distributed in space -- 2. Non-uniformly sampled in time -- 3. Irregular length or duration -- 4. Non-stationary statistics -- 5. Superposition of various physical processes -- ## Processes of interest -- 1. Low-frequency diffusive behavior -- 2. Intermediate-frequency coherent eddy structures -- 3. Higher-frequency inertial oscillations -- ##Strategy -- Construct particular analysis methods aimed at isolating various physical processes. -- Initially, work within simplified dynamics of numerical models. --- class: center, middle # Part I: Coherent Eddies as Modulated Oscillations --- 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 ### Analyzing Modulated Multivariate Oscillations -- We have developed a method for analyzing and interpreting *time-varying* properties of quasi-periodic or quasi-oscillatory signals. * 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 ## Conclusions, Part I 1. Coherent eddy velocity signals can be recovered 1. Polarization and frequency can distringuish eddies from waves 1. Random variability may sometimes appear as eddy signals 1. Excluding random variability recovers major eddy pathways To be continued... --- class: center, middle # Part II: Inertial Oscillations and Dispersion --- class: left ###Goal The **goal** is to analyze the detailed properties of inertial oscillations—including amplitude and damping time scale—from the global surface drifter dataset. -- This is **important** for understanding the input of mechanical energy to the ocean by the wind, its loss to the internal wave field, and eventual dissipation. -- Related work includes .cite[Chaigneau et. al (2008)], .cite[Park et. al (2009)], and .cite[Elipot et. al (2010)]. --- class: left ##Method This problem is **difficult** for two reasons: 1. Inertial oscillations have an unknown **bandwidth**. 2. **Interference** from low-frequency and tidal energy. -- So, you can't just bandpass the data to get globally meaningful estimates of amplitude, and especially, the damping timescale. -- The eddy-extraction method of .cite[Lilly, Scott, and Olhede (2011)], intended for the case of unknown **frequency**, doesn't work here. -- Actually, this turns out to be a really difficult problem. The basic approach is to represent the drifter trajectories as if they were **stochastic processes**, and then to estimate the parameters controlling those processes. -- Method papers: * .cite[Sykulski, Olhede, Lilly, and Danioux (2016), JRSSC. ] * .cite[Sykulski, Olhede, Lilly, and Early (2016), in prep. ] * .cite[Lilly, Sykulski, Early, and Olhede (2016), in prep. ] --- class: left ## Stochastic Model for Lagrangian Velocities We represent Lagrangian velocities as the sum of inertial oscillations, plus a background of geostrophic turbulence + tides: \begin{equation} z(t)=u(t)+\mathrm{i}v(t) = \overset{\mathrm{inertial}}{\overbrace{z_o(t)}}\, + \overset{\mathrm{background}}{\overbrace{z_b(t)}} +\color{red}{\overset{\mathrm{semidiurnal}}{\overbrace{z_s(t)}}} \end{equation} The inertial oscillation component can be approximated by the
damped-slab
model of .cite[Pollard and Millard (1970)], \begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} z_o(t) + \mathrm{i} f_o z_o(t) +\lambda z_o(t) = \tau(t) /\left(\rho H(t)\right) \end{equation} where `\(f_o=2\Omega \sin \phi\)` is the signed Corilios frequency, `\(\tau=\tau_x+\mathrm{i} \tau_y\)` is the wind stress, and `\(H\)` is the effective mixed layer depth. Our goal is to
infer
the inertial oscillation properties—the typical energy `\(\sigma^2\equiv\langle |z_o(t)|^2\rangle\)` and damping parameter `\(\lambda\)`—from the data. --- class: center ## Rotary Spectrum of the Conceptual Model
\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}}}+ \color{red}{\overset{\mathrm{semidiurnal\,\,tide}}{\overbrace{\frac{A_s^2}{\left(\omega-f_s\right)^2+\lambda_s^2}}}} \end{equation} --- class: center ## Density of Drifter Position Measurements
Total of 12,122 trajectories spanning 73,117,043 data points. Irregular temporal sampling. Median sample interval of 1.22 hours. --- class: center ## Density of 60-Inertial Period Segments
Total of 119,505 segments with 50% overlap. Each is 60 inertial periods in duration. Locations of segment-averaged latitude and longitude are shown. QC'd for gaps, data sparsity, etc. --- class: center ## Rotary Spectra of Drifter Velocities
Left-hand side is after Elipot and Lumpkin (2008). Right-hand side has frequency normalized as `\(\omega/f_o=1\)`, and inertial (anticyclonic) rotations flipped to positive frequencies. Tides become curved. --- class: center ## Rotary Spectrum of the Conceptual Model
\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} --- class: center ## Particle Trajectories in 2D Turbulence
Shading is the speed. Of 1024 particle trajectories, only 512 are shown, excluding those that exhibit trapping in eddies. --- class: center ## Modeling the Turbulent Background Flow
One of these is 2D turbulence, one is a stochastic model. We have identified a suitable, and very simple, stochastic model for representing Lagrangian trajectories in geostrophic turbulence. This is known as the
Matérn
process, which we show to be equivalent to **damped fractional Brownian motion**. .center[.cite[Lilly, Sykulski, Early, and Olhede (2016), in prep. ]] --- class: center ### Spectra in Latitude Bands
Averages in `\(\pm 5 ^\circ \)` bins. Frequency is nondimensionalized as `\(\omega/f_o\)`. Inertial peak broadens equatorward. --- class: center ### Comparison with Spectral Fit
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
“Background” fit to anticyclonic (inertial side), matches well the spectrum on the cyclonic (anti-inertial) side, apart from tides. --- class: center ### Wind Stress
Average wind stress and stress magnitude from NCEP 6-hourly product. Gray line separates regions of westward motion (e.g., the Trade Winds) from eastward motion (the subtropical Westerlies). --- class: center ### Inertial Amplitude `\(\sigma\)` from Spectral Fit
In general, higher amplitudes within Westerlies than in Trades, apart from an energetic equatorial band. Northern and southwestern North Pacific also have particularly large amplitudes. --- class: center ### Inertial `\(\sigma\)` weighted by Mixed Layer Depth
Using the mixed layer product of .cite[de Boyer Montégut et. al (2004)]. Quantity plotted is `\(\sigma \times\)` MLD / 20 m. Storm tracks become even more readily apparent. Smaller values and patchier appearance in the Southern Ocean? --- class: center ### Where is Inertial Energy Dominant?
Inertial energy is a large fraction of the total surface current eddy kinetic energy in the North Pacific, and a small fraction over the western boundary currents. --- class: center ### Semidiurnal Amplitude `\(\sigma\)` from Spectral Fit
Semidiurnal energy is highly localized geographically. Not explained by a barotropic tidal model. Equatorial artifact is due to semidiurnal tide falling out of fit band. --- class: center ### M2 Generation in a Numerical Simulation
From .cite[Simmons, Hallberg, and Arbic (2004), Deep Sea Research.] The surface drifters are picking up semidiurnal variability from M2 internal tidal generation as a major feature. --- class: center ### Damping Parameter `\(\lambda\)` from Spectral Fit
Damping timescale does **not** vary with latitude, unlike as in .cite[Park et. al (2009)]. Instead, it presents strong **geographic** variability. Typical damping timescales are `\(1/\lambda=\)`3–5 days. Over strong currents, `\(\approx\)` 1 day, in North Pacific, `\(\approx\)` 10 days. Compare with Fig. 10 of .cite[Elipot et. al (2010)]. --- class: center ### Damping Parameter `\(\lambda\)` without Tidal Model
Artificially inflated damping in M2 internal tide generation regions, as inertial peak tries to broaden to capture semidiurnal variability. --- class: center ### Annual Cycle of Inertial Amplitude
Inertial energy is large when the winds are becoming large, but the mixed layer depths are still shallow, but... --- class: center ### Weighting by Mixed Layer Depth
Weighting by the mixed layer depth aligns the timing of the vertically-integrated inertial energy with the wind stress. --- class: center ### Annual Cycle of Damping
Damping timescale is linked to inertial oscillation amplitude, not vertical integral. Weaker damping for larger amplitudes. --- 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), in prep.] -- * 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), in prep.] -- * Parameter estimation must be corrected for bias due to small sample size effects. .cite[Sykulski, Olhede, Lilly, and Early (2016b), in prep.] --- class: left ## Conclusions, Part II 1. A parametric model fit to the velocity spectra from the global drifter dataset provides a very good match. 1. The M2 tide is a surprisingly strong feature, dominating areas of predicted internal tide generation. 1. Inertial amplitude is larger in regions of Westerlies, and also in the equatorial band. 1. Damping timescales of 3–5 days are seen, with no systematic meridional variability. 1. Damping presents significant and sometimes unexplained geographic variability, suggesting role of stratification. 1. Conjecture that the “inertial peak” in quiescent times and locations actually reflects the internal wave spectrum. --- 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).]]