📖 Preventing stuck particles

In another notebook, we have shown how particles may end up getting stuck on land, especially in A gridded velocity fields. Here we show how you can work around this problem and how large the effects of the solutions on the trajectories are.

Common solutions are:

1. Delete the particles

2. Displace the particles when they are within a certain distance of the coast.

3. Implement free-slip or partial-slip boundary conditions

In the first two of these solutions, kernels are used to modify the trajectories near the coast. The kernels all consist of two parts:

1. Flag particles whose trajectory should be modified

2. Modify the trajectory accordingly

In the third solution, the interpolation method is changed; this has to be done when creating the FieldSet.

This notebook is mainly focused on comparing the different modifications to the trajectory. The flagging of particles is also very relevant however and further discussion on this is encouraged. Some options shown here are:

1. Flag particles within a specific distance to the shore

2. Flag particles in any gridcell that has a shore edge

As argued in the previous notebook, it is important to accurately plot the grid discretization, in order to understand the motion of particles near the boundary. The velocity fields can best be depicted using points or arrows that define the velocity at a single position. Four of these nodes then form gridcells that can be shown using tiles, for example with matplotlib.pyplot.pcolormesh.

from datetime import timedelta

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from scipy import interpolate

import parcels
import parcels.tutorial

1. Particle deletion

The simplest way to avoid trajectories that interact with the coastline is to remove them entirely. You can do this for example by a Kernel that samples the velocity field at the particle position and checks whether the particle is on land or not. If it is, the particle is deleted. This is shown in the following example.

def DeleteOnLand(particles, fieldset):
    u, v = fieldset.UV[particles]
    on_land = (u == 0) & (v == 0)
    particles[on_land].state = parcels.StatusCode.Delete

2. Displacement

A simple concept to avoid particles moving onto shore is displacing them towards the ocean as they get close to shore. This is for example done in Kaandorp et al. (2020) and Delandmeter and van Sebille (2018). To do so, a particle must be ‘aware’ of where the shore is and displaced accordingly. In Parcels, we can do this by adding a ‘displacement’ Field to the Fieldset, which contains vectors pointing away from shore.

Step 1 Import a velocity field - the A gridded CopernicusMarine velocity field

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)
fieldset = fieldset.to_windowed_arrays()

# 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"])

Step 2: Make a landmask where land = 1 and ocean = 0.

# Interpolate the landmask to the cell centers
# only cells with 4 neighboring land points will be land
landmask = np.ma.masked_invalid(ds["uo"][0, 0, :, :], 0)
landmask[landmask != 0] = 1

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

l_corners = fl(y_incorners[:, 0], x_incorners[0, :])

# land when interpolated value == 1
lmask = np.ma.masked_values(l_corners, 0)
lmask[lmask != 0] = 1

Show code cell source

Hide code cell source

