📖 Why particles can get stuck

In some Parcels simulations, particles end up moving onto parts of the grid where the velocity field is not defined, such as for example onto land. In this tutorial we look at how this depends on the grid structure.

Short conclusion: Particles can get stuck on Arakawa A grids and B grids but should not get stuck on C grids.

This tutorial first looks at how velocity fields are structured and what that means for where the ocean-land boundaries are located. Then we look at how particles may end up getting ‘stuck’ on these boundaries. How to make these particles become ‘unstuck’ again in Parcels will be discussed in this notebook.

• Introduction

• A grid interpolated velocity fields

• B grids

• C grid numerical model - NEMO

• Diffusion

Introduction

Parcels can handle several different types of velocity fields, which makes it widely applicable. This also means that the underlying code and therefore the accuracy of the calculated trajectories can differ, depending on the velocity data input. Even when Parcels runs smoothly with the velocity fields you use, it is good to realize how your velocity fields are structured and how those velocities are generated in the first place.

Horizontal velocity data may be structured on a staggered (Arakawa-C) or unstaggered grid (Arakawa-A). The implementation of these grids is covered in this explainer.

The source of your velocity data will influence how accurate and physically consistent Parcels can calculate trajectories. Common sources are numerical models and data assimilation products, interpolations of those products or discrete observations. The staggering of variables determines how the ocean-land boundaries are defined in models and how their boundary conditions can be satisfied. The Parcels interpolation scheme differs per grid configuration and makes some underlying assumptions that determine the trajectory within a grid cell. Whether the source of your velocity data has physically consistent boundary conditions and how well these align with the assumptions made in Parcels determines how accurately the trajectories move within grid cells. This is especially important near ocean-land boundaries, where particles may end up ‘stuck’ on land.

Here we will look at two examples of velocity fields in Parcels. We visualize the structure of the velocity field, briefly discuss the Parcels implementation and look at how particles get stuck. Then we shortly discuss Arakawa B grids and additional sources of movement that may not respect the ocean-land boundary, such as particle diffusion using a random walk.

from copy import copy
from datetime import timedelta

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import ListedColormap
from matplotlib.lines import Line2D
from scipy import interpolate

import parcels
import parcels.tutorial

1. A grids

Arakawa A grids are unstaggered grids where the velocities \(u\), \(v\) (and \(w\)), pressure and other tracers are defined at the same position (on so-called nodes). In numerical models, these nodes can be interpreted to be located at the corner or at the centre of the grid cells. This means that the cell boundaries, and therefore the solid-fluid boundaries can either be located at the nodes (figure 1A) or at 0.5 dx distance from the nodes (figure 1B) respectively.

Many ocean models natively run on a C grid, because boundary conditions are easier to implement there (see C grid). Sometimes, the C-grid output of these models is interpolated onto an A grid. This is the case for all(?) the data available from Copernicus Marine Data Store, which provide the data on a rectilinear A grid velocity field.

To visualize this, in figure 1 we show the nodes and cells of a coastal region. The ocean cells and nodes are in red and the land cells and nodes are in white.

ds = parcels.tutorial.open_dataset(
    "CopernicusMarine_data_for_stuck_particles_tutorial/data"
)

fields = {"U": ds["uo"], "V": ds["vo"]}
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)

Show code cell source

Hide code cell source

# Define meshgrid coordinates to plot velocity field with matplotlib pcolormesh
dlon = ds["longitude"][1] - ds["longitude"][0]
dlat = ds["latitude"][1] - ds["latitude"][0]

# Outside corner coordinates - coordinates + 0.5 dx
x_outcorners, y_outcorners = np.meshgrid(
    np.append(
        (ds["longitude"] - 0.5 * dlon),
        (ds["longitude"][-1] + 0.5 * dlon),
    ),
    np.append(
        (ds["latitude"] - 0.5 * dlat),
        (ds["latitude"][-1] + 0.5 * dlat),
    ),
)

# Inside corner coordinates - coordinates + 0.5 dx
# needed to plot cells between velocity field nodes
x_incorners, y_incorners = np.meshgrid(
    (ds["longitude"] + 0.5 * dlon)[:-1],
    (ds["latitude"] + 0.5 * dlat)[:-1],
)

# Center coordinates
x_centers, y_centers = np.meshgrid(ds["longitude"], ds["latitude"])

# --------- Velocity fields ---------

