Cabo Verde hyperbolic LCS¶

The FTLE map (see cabo_verde_ftle) shows where the flow stretches. It does not give the material curves along which the stretching happens. Those curves are tensor lines of the strain tensor. Haller (2015, §5.1 / Table 1, doi:10.1146/annurev-fluid-010313-141322) constructs those curves directly from the Cauchy–Green strain tensor \(C = (\nabla F)^\top \nabla F\), whose eigenpairs satisfy \(C\,\xi_i = \lambda_i\,\xi_i\) with \(0 < \lambda_1 \le \lambda_2\) and \(\xi_1 \perp \xi_2\). A repelling LCS is a shrink line — a curve tangent to \(\xi_1\), i.e. orthogonal to the strong-stretch direction \(\xi_2\) that the FTLE ridge marks. It solves the ODE \(\dot r = \xi_1(r)\).

Attracting LCS come from the forward–backward duality (Haller & Sapsis 2011, doi:10.1063/1.3579597): an attracting LCS is a repelling LCS of the backward flow. The notebook therefore advects the same seed grid both ways and takes the \(\xi_1\) shrink lines of each flow map: the forward map gives the repelling LCS, the backward map the attracting ones.

# Importing Parcels pulls in the holoviews/bokeh bootstrap and prints an
# alpha-version notice, neither of which belongs on the rendered page.
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from parcels import FieldSet, Particle, ParticleSet, StatusCode
from parcels.convert import copernicusmarine_to_sgrid
from parcels.kernels import AdvectionRK4

from lcs_parcels import (
    NeighborSeedGrid,
    ftle_ridge_seeds,
    prune_shrink_lines,
    shrink_lines,
)

Parcels v4 field set¶

copernicusmarine_to_sgrid tags the CMEMS A-grid with SGRID metadata; from_sgrid_conventions wraps it as a spherical FieldSet.

sgrid = copernicusmarine_to_sgrid(fields={"U": currents["uo"], "V": currents["vo"]})
fieldset = FieldSet.from_sgrid_conventions(sgrid, mesh="spherical")
z_surface = float(currents["depth"].values[0])

Seed grid and window¶

\(t_0\) sits at the middle of the window the local file covers, and the runs are \(\pm 5\) d, so the forward and the backward advection both stay inside the data. A rectilinear NeighborSeedGrid puts one particle on every diagnostic grid point; \(\nabla F\) is then differenced against the grid neighbours.

t0 = np.datetime64("2025-08-06")
T = np.timedelta64(5, "D")
resolution_deg = 1 / 25
seed_lon, seed_lat = (-27.0, -21.0), (13.5, 18.5)

lon_axis = np.arange(seed_lon[0], seed_lon[1] + 1e-9, resolution_deg)
lat_axis = np.arange(seed_lat[0], seed_lat[1] + 1e-9, resolution_deg)
seed = NeighborSeedGrid.from_axes(lon=lon_axis, lat=lat_axis)
seed

<NeighborSeedGrid 151x126 grid, lon -27.00..-21.00, lat 13.50..18.50>

Recovery kernel¶

Particles that leave the domain or hit land are turned into NaN in place (Parcels would otherwise abort the run), so losses propagate as NaN through every diagnostic below. StatusCode.EndofLoop rather than StatusCode.Delete: deleting shrinks the particle array and breaks the alignment with the seed order.

def set_lost_to_nan(particles, fieldset):
    lost = particles.state >= StatusCode.Error
    particles.x = np.where(lost, np.nan, particles.x)
    particles.y = np.where(lost, np.nan, particles.y)
    particles.state = np.where(lost, StatusCode.EndofLoop, particles.state)

Advect forward¶

to_parcels_pset emits the flat (lon, lat) pair; the finals go back in via pset_to_flowmap, which takes both times and derives the signed window \(T = t_1 - t_0\).