fig, axes = plt.subplots(1, 2, figsize=(10, 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",
]

x_vec = [x_centers, x_outcorners]
y_vec = [y_centers, y_outcorners]
mesh = [lmask, landmask]
edgecolors = ["orange", "k"]

for i, ax in enumerate(axes):
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.set_xlim(x_incorners[0, 0], x_incorners[0, -1])
    ax.set_ylim(y_incorners[0, 0], y_incorners[-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],
    )
    ax.scatter(
        x_centers,
        y_centers,
        c=landmask.mask.astype(int),
        s=50,
        cmap="coolwarm",
        vmin=-1,
        vmax=1,
        edgecolors="k",
    )

Figure 1 shows why it is important to be precise when visualizing the model land and ocean. Parcels trajectories should not cross the land boundary between two land nodes as seen in 1B.

Step 3: Detect the coast

We can detect the edges between land and ocean nodes by computing the Laplacian with the 4 nearest neighbors [i+1,j], [i-1,j], [i,j+1] and [i,j-1]:

\[\nabla^2 \text{landmask} = \partial_{xx} \text{landmask} + \partial_{yy} \text{landmask},\]

and filtering the positive and negative values. This gives us the location of coast nodes (ocean nodes next to land) and shore nodes (land nodes next to the ocean).

Additionally, we can find the nodes that border the coast/shore diagonally by considering the 8 nearest neighbors, including [i+1,j+1], [i-1,j+1], [i-1,j+1] and [i-1,j-1].

def get_coastal_nodes(landmask):
    """Function that detects the coastal nodes, i.e. the ocean nodes directly
    next to land. Computes the Laplacian of landmask.

    - landmask: the land mask where land cell = 1 and ocean cell = 0.

    Output: 2D array array containing the coastal nodes, the coastal nodes are
            equal to one, and the rest is zero.
    """
    mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
    mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
    mask_lap -= 4 * landmask
    coastal = np.ma.masked_array(landmask, mask_lap > 0)
    coastal = coastal.mask.astype("int")

    return coastal


def get_shore_nodes(landmask):
    """Function that detects the shore nodes, i.e. the land nodes directly
    next to the ocean. Computes the Laplacian of landmask.

    - landmask: the land mask where land cell = 1 and ocean cell = 0.

    Output: 2D array array containing the shore nodes, the shore nodes are
            equal to one, and the rest is zero.
    """
    mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
    mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
    mask_lap -= 4 * landmask
    shore = np.ma.masked_array(landmask, mask_lap < 0)
    shore = shore.mask.astype("int")

    return shore

def get_coastal_nodes_diagonal(landmask):
    """Function that detects the coastal nodes, i.e. the ocean nodes where
    one of the 8 nearest nodes is land. Computes the Laplacian of landmask
    and the Laplacian of the 45 degree rotated landmask.

    - landmask: the land mask built using `make_landmask`, where land cell = 1
                and ocean cell = 0.

    Output: 2D array array containing the coastal nodes, the coastal nodes are
            equal to one, and the rest is zero.
    """
    mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
    mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
    mask_lap += np.roll(landmask, (-1, 1), axis=(0, 1)) + np.roll(
        landmask, (1, 1), axis=(0, 1)
    )
    mask_lap += np.roll(landmask, (-1, -1), axis=(0, 1)) + np.roll(
        landmask, (1, -1), axis=(0, 1)
    )
    mask_lap -= 8 * landmask
    coastal = np.ma.masked_array(landmask, mask_lap > 0)
    coastal = coastal.mask.astype("int")

    return coastal


def get_shore_nodes_diagonal(landmask):
    """Function that detects the shore nodes, i.e. the land nodes where
    one of the 8 nearest nodes is ocean. Computes the Laplacian of landmask
    and the Laplacian of the 45 degree rotated landmask.

    - landmask: the land mask built using `make_landmask`, where land cell = 1
                and ocean cell = 0.

    Output: 2D array array containing the shore nodes, the shore nodes are
            equal to one, and the rest is zero.
    """
    mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
    mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
    mask_lap += np.roll(landmask, (-1, 1), axis=(0, 1)) + np.roll(
        landmask, (1, 1), axis=(0, 1)
    )
    mask_lap += np.roll(landmask, (-1, -1), axis=(0, 1)) + np.roll(
        landmask, (1, -1), axis=(0, 1)
    )
    mask_lap -= 8 * landmask
    shore = np.ma.masked_array(landmask, mask_lap < 0)
    shore = shore.mask.astype("int")

    return shore

coastal = get_coastal_nodes_diagonal(landmask)
shore = get_shore_nodes_diagonal(landmask)

Show code cell source

Hide code cell source

fig, axes = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)
fig.suptitle("Figure 2. Coast and shore", fontsize=16)

titles = [
    "A) Coastal points",
    "B) Shore points",
]

cmap = [coastal, shore]
for i, ax in enumerate(axes):
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    ax.set_xlim(x_incorners[0, 0], x_incorners[0, -1])
    ax.set_ylim(y_incorners[0, 0], y_incorners[-1, 0])
    ax.set_title(titles[i])
    ax.pcolormesh(
        x_centers,
        y_centers,
        lmask.mask,
        cmap="Blues_r",
    )
    inds = np.where(cmap[i] == 1)
    ax.plot(
        x_centers[inds],
        y_centers[inds],
        "o",
        color="y",
        markersize=10,
        markeredgecolor="k",
    )