# Empty cells between coordinate nodes - essentially on inside corners
cells = np.zeros((len(ds["latitude"]), len(ds["longitude"])))

# # Masking the flowfield where U = NaN
umask = np.ma.masked_invalid(ds["uo"][0, 0, :, :])

# Velocity field with NaN -> zero to be able to use in RectBivariateSpline
u_zeros = ds["uo"][0, 0, :, :].fillna(0)
# Interpolate U
fu = interpolate.RectBivariateSpline(
    ds["latitude"],
    ds["longitude"],
    u_zeros,
    kx=1,
    ky=1,
)

# Velocity field interpolated on the inside corners
u_corners = fu(y_incorners[:, 0], x_incorners[0, :])

# Masking the interpolated flowfield where U = 0
udmask = np.ma.masked_values(u_corners, 0)

# Velocity field with NaN -> zero to be able to use in RectBivariateSpline
v_zeros = ds["vo"][0, 0, :, :].fillna(0)

# Interpolate V
fv = interpolate.RectBivariateSpline(ds["latitude"], ds["longitude"], v_zeros)

# --------- Plotting domain ---------
# Select velocity domain to plot
U = ds["uo"][0, 0, :, :].fillna(0)
V = ds["vo"][0, 0, :, :].fillna(0)

### Create seethrough colormap to show different grid interpretations
cmap = plt.get_cmap("Blues")
my_cmap = cmap(np.arange(cmap.N))
my_cmap[:, -1] = 0  # set alpha to zero
my_cmap = ListedColormap(my_cmap)

Show code cell source

Hide code cell source

fig, axes = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)
fig.suptitle("Figure 1", fontsize=16)

titles = [
    "A) Node at corner cell\nAgreement with Parcels interpolation",
    "B) Node at center cell\nNot according to Parcels interpolation",
    "C) Different boundaries\nParcels interpolates according to orange grid",
]

x_vec = [x_centers, x_outcorners]
y_vec = [y_centers, y_outcorners]
mesh = [udmask.mask, umask.mask]
edgecolors = ["orange", "k"]

for i, ax in enumerate(axes):
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.set_xlim(x_outcorners[0, 0], x_outcorners[0, -1])
    ax.set_ylim(y_outcorners[0, 0], y_outcorners[-1, 0])
    ax.set_title(titles[i])
    ax.pcolormesh(
        x_vec[i % 2],
        y_vec[i % 2],
        mesh[i % 2],
        cmap="Blues_r",
        edgecolors=edgecolors[i % 2],
    )
    if i == 2:
        ax.pcolormesh(
            x_outcorners,
            y_outcorners,
            cells,
            cmap=my_cmap,
            edgecolors="black",
        )
    ax.scatter(
        x_centers,
        y_centers,
        s=50,
        c=U,
        cmap="coolwarm",
        vmin=-0.05,
        vmax=0.05,
        edgecolors="k",
    )
    ax.quiver(
        x_centers,
        y_centers,
        U,
        V,
        angles="xy",
        scale_units="xy",
        scale=2,
        color="w",
        edgecolor="k",
    )

Figure 1 shows how land grid points and boundaries can be interpreted, depending on whether you assume a grid node is at the centre or corner of a cell. When they are at the corner of a cell, boundaries are interpreted as the edges between two nodes. This is especially visible in the upper centre of the figures, where the white nodes with surrounding ocean can be assumed to only have a line of land between them (figure 1A), or as entire land cells (figure 1B).

1.1 Parcels bilinear interpolation

On Arakawa A grids, Parcels uses a simple bilinear interpolation to find the particle velocity at a specific location in a cell. It finds the four nearest velocity components in a 2D field and interpolates between those. This means that the velocity field is essentially divided into cells by Parcels as in figure 1A. The boundaries of cells in this case are located between the nodes of the velocity field and therefore the ocean-land boundary lies in between the nodes of the land. Since both velocity components are defined at all four corner nodes, they equally persist in the limit toward the boundary as they go to zero.

1.2 How do particles get stuck?

Here we create a Parcels simulation of trajectories of particles released near and on the land cells in the velocity field and see how the bilinear interpolation causes particles to get stuck.

npart = 3  # number of particles to be released
lon = np.linspace(4.95, 5.2, npart)
lat = np.linspace(52.95, 53.2, npart)
X, Y = np.meshgrid(lon, lat)

pset = parcels.ParticleSet(fieldset=fieldset, x=X, y=Y)