lon, lat = seed.to_parcels_pset()
z = np.full(len(lon), z_surface)
pset = ParticleSet(fieldset, pclass=Particle, x=lon, y=lat, z=z, t=t0)
pset.execute(
    [AdvectionRK4, set_lost_to_nan],
    dt=np.timedelta64(1, "h"),
    runtime=T,
    verbose_progress=False,
)
forward = seed.pset_to_flowmap(lon=pset.x, lat=pset.y, t0=t0, t1=t0 + T)
forward

<NeighborFlowMap 151x126 grid, lon -27.00..-21.00, lat 13.50..18.50,
                 t0 2025-08-06T00:00:00, T +5.0 days>

Advect backward¶

The same seed grid, same \(t_0\), negative dt, and \(t_1 = t_0 - T\) so the stored window is negative.

lon, lat = seed.to_parcels_pset()
z = np.full(len(lon), z_surface)
pset = ParticleSet(fieldset, pclass=Particle, x=lon, y=lat, z=z, t=t0)
pset.execute(
    [AdvectionRK4, set_lost_to_nan],
    dt=-np.timedelta64(1, "h"),
    runtime=T,
    verbose_progress=False,
)
backward = seed.pset_to_flowmap(lon=pset.x, lat=pset.y, t0=t0, t1=t0 - T)
backward

<NeighborFlowMap 151x126 grid, lon -27.00..-21.00, lat 13.50..18.50,
                 t0 2025-08-06T00:00:00, T -5.0 days>

Repelling LCS over the forward FTLE¶

The pruned set of curves. The seed points are dotted, which shows whether they sit on the ridge tops and stay separated. The FTLE gets a greyscale so the coloured curves stand out against it. The dots with no curve through them are the untraceable seeds and the seeds whose line was pruned.

fig, ax = plt.subplots()
ftle_forward.plot.pcolormesh(x="lon_grid", y="lat_grid", ax=ax, cmap="Greys")
ax.plot(
    repelling_lcs_pruned["lon"].values.T,
    repelling_lcs_pruned["lat"].values.T,
    color="tab:red",
    lw=0.8,
)
ax.scatter(repelling_seeds["lon"], repelling_seeds["lat"], s=8, color="tab:red")
plt.show()

../_images/8ecd256317f7c7f64fb43e63f612345f9f3563f975c77a3248031f937e18f10d.png

Attracting LCS over the backward FTLE¶

Same picture for the backward flow map, straight out of the one dataset hyperbolic_lcs returned.

fig, ax = plt.subplots()
attracting_lcs["ftle"].plot.pcolormesh(x="lon_grid", y="lat_grid", ax=ax, cmap="Greys")
ax.plot(
    attracting_lcs["lon"].values.T,
    attracting_lcs["lat"].values.T,
    color="tab:blue",
    lw=0.8,
)
plt.show()

../_images/e556fb00ed7ab771cf5a68fc5947e9b80db54599e5954c31957f0da464a540de.png

Both families together¶

Both families over the forward FTLE in one frame. Where a repelling and an attracting curve cross, the flow is locally hyperbolic.

fig, ax = plt.subplots()
ftle_forward.plot.pcolormesh(x="lon_grid", y="lat_grid", ax=ax, cmap="Greys")
ax.plot(
    repelling_lcs_pruned["lon"].values.T,
    repelling_lcs_pruned["lat"].values.T,
    color="tab:red",
    lw=0.8,
)
ax.plot(
    attracting_lcs["lon"].values.T,
    attracting_lcs["lat"].values.T,
    color="tab:blue",
    lw=0.8,
)
plt.show()

../_images/09a4460e384c9f83f0f1bd5c94057d38e999c8681d0bd452886108b8ed2325ca.png

Cabo Verde hyperbolic LCS¶

Currents¶

Parcels v4 field set¶

Seed grid and window¶

Recovery kernel¶

Advect forward¶

Advect backward¶

Repelling LCS, step by step¶

Pruning the bundles¶

Attracting LCS, in one call¶

Repelling LCS over the forward FTLE¶

Attracting LCS over the backward FTLE¶

Both families together¶

LCS-Parcels

Navigation

Related Topics