Step 4: Assigning coastal velocities

For the displacement kernel we define a velocity field that pushes the particles back to the ocean. This velocity is a vector normal to the shore.

For the shore nodes directly next to the ocean, we can take the simple derivative of landmask and project the result to the shore array, this will capture the orientation of the velocity vectors.

For the shore nodes that only have a diagonal component, we need to take into account the diagonal nodes also and project the vectors only onto the inside corners that border the ocean diagonally.

Then to make the vectors unitary, we normalise them by their magnitude.

def create_displacement_field(landmask, double_cell=False):
    """Function that creates a displacement field 1 m/s away from the shore.

    - landmask: the land mask
    - double_cell: Boolean for determining if you want a double cell.
      Default set to False.

    Output: two 2D arrays, one for each component of the velocity.
    """
    shore = get_shore_nodes(landmask)

    # nodes bordering ocean directly and diagonally
    shore_d = get_shore_nodes_diagonal(landmask)
    # corner nodes that only border ocean diagonally
    shore_c = shore_d - shore

    # Simple derivative
    Ly = np.roll(landmask, -1, axis=0) - np.roll(landmask, 1, axis=0)
    Lx = np.roll(landmask, -1, axis=1) - np.roll(landmask, 1, axis=1)

    Ly_c = np.roll(landmask, -1, axis=0) - np.roll(landmask, 1, axis=0)
    # Include y-component of diagonal neighbors
    Ly_c += np.roll(landmask, (-1, -1), axis=(0, 1)) + np.roll(
        landmask, (-1, 1), axis=(0, 1)
    )
    Ly_c += -np.roll(landmask, (1, -1), axis=(0, 1)) - np.roll(
        landmask, (1, 1), axis=(0, 1)
    )

    Lx_c = np.roll(landmask, -1, axis=1) - np.roll(landmask, 1, axis=1)
    # Include x-component of diagonal neighbors
    Lx_c += np.roll(landmask, (-1, -1), axis=(1, 0)) + np.roll(
        landmask, (-1, 1), axis=(1, 0)
    )
    Lx_c += -np.roll(landmask, (1, -1), axis=(1, 0)) - np.roll(
        landmask, (1, 1), axis=(1, 0)
    )

    v_x = -Lx * (shore)
    v_y = -Ly * (shore)

    v_x_c = -Lx_c * (shore_c)
    v_y_c = -Ly_c * (shore_c)

    v_x = v_x + v_x_c
    v_y = v_y + v_y_c

    magnitude = np.sqrt(v_y**2 + v_x**2)
    # the coastal nodes between land create a problem. Magnitude there is zero
    # I force it to be 1 to avoid problems when normalising.
    ny, nx = np.where(magnitude == 0)
    magnitude[ny, nx] = 1

    v_x = v_x / magnitude
    v_y = v_y / magnitude

    return v_x, v_y

v_x, v_y = create_displacement_field(landmask.mask.astype(int))

Show code cell source

Hide code cell source

fig, ax = plt.subplots(1, 1, figsize=(7, 6), constrained_layout=True)
fig.suptitle("Figure 3. Displacement field", fontsize=16)

ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_xlim(x_incorners[0, 0], x_incorners[0, -1])
ax.set_ylim(y_incorners[0, 0], y_incorners[-1, 0])
land0 = ax.pcolormesh(
    x_centers,
    y_centers,
    lmask.mask,
    cmap="Blues_r",
)
inds = np.where(shore == 1)
coa = ax.plot(
    x_centers[inds], y_centers[inds], "o", color="y", markersize=10, markeredgecolor="k"
)
quiv = ax.quiver(
    x_centers,
    y_centers,
    v_x,
    v_y,
    color="orange",
    angles="xy",
    scale_units="xy",
    scale=25,
    width=0.005,
)
plt.show()

Step 5: Calculate the distance to the shore

In this tutorial, we will only displace particles that are within some distance (smaller than the grid size) to the shore.