output_file = parcels.ParticleFile(
    "stuck_particles.parquet",
    outputdt=timedelta(hours=2),
    mode="w",
)

pset.execute(
    parcels.kernels.AdvectionRK2,
    runtime=timedelta(days=3),
    dt=timedelta(minutes=12),
    output_file=output_file,
    verbose_progress=False,
)

INFO: Output files are stored in stuck_particles.parquet

Particles moving toward the boundary will keep slowing down as long as the cross-boundary component of one of the 4 nearest velocity vectors is directed toward the boundary. The distance travelled each dt will decrease to unrealistic values in the absence of local forces directing the flow along the boundary. If we define a particle to be stuck when it moves less than approximately 100 meter every hour, we can see how particles get stuck on the model boundary.

df = parcels.read_particlefile("stuck_particles.parquet")

# Calculating when a particle is stuck
distance = []
for traj in df.partition_by("particle_id", maintain_order=True):
    pdx = np.diff(traj["x"], prepend=traj["x"][0])
    pdy = np.diff(traj["y"], prepend=traj["y"][0])
    distance.append(np.sqrt(np.square(pdx) + np.square(pdy)))
distance = np.array(distance)
stuck = distance < 1e-3

Show code cell source

Hide code cell source

fig, axes = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)
fig.suptitle("Figure 2. CopernicusMarine A-grid", fontsize=16)
gs = gridspec.GridSpec(ncols=2, nrows=1, figure=fig)

titles = [
    "A) Parcels bilinear interpolation",
    "B) Incompatible A grid interpretation",
]

x_vec = [x_centers, x_outcorners]
y_vec = [y_centers, y_outcorners]
mesh = [udmask.mask, umask.mask]
edgecolors = ["orange", "k"]

for i, ax in enumerate(axes):
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.set_xlim(x_outcorners[0, 0], x_outcorners[0, -1])
    ax.set_ylim(y_outcorners[0, 0], y_outcorners[-1, 0])
    ax.set_title(titles[i])
    ax.pcolormesh(
        x_vec[i % 2],
        y_vec[i % 2],
        mesh[i % 2],
        cmap="Blues_r",
        edgecolors=edgecolors[i % 2],
    )
    if i == 2:
        ax.pcolormesh(
            x_outcorners,
            y_outcorners,
            cells,
            cmap=my_cmap,
            edgecolors="black",
        )
    ax.scatter(
        x_centers,
        y_centers,
        s=50,
        c=U,
        cmap="coolwarm",
        vmin=-0.05,
        vmax=0.05,
        edgecolors="k",
    )
    for j, traj in enumerate(df.partition_by("particle_id", maintain_order=True)):
        ax.plot(traj["x"], traj["y"], linewidth=3, zorder=1)
        ax.scatter(traj["x"], traj["y"], c=stuck[j, :], cmap="viridis_r", zorder=2)

color_stuck = copy(plt.get_cmap("viridis"))(0)
color_moving = copy(plt.get_cmap("viridis"))(256)
custom_lines = [
    Line2D([0], [0], c=color_moving, marker="o", markersize=10),
    Line2D([0], [0], c=color_stuck, marker="o", markersize=10),
]
axes[1].legend(
    custom_lines,
    ["Moving", "Stuck"],
    bbox_to_anchor=(1.05, 1),
    loc="upper left",
    borderaxespad=0.0,
    framealpha=1,
)
plt.show()

In figure 2A you can see how particles released in a 3x3 grid keep moving toward the boundary between two land nodes. The ratio of the \(u\) and \(v\) components stays at a similar value due to both being linearly interpolated to zero at the boundary. Note that the interpretation of nodes at the centre of grid cells (figure 2B) is clearly incompatible with Parcels interpolation.

2. B grids

On Arakawa B grids, \(u\) and \(v\) are at the same location, while \(w\) and scalar variables like pressure are staggered. In 2 dimensional flow, Parcels therefore uses the same bilinear interpolation as for A grids to find the horizontal velocity, since the \(u\) and \(v\) components are similarly collocated. This can cause particles to get stuck in the same way as on A grids.

3. C grids

On staggered grids, different types of boundary conditions can be satisfied simultaneously. On a C-grid, the velocities are defined on the cell-edges normal to the velocity-component and pressure, temperature and tracers are defined at the cell centres. This way, a Dirichlet boundary condition can be used for the velocities, while a Neumann boundary condition can be satisfied for the gradient of pressure.

