🖥️ FESOM tutorial#

Parcels v4 supports unstructured-grid model output via uxarray. This tutorial walks through the minimum steps to advect particles in real FESOM2 output. The recipe is:

  1. Open the FESOM grid and data files with uxarray.

  2. Rename FESOM-specific dimensions to Parcels’ UGRID conventions with parcels.convert.fesom_to_ugrid().

  3. Build a parcels.FieldSet with parcels.FieldSet.from_ugrid_conventions().

  4. Run the simulation as on any structured grid.

If you have not done so already, work through the quickstart tutorial first to get familiar with ParticleSet, Kernel, and ParticleFile.

from pathlib import Path

import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np
import uxarray as ux

import parcels
import parcels.tutorial

Get the FESOM tutorial dataset#

We use a small periodic-channel snapshot from a FESOM2 simulation that ships with Parcels’ tutorial data registry. As in the quickstart, parcels.tutorial.open_dataset() downloads the files into a local cache on first use; subsequent calls just return the cached copy.

uxarray expects file paths rather than an in-memory dataset, so we trigger the downloads and then point ux.open_mfdataset at the cached files:

for name in [
    "FESOM_periodic_channel/fesom_channel",  # grid description
    "FESOM_periodic_channel/u.fesom_channel",  # zonal velocity (face-registered)
    "FESOM_periodic_channel/v.fesom_channel",  # meridional velocity (face-registered)
    "FESOM_periodic_channel/w.fesom_channel",  # vertical velocity (node-registered)
]:
    parcels.tutorial.open_dataset(name)

from parcels._datasets.remote import _DATA_HOME

data_dir = Path(_DATA_HOME) / "data" / "FESOM_periodic_channel"

grid_path = str(data_dir / "fesom_channel.nc")
data_paths = [
    str(data_dir / "u.fesom_channel.nc"),
    str(data_dir / "v.fesom_channel.nc"),
    str(data_dir / "w.fesom_channel.nc"),
]

Note

The FESOM_periodic_channel dataset is a single-time-step snapshot of an idealised channel (~4.4° × ~18° wide, 40 vertical layers). Working with multi-time FESOM output is identical, except you pass a glob or list of data files to ux.open_mfdataset.

Open the data with uxarray#

ux.open_mfdataset(grid_path, data_paths) reads the FESOM grid description and joins the velocity files on its grid. FESOM names its velocity variables u, v, w — we rename them to U, V, W so that parcels.FieldSet.from_ugrid_conventions recognises them as the velocity components:

ds = ux.open_mfdataset(grid_path, data_paths).rename_vars(
    {"u": "U", "v": "V", "w": "W"}
)
ds

<xarray.UxDataset> Size: 8MB
Dimensions:  (nz1: 40, time: 1, n_face: 21120, nz: 41, n_node: 10608)
Coordinates:
  * nz1      (nz1) float64 320B 4.519 14.01 24.45 ... 3.459e+03 3.814e+03
  * time     (time) datetime64[ns] 8B 2005-01-01T23:50:00
  * nz       (nz) float64 328B 0.0 9.038 18.98 ... 3.29e+03 3.628e+03 4e+03
Dimensions without coordinates: n_face, n_node
Data variables:
    U        (time, nz1, n_face) float32 3MB dask.array<chunksize=(1, 40, 21120), meta=np.ndarray>
    V        (time, nz1, n_face) float32 3MB dask.array<chunksize=(1, 40, 21120), meta=np.ndarray>
    W        (time, nz, n_node) float32 2MB dask.array<chunksize=(1, 41, 10608), meta=np.ndarray>