For this we map the distance of the coastal nodes to the shore: Coastal nodes directly neighboring the shore are \(1dx\) away. Diagonal neighbors are \(\sqrt{2}dx\) away. The particles can then sample this field and will only be displaced when closer than a threshold value. This gives a crude estimate of the distance.

def distance_to_shore(landmask, dx=1):
    """Function that computes the distance to the shore. It is based in the
    the `get_coastal_nodes` algorithm.

    - landmask: the land mask dUilt using `make_landmask` function.
    - dx: the grid cell dimension. This is a crude approxsimation of the real
    distance (be careful).

    Output: 2D array containing the distances from shore.
    """
    ci = get_coastal_nodes(landmask)  # direct neighbors
    dist = ci * dx  # 1 dx away

    ci_d = get_coastal_nodes_diagonal(landmask)  # diagonal neighbors
    dist_d = (ci_d - ci) * np.sqrt(2 * dx**2)  # sqrt(2) dx away

    return dist + dist_d

d_2_s = distance_to_shore(landmask.mask.astype(int))

Show code cell source

Hide code cell source

fig, ax = plt.subplots(1, 1, figsize=(6, 5), constrained_layout=True)

fig.suptitle("Figure 4. Distance to shore", fontsize=16)
ax.set_ylabel("Latitude [degrees]")
ax.set_xlabel("Longitude [degrees]")

ax.pcolormesh(
    x_centers,
    y_centers,
    lmask.mask,
    cmap="Blues_r",
)
d2s = ax.scatter(x_centers, y_centers, c=d_2_s, edgecolors="k")

plt.colorbar(d2s, ax=ax, label="Distance [gridcells]")
plt.show()

Step 6: Particle and Kernels

The distance to shore, used to flag whether a particle must be displaced, is stored in a particle Variable d2s. To visualize the displacement, the zonal and meridional displacements are stored in the variables dU and dV.

To write the displacement vector to the output before displacing the particle, the SetDisplacement kernel is invoked after the advection kernel.

DisplacementParticle = parcels.Particle.add_variable(
    [
        parcels.Variable("dU"),
        parcels.Variable("dV"),
        parcels.Variable("d2s", initial=1e3),
    ]
)


def SetDisplacement(particles, fieldset):
    particles.d2s = fieldset.distance2shore[particles]
    ptcls_nearshore = particles[particles.d2s < 0.5]

    dU, dV = fieldset.dispUV[ptcls_nearshore]

    ptcls_nearshore.dx += dU * ptcls_nearshore.dt
    ptcls_nearshore.dy += dV * ptcls_nearshore.dt

Step 7: Simulation

Let us first do a simulation with the default AdvectionRK2 kernel for comparison later

And we use the following set of 9 particles

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)

runtime = timedelta(hours=100)
dt = timedelta(minutes=10)

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

kernels = parcels.kernels.AdvectionRK2

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

pset.execute(
    kernels, runtime=runtime, dt=dt, output_file=output_file, verbose_progress=False
)

INFO: Output files are stored in stuck_particles.parquet

Now let’s add the Fields we created above to the FieldSet and do a simulation to test the displacement of the particles as they approach the shore.

ds_2d = xr.Dataset(
    data_vars={
        "dispU": (["lat", "lon"], v_x),
        "dispV": (["lat", "lon"], v_y),
        "distance2shore": (["lat", "lon"], d_2_s),
    },
    coords={
        "lon": (
            ["lon"],
            ds["longitude"].values,
            {"axis": "X", "units": "degrees_east"},
        ),
        "lat": (
            ["lat"],
            ds["latitude"].values,
            {"axis": "Y", "units": "degrees_north"},
        ),
    },
)
fields = {
    "dispU": ds_2d["dispU"],
    "dispV": ds_2d["dispV"],
    "distance2shore": ds_2d["distance2shore"],
}
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)

# Combine the original fieldset with the new 2D displacement fieldset
fieldset_land = parcels.FieldSet.from_sgrid_conventions(
    ds_fset, vector_fields={"dispUV": ("dispU", "dispV")}
)
fieldset_with_land = fieldset + fieldset_land
fieldset_with_land.describe()