Here we investigate how Parcels interprets the the boundaries on a C grid. For background information, see Delanmeter & Van Sebille (2019). First we show how the velocities are staggered and how the velocity input necessary to create a FieldSet results in the definition of boundaries in Parcels. This example uses a NEMO dataset but many of the relevant C grid assumptions are similar in other models, such as MITgcm.

cufields = parcels.tutorial.open_dataset("NemoNorthSeaORCA025-N006_data/U")
cvfields = parcels.tutorial.open_dataset("NemoNorthSeaORCA025-N006_data/V")
cwfields = parcels.tutorial.open_dataset("NemoNorthSeaORCA025-N006_data/W")
ds_coords = parcels.tutorial.open_dataset("NemoNorthSeaORCA025-N006_data/mesh_mask")[
    ["glamf", "gphif"]
]

xu_corners, yu_corners = np.meshgrid(
    np.arange(cufields["x"].values[0], cufields["x"].values[-1] + 1, 1),
    np.arange(cufields["y"].values[0] - 0.5, cufields["y"].values[-1] + 0.5, 1),
)
xv_corners, yv_corners = np.meshgrid(
    np.arange(cvfields["x"].values[0] - 0.5, cvfields["x"].values[-1] + 0.5, 1),
    np.arange(cvfields["y"].values[0], cvfields["y"].values[-1] + 1, 1),
)
cx_centers, cy_centers = np.meshgrid(
    np.arange(cvfields["x"].values[0] - 0.5, cvfields["x"].values[-1] + 1.5, 1),
    np.arange(cvfields["y"].values[0] - 0.5, cvfields["y"].values[-1] + 1.5, 1),
)
fx_corners, fy_corners = np.meshgrid(
    np.arange(cufields["x"].values[0] - 1, cufields["x"].values[-1] + 1, 1),
    np.arange(cufields["y"].values[0] - 1, cufields["y"].values[-1] + 1, 1),
)
c_cells = np.zeros((len(cufields["y"]), len(cufields["x"])))

# Velocity field with NaN -> zero to be able to use in RectBivariateSpline
cu_zeros = np.nan_to_num(cufields["uos"][0])
f = interpolate.RectBivariateSpline(
    yu_corners[:, 0], xu_corners[0, :], cu_zeros, kx=1, ky=1
)

# Velocity field interpolated on the T-points - center
cu_centers = f(cy_centers[:-1, 0], cx_centers[0, :-1])

# Masking the interpolated flowfield where U = 0
cudmask = np.ma.masked_values(cu_centers, 0)

Show code cell source

Hide code cell source

fig = plt.figure(figsize=(8, 6))
fig.suptitle("Figure 3 - C grid structure", fontsize=16)
ax1 = plt.axes()

ax1.set_xlim(105, 112)
ax1.set_ylim(50, 54)
ax1.pcolormesh(
    fx_corners, fy_corners, cudmask, cmap="Blues", edgecolors="k", linewidth=1
)
ax1.scatter(
    xu_corners,
    yu_corners,
    s=80,
    c=cufields["uos"][0],
    cmap="seismic",
    vmin=-0.1,
    vmax=0.1,
    edgecolor="k",
    label="U",
)
ax1.scatter(
    xv_corners,
    yv_corners,
    s=80,
    c=cvfields["vos"][0],
    cmap="PRGn",
    vmin=-0.1,
    vmax=0.1,
    edgecolor="k",
    label="V",
)
ax1.scatter(cx_centers, cy_centers, s=80, c="orange", edgecolor="k", label="T")
ax1.quiver(
    xu_corners,
    yu_corners,
    cufields["uos"][0],
    np.zeros(xu_corners.shape),
    angles="xy",
    scale_units="xy",
    scale=0.1,
    width=0.007,
)
ax1.quiver(
    xv_corners,
    yv_corners,
    np.zeros(xv_corners.shape),
    cvfields["vos"][0],
    angles="xy",
    scale_units="xy",
    scale=0.3,
    width=0.007,
)

custom_lines = [
    Line2D([0], [0], marker="o", color="r", lw=0),
    Line2D([0], [0], marker="o", color="g", lw=0),
    Line2D([0], [0], marker="o", color="orange", lw=0),
]

ax1.legend(
    custom_lines,
    ["U", "V", "T"],
    bbox_to_anchor=(1.05, 1),
    loc="upper left",
    borderaxespad=0.0,
);