Attributes: (12/20)
    Conventions:                         UGRID-1.0
    FESOM_model:                         FESOM2
    FESOM_website:                       fesom.de
    FESOM_git_SHA:                       ab6d892b
    FESOM_MeshPath:                      /albedo/work/projects/p_clidyn_work/...
    FESOM_mesh_representative_checksum:  
    ...                                  ...
    FESOM_evp_rheol_steps:               120
    FESOM_opt_visc:                      5
    FESOM_use_wsplit:                    0
    FESOM_autorotate_back_to_geo:        0
    history:                             Tue Dec  3 11:56:05 2024: ncks -d ti...
    NCO:                                 netCDF Operators version 5.2.9 (Home...

Note the FESOM-specific dimension names: elem (number of triangular faces), nod2 (number of nodes), nz1 (vertical layer centres), and nz (layer interfaces). The horizontal velocities U and V live on face centres along nz1; the vertical velocity W lives on nodes along nz.

Convert to UGRID conventions#

Parcels works with a small UGRID-compliant dialect: nodes are n_node, faces are n_face, vertical centres are zc, and vertical interfaces are zf. The helper parcels.convert.fesom_to_ugrid does this rename in one call:

ds = parcels.convert.fesom_to_ugrid(ds)
print("dims:", dict(ds.sizes))

INFO: Using known vertical dimension mapping: 'nz' (interfaces) and 'nz1' (centers).
INFO: Renaming vertical dimensions: {'nz': 'zf', 'nz1': 'zc'}
dims: {'zc': 40, 'time': 1, 'n_face': 21120, 'zf': 41, 'n_node': 10608}

Build the FieldSet#

With UGRID-compliant dimensions in place, parcels.FieldSet.from_ugrid_conventions() builds the FieldSet. It detects U, V, W, assigns a UxGrid to each field, and picks an appropriate interpolator based on each variable’s coordinate location (face- vs. node-registered, centre vs. interface).

fieldset = parcels.FieldSet.from_ugrid_conventions(ds)

fieldset.describe()

| Name   | Type        |   Grid number | Interp method / value         | Parcels backend   |
|:-------|:------------|--------------:|:------------------------------|:------------------|
| U      | Field       |             0 | UxConstantFaceConstantZC(...) | Dask              |
| V      | Field       |             0 | UxConstantFaceConstantZC(...) | Dask              |
| W      | Field       |             0 | UxLinearNodeLinearZF(...)     | Dask              |
| UV     | VectorField |             0 | Ux_Velocity(...)              | -                 |
| UVW    | VectorField |             0 | Ux_Velocity(...)              | -                 |

mesh: SphericalMesh(radius=6366707.019493707)
time interval: None

U and V get face-registered interpolation, W gets node-registered linear interpolation. The combined vector fields UV and UVW are assembled automatically.

Release particles and advect#

We seed particles on a grid of four latitudes spanning the channel and ten longitudes, and integrate for two days with RK4. Because this snapshot has only a single time level, fieldset.time_interval is None and we omit the t= argument so that Parcels treats the flow as constant in time:

lon_grid, lat_grid = np.meshgrid(
    np.linspace(0.5, 4.0, 10),
    np.linspace(3.0, 15.0, 4),
)
lon = lon_grid.ravel()
lat = lat_grid.ravel()
z = np.full(lon.size, 50.0)  # release at 50 m depth

pset = parcels.ParticleSet(
    fieldset=fieldset,
    pclass=parcels.Particle,
    x=lon,
    y=lat,
    z=z,
)

output_file = parcels.ParticleFile(
    "output-fesom.parquet", outputdt=np.timedelta64(1, "h")
)

pset.execute(
    [parcels.kernels.AdvectionRK4],
    runtime=np.timedelta64(2, "D"),
    dt=np.timedelta64(5, "m"),
    output_file=output_file,
    verbose_progress=False,
)

INFO: Output files are stored in output-fesom.parquet

Plot the velocity field and trajectories#

We plot the particle paths on top of the velocity field they advect through: triangle colour shows the speed at the release depth (≈50 m), black arrows show the velocity at face centres (drawn in lon/lat space, length proportional to speed), grey lines trace each particle’s path, and the coloured dots mark the positions over time. Drawing the arrows with angles="xy" keeps them aligned with the trajectories, so you can see the particles streak along the fast jets and barely move in the quiet bands between them:

df = parcels.read_particlefile("output-fesom.parquet")

triang = mtri.Triangulation(
    ds.uxgrid.node_lon.values,
    ds.uxgrid.node_lat.values,
    triangles=ds.uxgrid.face_node_connectivity.values,
)

depth_idx = int(np.argmin(np.abs(ds.zc.values - 50.0)))
U_face = np.asarray(ds["U"].isel(zc=depth_idx)).squeeze()
V_face = np.asarray(ds["V"].isel(zc=depth_idx)).squeeze()
speed = np.hypot(U_face, V_face)

fig, ax = plt.subplots(figsize=(11, 5))

# Background: speed at the release depth.
tpc = ax.tripcolor(triang, facecolors=speed, shading="flat", cmap="Blues")
fig.colorbar(tpc, ax=ax, label="speed [m/s]", location="left", shrink=0.85)

# Velocity arrows at face centres. Drawing in lon/lat space (angles/scale_units
# "xy") keeps the arrows aligned with the trajectories; length tracks speed.
step = max(1, U_face.size // 400)
xf = ds.uxgrid.face_lon.values
yf = ds.uxgrid.face_lat.values
max_speed = float(np.nanmax(speed))
q = ax.quiver(
    xf[::step],
    yf[::step],
    U_face[::step],
    V_face[::step],
    angles="xy",
    scale_units="xy",
    scale=max_speed / 0.3,
    color="k",
    width=0.0018,
    pivot="tail",
)
ax.quiverkey(
    q, 0.86, 1.04, max_speed, f"{max_speed:.2f} m/s", labelpos="E", coordinates="axes"
)

# Particle paths (grey lines) and positions coloured by time.
for traj in df.sort("t").partition_by("particle_id"):
    ax.plot(traj["x"], traj["y"], color="0.4", linewidth=0.6, alpha=0.7, zorder=2)
ax.scatter(lon, lat, facecolors="none", edgecolors="k", s=60, label="release", zorder=3)
sc = ax.scatter(
    df["x"],
    df["y"],
    c=df["t"].dt.total_seconds(),
    s=6,
    cmap="viridis",
    zorder=3,
)
fig.colorbar(sc, ax=ax, label="time since release [s]", shrink=0.85)

ax.set_xlabel("Longitude [deg E]")
ax.set_ylabel("Latitude [deg N]")
ax.set_title(
    f"FESOM2 velocity at z ≈ {float(ds.zc.values[depth_idx]):.1f} m with particle trajectories"
)
ax.legend(loc="upper right")
plt.show()
../../_images/8598c829bbe3c341ae93ca7fd863b793db2c1a87da3ac391c0cc8c2b92d0414b.png

The particles drift through the channel following the FESOM2 velocity field. From here, the rest of Parcels — custom kernels, sampling fields onto particles, writing your own interpolators — works identically to structured grids. See the interpolation tutorial for the available Ux* interpolators and how to write a custom one.