INFO: No depth dimension found in dataset. Assuming surface data and adding a singleton depth dimension.
| Name           | Type        |   Grid number | Interp method / value   | Parcels backend   |
|:---------------|:------------|--------------:|:------------------------|:------------------|
| U              | Field       |             0 | XLinear(...)            | WindowedArray     |
| V              | Field       |             0 | XLinear(...)            | WindowedArray     |
| UV             | VectorField |             0 | XLinear_Velocity(...)   | -                 |
| dispU          | Field       |             1 | XLinear(...)            | NumPy             |
| dispV          | Field       |             1 | XLinear(...)            | NumPy             |
| distance2shore | Field       |             1 | XLinear(...)            | NumPy             |
| dispUV         | VectorField |             1 | XLinear_Velocity(...)   | -                 |

mesh: SphericalMesh(radius=6366707.019493707)
time interval: (np.datetime64('2019-07-01T00:00:00.000000000'), np.datetime64('2019-07-08T00:00:00.000000000'))

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

kernels = [parcels.kernels.AdvectionRK2, SetDisplacement]

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

pset.execute(
    kernels, runtime=runtime, dt=dt, output_file=output_file, verbose_progress=False
)

INFO: Output files are stored in displace_particles.parquet

Step 8: Output

To visualize the effect of the displacement, the particle trajectory output can be compared to the simulation without the displacement kernel.

df_stuck = parcels.read_particlefile("stuck_particles.parquet")
df_displace = parcels.read_particlefile("displace_particles.parquet")

Show code cell source

Hide code cell source

fig, ax = plt.subplots(1, 1, figsize=(6, 5), constrained_layout=True)

fig.suptitle("Figure 5. Trajectory difference", fontsize=16)
ax.set_ylabel("Latitude [degrees]")
ax.set_xlabel("Longitude [degrees]")

ax.pcolormesh(
    x_centers,
    y_centers,
    lmask.mask,
    cmap="Blues_r",
)
for j, traj in enumerate(df_stuck.partition_by("particle_id", maintain_order=True)):
    label = "stuck" if j == 0 else None
    ax.plot(
        traj["x"],
        traj["y"],
        "-",
        linewidth=2,
        zorder=3,
        color="y",
        alpha=0.7,
        label=label,
    )
for j, traj in enumerate(df_displace.partition_by("particle_id", maintain_order=True)):
    label = "displaced" if j == 0 else None
    ax.plot(
        traj["x"],
        traj["y"],
        "-",
        linewidth=2,
        zorder=2,
        color="c",
        alpha=0.7,
        label=label,
    )
ax.plot(
    X.flatten(), Y.flatten(), "o", color="w", markersize=5, label="initial positions"
)
plt.legend()
plt.show()

Conclusion

Figure 5 shows how particles are prevented from approaching the coast in a 5 day simulation. Note that to show each computation, the integration timestep (dt) is equal to the output timestep (outputdt): 1 hour. This is relatively large, and causes the displacement to be on the order of 4 km and be relatively infrequent. It is advised to use smaller dt in real simulations.

3. Slip boundary conditions

The reason trajectories do not neatly follow the coast in A grid velocity fields is that the lack of staggering causes both velocity components to go to zero in the same way towards the cell edge. This no-slip condition can be turned into a free-slip or partial-slip condition by separately considering the cross-shore and along-shore velocity components as in the C-grids part of the Stuck particles tutorial. Each interpolation of the velocity field must then be corrected with a factor depending on the direction of the boundary.

Note

While the code below has been developed for A-grids, we have experience that it can also work well with B-grids. So if you encounter stuck particles in a B-grid model then it might be worth to use the parcels.interpolators.XPartialslip() or parcels.interpolators.XFreeslip() described below.

These boundary conditions have been implemented in Parcels as parcels.interpolators.XPartialslip() and parcels.interpolators.XFreeslip(), which we will show in the plot below

Show code cell source