In figure 3 you can see how a boundary can be traced along the cell edges, where the velocity is zero, because the normal velocities are defined at the edges in a C-grid. This ensures that the most important boundary condition in many models is satisfied: the normal-component of velocity is zero at the boundary. Within a cell, Parcels interpolates the velocities on a C-grid linearly in the direction of that component and constant throughout a cell in the other directions. This means that the cross-boundary component is linearly interpolated to zero at the boundary and the along-boundary component is constant towards the boundary. Note that this is the same interpolation as used in the ‘Analytical Advection’ Scheme of e.g. TRACMASS and Ariane. TODO add link to Analytical Advection Scheme.

3.1 Consistency with NEMO configuration

Here we look at how the lateral boundary conditions are implemented in NEMO. This is documented on page 23 of the NEMO documentation.

The cross-boundary component of velocity, defined at the cell faces where the boundary is defined, is set to zero at all boundaries. See the figure below. This corresponds exactly with the velocity field in a Parcels FieldSet shown in figure 3.

There are different options for the along-boundary velocity in NEMO: a free-slip condition, a partial-slip condition and a no-slip condition. Since the tangential velocity is not defined at the boundary, this boundary condition is defined by a Neumann boundary condition: the normal derivative of the along-boundary velocity is specified. This derivative is schematically represented by a “ghost” velocity on the adjacent land node. The specified derivative is equivalent to what would result from the central difference between the along-boundary velocity at the nearest ocean cell and this “ghost” velocity. The type of boundary condition defines the direction and magnitude of this “ghost” velocity relative to the along-boundary velocity in the fluid domain. In Parcels, these “ghost” velocities may be used to determine how the velocity should be interpolated near the coast. By default, Parcels interpolates piecewise-constant in the direction normal to the velocity component. This means that the along-boundary velocity is the same for any distance away from the boundary and therefore equivalent to the free-slip boundary condition shown in subfigure (a) below.

3.2 Particle trajectories near the boundary in a C grid

Now let’s see how particles move near the boundary in a FieldSet derived from the mass-conserving NEMO C grid.

ds_fset = parcels.convert.nemo_to_sgrid(
    fields={"U": cufields["uo"], "V": cvfields["vo"], "W": cwfields["wo"]},
    coords=ds_coords,
)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
fieldset = fieldset.to_windowed_arrays()

npart = 10  # number of particles to be released
lon = np.linspace(3, 4, npart)
lat = 51.5 * np.ones(npart)

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

output_file = parcels.ParticleFile(
    "Cgrid-stuck.parquet",
    outputdt=timedelta(hours=1),
    mode="w",
)

pset.execute(
    parcels.kernels.AdvectionRK2,
    runtime=timedelta(days=10),
    dt=timedelta(minutes=5),
    output_file=output_file,
    verbose_progress=False,
)

INFO: Output files are stored in Cgrid-stuck.parquet

df_C = parcels.read_particlefile("Cgrid-stuck.parquet")

distance_C = []
for traj in df_C.partition_by("particle_id", maintain_order=True):
    pdx = np.diff(traj["x"], prepend=traj["x"][0])
    pdy = np.diff(traj["y"], prepend=traj["y"][0])
    distance_C.append(np.sqrt(np.square(pdx) + np.square(pdy)))
distance_C = np.array(distance_C)
stuck_C = distance_C < 1e-3

Show code cell source

Hide code cell source

fig, axes = plt.subplots(1, 2, figsize=(14, 6), constrained_layout=True)
fig.suptitle(
    "Figure 4. - NEMO model movement along solid cells in a C grid", fontsize=16
)

xlim = ((-5, 5), (2.9, 4.1))
ylim = ((48, 58), (51.4, 52.2))
titles = ["North Sea domain", "North Sea domain - zoomed in"]
edgecolors = [None, "k"]
for i, ax in enumerate(axes):
    ax.set_ylabel("Latitude [degrees]")
    ax.set_xlabel("Longitude [degrees]")
    ax.set_title(titles[i])
    ax.set_xlim(xlim[i])
    ax.set_ylim(ylim[i])

    ax.pcolormesh(
        ds_coords["glamf"],
        ds_coords["gphif"],
        cvfields["vo"][0, 0, 1:, 1:],
        cmap="RdBu",
        vmin=-0.1,
        vmax=0.1,
        edgecolors=edgecolors[i],
    )
    for i, traj in enumerate(df_C.partition_by("particle_id", maintain_order=True)):
        ax.plot(traj["x"], traj["y"], linewidth=3, zorder=1)
        ax.scatter(
            traj["x"],
            traj["y"],
            c=stuck_C[i, :],
            cmap="viridis_r",
            zorder=2,
            vmin=0,
            vmax=1,
        )