Hide code cell source

cells_x = np.array([[0, 0], [1, 1], [2, 2]])
cells_y = np.array([[0, 1], [0, 1], [0, 1]])
U0 = 1
V0 = 1
U = np.array([U0, U0, 0, 0, 0, 0])
V = np.array([V0, V0, 0, 0, 0, 0])
xsi = np.linspace(0.001, 0.999)

u_interp = U0 * (1 - xsi)
v_interp = V0 * (1 - xsi)

u_freeslip = u_interp
v_freeslip = v_interp / (1 - xsi)

u_partslip = u_interp
v_partslip = v_interp * (1 - 0.5 * xsi) / (1 - xsi)


fig = plt.figure(figsize=(15, 4), constrained_layout=True)
fig.suptitle("Figure 6. Boundary conditions", fontsize=16)
gs = gridspec.GridSpec(ncols=3, nrows=1, figure=fig)

ax0 = fig.add_subplot(gs[0, 0])

ax0.pcolormesh(cells_x, cells_y, np.array([[0], [1]]), cmap="Greys", edgecolor="k")
ax0.scatter(cells_x, cells_y, c="w", edgecolor="k")
ax0.quiver(cells_x, cells_y, U, V, scale=15)

ax0.plot(xsi, u_interp, linewidth=5, label="u_interpolation")
ax0.plot(xsi, v_interp, linestyle="dashed", linewidth=5, label="v_interpolation")
ax0.set_xlim(-0.3, 2.3)
ax0.set_ylim(-0.5, 1.5)
ax0.set_ylabel("u - v [-]", fontsize=14)
ax0.set_xlabel(r"$\xi$", fontsize=14)
ax0.set_title("A) Bilinear interpolation")
ax0.legend(loc="lower right")


ax1 = fig.add_subplot(gs[0, 1])

ax1.pcolormesh(cells_x, cells_y, np.array([[0], [1]]), cmap="Greys", edgecolor="k")
ax1.scatter(cells_x, cells_y, c="w", edgecolor="k")
ax1.quiver(cells_x, cells_y, U, V, scale=15)

ax1.plot(xsi, u_freeslip, linewidth=5, label="u_freeslip")
ax1.plot(xsi, v_freeslip, linestyle="dashed", linewidth=5, label="v_freeslip")
ax1.set_xlim(-0.3, 2.3)
ax1.set_ylim(-0.5, 1.5)
ax1.set_xlabel(r"$\xi$", fontsize=14)
ax1.text(0.0, 1.3, r"$v_{freeslip} = v_{interpolation}*\frac{1}{1-\xi}$", fontsize=18)
ax1.set_title("B) Free slip condition")
ax1.legend(loc="lower right")

ax2 = fig.add_subplot(gs[0, 2])

ax2.pcolormesh(cells_x, cells_y, np.array([[0], [1]]), cmap="Greys", edgecolor="k")
ax2.scatter(cells_x, cells_y, c="w", edgecolor="k")
ax2.quiver(cells_x, cells_y, U, V, scale=15)

ax2.plot(xsi, u_partslip, linewidth=5, label="u_partialslip")
ax2.plot(xsi, v_partslip, linestyle="dashed", linewidth=5, label="v_partialslip")
ax2.set_xlim(-0.3, 2.3)
ax2.set_ylim(-0.5, 1.5)
ax2.set_xlabel(r"$\xi$", fontsize=14)
ax2.text(
    0.0,
    1.3,
    r"$v_{partialslip} = v_{interpolation}*\frac{1-1/2\xi}{1-\xi}$",
    fontsize=18,
)
ax2.set_title("C) Partial slip condition")
ax2.legend(loc="lower right");

Consider a grid cell with a solid boundary to the right and vectors \((U0, V0)\) = \((1, 1)\) on the left-hand nodes, as in figure 6. Parcels bilinear interpolation will interpolate in the \(x\) and \(y\) directions. This cell is invariant in the \(y\)-direction, we will only consider the effect in the direction normal to the boundary. In the x-direction, both u and v will be interpolated along \(\xi\), the normalised \(x\)-coordinate within the cell. This is plotted with the blue and orange dashed lines in figure 6A.

A free slip boundary condition is defined with \(\frac{\delta v}{\delta \xi}=0\). This means that the tangential velocity is constant in the direction normal to the boundary. This can be achieved in a kernel after interpolation by dividing by \((1-\xi)\). The resulting velocity profiles are shown in figure 6B.

A partial slip boundary condition is defined with a tangential velocity profile that decreases toward the boundary, but not to zero. This can be achieved by multiplying the interpolated velocity by \(\frac{1-1/2\xi}{1-\xi}\). This is shown in figure 6C.

For each direction and boundary condition a different factor must be used (where \(\xi\) and \(\eta\) are the normalised x- and y-coordinates within the cell, respectively):

• Free slip

  1: \(f_u = \frac{1}{\eta}\)

  2: \(f_u = \frac{1}{(1-\eta)}\)

  4: \(f_v = \frac{1}{\xi}\)

  8: \(f_v = \frac{1}{(1-\xi)}\)

• Partial slip

  1: \(f_u = \frac{1/2+1/2\eta}{\eta}\)

  2: \(f_u = \frac{1-1/2\eta}{1-\eta}\)

  4: \(f_v = \frac{1/2+1/2\xi}{\xi}\)

  8: \(f_v = \frac{1-1/2\xi}{1-\xi}\)

We now simulate the three different boundary conditions by advecting the 9 particles from above in a time-evolving dataset from the Copernicus Marine Data Store.

First up is the partialslip interpolation

fieldset.UV.interp_method = parcels.interpolators.XPartialslip()

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

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

pset.execute(
    parcels.kernels.AdvectionRK2,
    runtime=runtime,
    dt=dt,
    output_file=output_file,
    verbose_progress=False,
)

INFO: Output files are stored in partialslip_particles.parquet

And then we also use the freeslip interpolation

fieldset.UV.interp_method = parcels.interpolators.XFreeslip()

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

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

pset.execute(
    parcels.kernels.AdvectionRK2,
    runtime=runtime,
    dt=dt,
    output_file=output_file,
    verbose_progress=False,
)

INFO: Output files are stored in freeslip_particles.parquet

Now we can load and plot the two interpolation_methods

Show code cell source

Hide code cell source

df_freeslip = parcels.read_particlefile("freeslip_particles.parquet")
df_partialslip = parcels.read_particlefile("partialslip_particles.parquet")

fig, ax = plt.subplots(1, 1, figsize=(6, 5), constrained_layout=True)

fig.suptitle("Figure 7. Solution comparison", fontsize=16)
ax.set_ylabel("Latitude [degrees]")
ax.set_xlabel("Longitude [degrees]")

ax.pcolormesh(
    x_centers,
    y_centers,
    lmask.mask,
    cmap="Blues_r",
)
for j, traj in enumerate(df_stuck.partition_by("particle_id", maintain_order=True)):
    label = "stuck" if j == 0 else None
    ax.plot(
        traj["x"],
        traj["y"],
        "-",
        linewidth=2,
        zorder=3,
        color="y",
        alpha=0.7,
        label=label,
    )
for j, traj in enumerate(
    df_partialslip.partition_by("particle_id", maintain_order=True)
):
    label = "partialslip" if j == 0 else None
    ax.plot(
        traj["x"],
        traj["y"],
        "-",
        linewidth=2,
        zorder=2,
        color="r",
        alpha=0.7,
        label=label,
    )
for j, traj in enumerate(df_freeslip.partition_by("particle_id", maintain_order=True)):
    label = "freelslip" if j == 0 else None
    ax.plot(
        traj["x"],
        traj["y"],
        "-",
        linewidth=2,
        zorder=1,
        color="m",
        alpha=0.7,
        label=label,
    )
ax.plot(
    X.flatten(), Y.flatten(), "o", color="w", markersize=5, label="initial positions"
)
plt.legend()
plt.show()

Figure 7 shows the influence of the different solutions on the particle trajectories near the shore.