color_stuck = copy(plt.get_cmap("viridis"))(0)
color_moving = copy(plt.get_cmap("viridis"))(256)
custom_lines = [
    Line2D([0], [0], c=color_moving, marker="o", markersize=10),
    Line2D([0], [0], c=color_stuck, marker="o", markersize=10),
]
axes[1].legend(
    custom_lines,
    ["Moving", "Stuck"],
    bbox_to_anchor=(1.05, 1),
    loc="upper left",
    borderaxespad=0.0,
    framealpha=1,
)
plt.show()

As we can see in figure 4, the particles that start in the ocean move along the boundary, since the cross-boundary component goes to zero, but the along-boundary component is equal to the value away from the coast.

4. Diffusion

In many Parcels simulations other motions are added. Examples are Brownian motion from a diffusivity field and Stokes drift based on a wind field. These sources of motion will not necessarily consider the same boundaries and, unless they are also based on a C grid, they may also result in particles getting stuck on land. To show how this happens, let us add a random walk to the advection in a C grid velocity field.

This random walk can be added using a diffusion kernel, as documented in this notebook. Since the particles will move randomly through the domain, without awareness of the solid-fluid boundaries in the velocity field, we cannot define stuck particles as having moved less than a tolerance value and we will instead check whether particles find themselves on land or not. To do this, we sample the local velocities.

LandParticle = parcels.Particle.add_variable(parcels.Variable("on_land"))


def SampleLand(particles, fieldset):
    u, v = fieldset.UV[particles]
    particles.on_land = (u == 0) & (v == 0)

fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
fieldset = fieldset.to_windowed_arrays()

fieldset.add_constant_field("Kh_zonal", 5)
fieldset.add_constant_field("Kh_meridional", 5)
fieldset.add_context("dres", 0.083)

npart = 10  # number of particles to be released
lon = np.linspace(3, 4, npart)
lat = 51.5 * np.ones(npart)

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

output_file = parcels.ParticleFile(
    "Cgrid-diffusion.parquet",
    outputdt=timedelta(hours=1),
    mode="w",
)

pset.execute(
    [parcels.kernels.AdvectionDiffusionEM, SampleLand],
    runtime=timedelta(days=10),
    dt=timedelta(minutes=5),
    output_file=output_file,
    verbose_progress=False,
)

INFO: Output files are stored in Cgrid-diffusion.parquet

df_C_diff = parcels.read_particlefile("Cgrid-diffusion.parquet")

Show code cell source

Hide code cell source

fig = plt.figure(figsize=(14, 8))
fig.suptitle("Figure 5. Diffusion added to C grid advection", fontsize=16)
ax = plt.axes()

ax.set_ylabel("Latitude [degrees]")
ax.set_xlabel("Longitude [degrees]")
ax.set_title("Particle getting stuck")
ax.set_xlim(2.9, 4.1)
ax.set_ylim(51.4, 52.2)

ax.pcolormesh(
    ds_coords["glamf"],
    ds_coords["gphif"],
    cvfields["vo"][0, 0, 1:, 1:],
    cmap="RdBu",
    vmin=-0.1,
    vmax=0.1,
    edgecolors="k",
    linewidth=1,
)
for traj in df_C_diff.partition_by("particle_id", maintain_order=True):
    ax.plot(traj["x"], traj["y"], linewidth=3, zorder=1)
    ax.scatter(
        traj["x"],
        traj["y"],
        c=traj["on_land"],
        cmap="viridis_r",
        zorder=2,
        vmin=0,
        vmax=1,
    )

color_stuck = copy(plt.get_cmap("viridis"))(0)
color_moving = copy(plt.get_cmap("viridis"))(256)

custom_lines = [
    Line2D([0], [0], c=color_stuck, marker="o", markersize=10, lw=0),
    Line2D([0], [0], c=color_moving, marker="o", markersize=10, lw=0),
]
ax.legend(
    custom_lines,
    ["On land", "In ocean"],
    bbox_to_anchor=(1.05, 1),
    loc="upper left",
    borderaxespad=0.0,
    framealpha=1,
)
plt.show()

Figure 5 shows why you should consider additional boundary conditions for each component of the motion added to the particle.