Tutorial on implementing boundary conditions in an A grid

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 all 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

This notebook is mainly focused on the second part, comparing the different modifications to the trajectory. The first part 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.

In [1]:
import numpy as np
import numpy.ma as ma
from netCDF4 import Dataset
import xarray as xr
from scipy import interpolate

from parcels import FieldSet, ParticleSet, JITParticle, ScipyParticle, AdvectionRK4, Variable, Field,GeographicPolar,Geographic
from datetime import timedelta as delta

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.colors import ListedColormap
from matplotlib.lines import Line2D
from copy import copy
import cmocean
INFO: Compiled ParcelsRandom ==> C:\Users\GEBRUI~1\AppData\Local\Temp\parcels-tmp\libparcels_random_041e4d11-367a-41de-bbaa-c9f5076b6e99.so

1. Particle deletion

The simplest way to avoid trajectories that interact with the coastline is to remove them entirely. To do this, all Particle objects have a delete function that can be invoked in a kernel using particle.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.

Import a velocity field - the A gridded SMOC product

In [2]:
file_path = "SMOC_20190704_R20190705.nc"
model = xr.open_dataset(file_path)
In [3]:
# --------- Define meshgrid coordinates to plot velocity field with matplotlib pcolormesh ---------
latmin = 1595
latmax = 1612
lonmin = 2235
lonmax = 2260

# Velocity nodes
lon_vals, lat_vals = np.meshgrid(model['longitude'], model['latitude'])
lons_plot = lon_vals[latmin:latmax,lonmin:lonmax]
lats_plot = lat_vals[latmin:latmax,lonmin:lonmax]

dlon = 1/12
dlat = 1/12

# Centers of the gridcells formed by 4 nodes = velocity nodes + 0.5 dx
x = model['longitude'][:-1]+np.diff(model['longitude'])/2
y = model['latitude'][:-1]+np.diff(model['latitude'])/2
lon_centers, lat_centers = np.meshgrid(x, y)

color_land = copy(plt.get_cmap('Reds'))(0)
color_ocean = copy(plt.get_cmap('Reds'))(128)

Make a landmask where land = 1 and ocean = 0.

In [4]:
def make_landmask(fielddata):
    """Returns landmask where land = 1 and ocean = 0
    fielddata is a netcdf file.
    """
    datafile = Dataset(fielddata)

    landmask = datafile.variables['uo'][0, 0]
    landmask = np.ma.masked_invalid(landmask)
    landmask = landmask.mask.astype('int')

    return landmask
In [5]:
landmask = make_landmask(file_path)
In [6]:
# Interpolate the landmask to the cell centers - only cells with 4 neighbouring land points will be land
fl = interpolate.interp2d(model['longitude'],model['latitude'],landmask)

l_centers = fl(lon_centers[0,:],lat_centers[:,0])  

lmask = np.ma.masked_values(l_centers,1) # land when interpolated value == 1
In [7]:
fig = plt.figure(figsize=(12,5))
fig.suptitle('Figure 1. Landmask', fontsize=18, y=1.01)
gs = gridspec.GridSpec(ncols=2, nrows=1, figure=fig)

ax0 = fig.add_subplot(gs[0, 0])
ax0.set_title('A) lazy use of pcolormesh', fontsize=11)
ax0.set_ylabel('Latitude [degrees]')
ax0.set_xlabel('Longitude [degrees]')

land0 = ax0.pcolormesh(lons_plot, lats_plot, landmask[latmin:latmax,lonmin:lonmax],cmap='Reds_r', shading='auto')
ax0.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=20,cmap='Reds_r',vmin=-0.05,vmax=0.05,edgecolors='k')

custom_lines = [Line2D([0], [0], c = color_ocean, marker='o', markersize=10, markeredgecolor='k', lw=0),
                Line2D([0], [0], c = color_land, marker='o', markersize=10, markeredgecolor='k', lw=0)]
ax0.legend(custom_lines, ['ocean point', 'land point'], bbox_to_anchor=(.01,.93), loc='center left', borderaxespad=0.,framealpha=1)

ax1 = fig.add_subplot(gs[0, 1])
ax1.set_title('B) correct A grid representation in Parcels', fontsize=11)
ax1.set_ylabel('Latitude [degrees]')
ax1.set_xlabel('Longitude [degrees]')

land1 = ax1.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
ax1.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=20,cmap='Reds_r',vmin=-0.05,vmax=0.05,edgecolors='k')

ax1.legend(custom_lines, ['ocean point', 'land point'], bbox_to_anchor=(.01,.93), loc='center left', borderaxespad=0.,framealpha=1)
Out[7]:
<matplotlib.legend.Legend at 0x1fb804682e8>

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.

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].

In [8]:
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 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 -= 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 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 -= 4*landmask
    shore = np.ma.masked_array(landmask, mask_lap < 0)
    shore = shore.mask.astype('int')

    return shore
In [9]:
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
In [10]:
coastal = get_coastal_nodes_diagonal(landmask)
shore = get_shore_nodes_diagonal(landmask)
In [11]:
fig = plt.figure(figsize=(10,4), constrained_layout=True)
fig.suptitle('Figure 2. Coast and Shore', fontsize=18, y=1.04)
gs = gridspec.GridSpec(ncols=2, nrows=1, figure=fig)


ax0 = fig.add_subplot(gs[0, 0])
land0 = ax0.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
coa = ax0.scatter(lons_plot,lats_plot, c=coastal[latmin:latmax,lonmin:lonmax], cmap='Reds_r', s=50)
ax0.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=20,cmap='Reds_r',vmin=-0.05,vmax=0.05)

ax0.set_title('Coast')
ax0.set_ylabel('Latitude [degrees]')
ax0.set_xlabel('Longitude [degrees]')

custom_lines = [Line2D([0], [0], c = color_ocean, marker='o', markersize=5, lw=0),
                Line2D([0], [0], c = color_ocean, marker='o', markersize=7, markeredgecolor='w', markeredgewidth=2, lw=0),
                Line2D([0], [0], c = color_land, marker='o', markersize=7, markeredgecolor='firebrick', lw=0)]
ax0.legend(custom_lines, ['ocean node', 'coast node', 'land node'], bbox_to_anchor=(.01,.9), loc='center left', borderaxespad=0.,framealpha=1, facecolor='silver')


ax1 = fig.add_subplot(gs[0, 1])
land1 = ax1.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
sho = ax1.scatter(lons_plot,lats_plot, c=shore[latmin:latmax,lonmin:lonmax], cmap='Reds_r', s=50)
ax1.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=20,cmap='Reds_r',vmin=-0.05,vmax=0.05)

ax1.set_title('Shore')
ax1.set_ylabel('Latitude [degrees]')
ax1.set_xlabel('Longitude [degrees]')

custom_lines = [Line2D([0], [0], c = color_ocean, marker='o', markersize=5, lw=0),
                Line2D([0], [0], c = color_land, marker='o', markersize=7, markeredgecolor='w', markeredgewidth=2, lw=0),
                Line2D([0], [0], c = color_land, marker='o', markersize=7, markeredgecolor='firebrick', lw=0)]
ax1.legend(custom_lines, ['ocean node', 'shore node', 'land node'], bbox_to_anchor=(.01,.9), loc='center left', borderaxespad=0.,framealpha=1, facecolor='silver')
Out[11]:
<matplotlib.legend.Legend at 0x1fb8051eda0>

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 normalize them by their magnitude.

In [12]:
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 dUilt using `make_landmask`.
    - double_cell: Boolean for determining if you want a double cell.
      Default set to False.

    Output: two 2D arrays, one for each camponent of the velocity.
    """
    shore = get_shore_nodes(landmask)
    shore_d = get_shore_nodes_diagonal(landmask) # bordering ocean directly and diagonally
    shore_c = shore_d - shore                    # corner nodes that only border ocean diagonally
    
    Ly = np.roll(landmask, -1, axis=0) - np.roll(landmask, 1, axis=0) # Simple derivative
    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)
    Ly_c += np.roll(landmask, (-1,-1), axis=(0,1)) + np.roll(landmask, (-1,1), axis=(0,1)) # Include y-component of diagonal neighbours
    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)
    Lx_c += np.roll(landmask, (-1,-1), axis=(1,0)) + np.roll(landmask, (-1,1), axis=(1,0)) # Include x-component of diagonal neighbours
    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 normalizing.
    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
In [13]:
v_x, v_y = create_displacement_field(landmask)
In [14]:
fig = plt.figure(figsize=(7,6), constrained_layout=True)
fig.suptitle('Figure 3. Displacement field', fontsize=18, y=1.04)
gs = gridspec.GridSpec(ncols=1, nrows=1, figure=fig)

ax0 = fig.add_subplot(gs[0, 0])
land = ax0.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
ax0.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=30,cmap='Reds_r',vmin=-0.05,vmax=0.05, edgecolors='k')
quiv = ax0.quiver(lons_plot,lats_plot,v_x[latmin:latmax,lonmin:lonmax],v_y[latmin:latmax,lonmin:lonmax],color='orange',angles='xy', scale_units='xy', scale=19, width=0.005)

ax0.set_ylabel('Latitude [degrees]')
ax0.set_xlabel('Longitude [degrees]')

custom_lines = [Line2D([0], [0], c = color_ocean, marker='o', markersize=10, markeredgecolor='k', lw=0),
                Line2D([0], [0], c = color_land, marker='o', markersize=10, markeredgecolor='k', lw=0)]
ax0.legend(custom_lines, ['ocean point', 'land point'], bbox_to_anchor=(.01,.93), loc='center left', borderaxespad=0.,framealpha=1)
Out[14]:
<matplotlib.legend.Legend at 0x1fb804fd860>

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.

In [15]:
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 neighbours
    dist = ci*dx                     # 1 dx away
    
    ci_d = get_coastal_nodes_diagonal(landmask) # diagonal neighbours
    dist_d = (ci_d - ci)*np.sqrt(2*dx**2)       # sqrt(2) dx away
        
    return dist+dist_d
In [16]:
d_2_s = distance_to_shore(landmask)
In [17]:
fig = plt.figure(figsize=(6,5), constrained_layout=True)

ax0 = fig.add_subplot()
ax0.set_title('Figure 4. Distance to shore', fontsize=18)
ax0.set_ylabel('Latitude [degrees]')
ax0.set_xlabel('Longitude [degrees]')

land = ax0.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
d2s = ax0.scatter(lons_plot,lats_plot, c=d_2_s[latmin:latmax,lonmin:lonmax])

plt.colorbar(d2s,ax=ax0, label='Distance [gridcells]')
Out[17]:
<matplotlib.colorbar.Colorbar at 0x1fb80ca6f60>

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 set_displacement kernel is invoked after the advection kernel. Then only in the next timestep are particles displaced by displace, before resuming the advection.

In [18]:
class DisplacementParticle(JITParticle):
    dU = Variable('dU')
    dV = Variable('dV')
    d2s = Variable('d2s', initial=1e3)
    
def set_displacement(particle, fieldset, time):
    particle.d2s = fieldset.distance2shore[time, particle.depth,
                               particle.lat, particle.lon]
    if  particle.d2s < 0.5:
        dispUab = fieldset.dispU[time, particle.depth, particle.lat,
                               particle.lon]
        dispVab = fieldset.dispV[time, particle.depth, particle.lat,
                               particle.lon]
        particle.dU = dispUab
        particle.dV = dispVab
    else:
        particle.dU = 0.
        particle.dV = 0.
    
def displace(particle, fieldset, time):    
    if  particle.d2s < 0.5:
        particle.lon += particle.dU*particle.dt
        particle.lat += particle.dV*particle.dt

Simulation

Now let us test the displacement of the particles as they approach the shore.

In [19]:
SMOCfile = 'SMOC_201907*.nc'
filenames = {'U': SMOCfile,
             'V': SMOCfile}

variables = {'U': 'uo',
             'V': 'vo'}

dimensions = {'U': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'},
              'V': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'}}

fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
In [20]:
u_displacement = v_x
v_displacement = v_y
fieldset.add_field(Field('dispU', data=u_displacement,
                         lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat,
                         mesh='spherical'))
fieldset.add_field(Field('dispV', data=v_displacement,
                         lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat,
                         mesh='spherical'))
fieldset.dispU.units = GeographicPolar()
fieldset.dispV.units = Geographic()
In [21]:
fieldset.add_field(Field('landmask', landmask,
                         lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat,
                         mesh='spherical'))
fieldset.add_field(Field('distance2shore', d_2_s,
                         lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat,
                         mesh='spherical'))
In [22]:
npart = 3  # number of particles to be released
lon = np.linspace(7, 7.2, npart, dtype=np.float32)
lat = np.linspace(53.45, 53.6, npart, dtype=np.float32)
lons, lats = np.meshgrid(lon,lat)
time = np.zeros(lons.size)

pset = ParticleSet(fieldset=fieldset, pclass=DisplacementParticle, lon=lons, lat=lats, time=time)

kernels = pset.Kernel(displace)+pset.Kernel(AdvectionRK4)+pset.Kernel(set_displacement)

output_file = pset.ParticleFile(name="SMOC-disp.nc", outputdt=delta(hours=1))

pset.execute(kernels, runtime=delta(hours=119), dt=delta(hours=1),
             output_file=output_file)
output_file.close()
INFO: Compiled ArrayDisplacementParticledisplaceAdvectionRK4set_displacement ==> C:\Users\GEBRUI~1\AppData\Local\Temp\parcels-tmp\libbc6871644b3922eabb2953053ca490e6_0.dll
INFO: Temporary output files are stored in out-QCOSPTBA.
INFO: You can use "parcels_convert_npydir_to_netcdf out-QCOSPTBA" to convert these to a NetCDF file during the run.
100% (428400.0 of 428400.0) |############| Elapsed Time: 0:01:30 Time:  0:01:30

Output

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

In [23]:
ds_SMOC = xr.open_dataset('SMOC.nc')
ds_SMOC_disp = xr.open_dataset('SMOC-disp.nc')
In [24]:
fig = plt.figure(figsize=(16,4), facecolor='silver', constrained_layout=True)
fig.suptitle('Figure 5. Trajectory difference', fontsize=18, y=1.06)
gs = gridspec.GridSpec(ncols=4, nrows=1, width_ratios=[1,1,1,0.3], figure=fig)

ax0 = fig.add_subplot(gs[0, 0])
ax0.set_ylabel('Latitude [degrees]')
ax0.set_xlabel('Longitude [degrees]')
ax0.set_title('A) No displacement', fontsize=14, fontweight = 'bold')
ax0.set_xlim(6.9, 7.6)
ax0.set_ylim(53.4, 53.8)

land = ax0.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
ax0.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=50,cmap='Reds_r',vmin=-0.05,vmax=0.05, edgecolors='k')
ax0.plot(ds_SMOC['lon'].T, ds_SMOC['lat'].T,linewidth=3, zorder=1)
ax0.scatter(ds_SMOC['lon'], ds_SMOC['lat'], color='limegreen', zorder=2)

n_p0 = 0
ax1 = fig.add_subplot(gs[0, 1])
ax1.set_ylabel('Latitude [degrees]')
ax1.set_xlabel('Longitude [degrees]')
ax1.set_title('B) Displacement trajectory '+str(n_p0), fontsize=14, fontweight = 'bold')
ax1.set_xlim(6.9, 7.3)
ax1.set_ylim(53.4, 53.55)

land = ax1.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
ax1.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=50,cmap='Reds_r',vmin=-0.05,vmax=0.05, edgecolors='k')
quiv = ax1.quiver(lons_plot,lats_plot,v_x[latmin:latmax,lonmin:lonmax],v_y[latmin:latmax,lonmin:lonmax],color='orange', scale=19, width=0.005)
ax1.plot(ds_SMOC_disp['lon'][n_p0].T, ds_SMOC_disp['lat'][n_p0].T,linewidth=3, zorder=1)
ax1.scatter(ds_SMOC['lon'][n_p0], ds_SMOC['lat'][n_p0], color='limegreen', zorder=2)
ax1.scatter(ds_SMOC_disp['lon'][n_p0], ds_SMOC_disp['lat'][n_p0], cmap='viridis_r', zorder=2)
ax1.quiver(ds_SMOC_disp['lon'][n_p0], ds_SMOC_disp['lat'][n_p0],ds_SMOC_disp['dU'][n_p0], ds_SMOC_disp['dV'][n_p0], color='w',angles='xy', scale_units='xy', scale=2e-4, zorder=3)

n_p1 = 4
ax2 = fig.add_subplot(gs[0, 2])
ax2.set_ylabel('Latitude [degrees]')
ax2.set_xlabel('Longitude [degrees]')
ax2.set_title('C) Displacement trajectory '+str(n_p1), fontsize=14, fontweight = 'bold')
ax2.set_xlim(7., 7.6)
ax2.set_ylim(53.4, 53.8)

land = ax2.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
ax2.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=50,cmap='Reds_r',vmin=-0.05,vmax=0.05, edgecolors='k')
q1 = ax2.quiver(lons_plot,lats_plot,v_x[latmin:latmax,lonmin:lonmax],v_y[latmin:latmax,lonmin:lonmax],color='orange', scale=19, width=0.005)
ax2.plot(ds_SMOC_disp['lon'][n_p1].T, ds_SMOC_disp['lat'][n_p1].T,linewidth=3, zorder=1)
ax2.scatter(ds_SMOC['lon'][n_p1], ds_SMOC['lat'][n_p1], color='limegreen', zorder=2)
ax2.scatter(ds_SMOC_disp['lon'][n_p1], ds_SMOC_disp['lat'][n_p1], cmap='viridis_r', zorder=2)
q2 = ax2.quiver(ds_SMOC_disp['lon'][n_p1], ds_SMOC_disp['lat'][n_p1],ds_SMOC_disp['dU'][n_p1], ds_SMOC_disp['dV'][n_p1], color='w',angles='xy', scale_units='xy', scale=2e-4, zorder=3)


ax3 = fig.add_subplot(gs[0, 3])
ax3.axis('off')
custom_lines = [Line2D([0], [0], c = 'tab:blue', marker='o', markersize=10),
                Line2D([0], [0], c = 'limegreen', marker='o', markersize=10),
                Line2D([0], [0], c = color_ocean, marker='o', markersize=10, markeredgecolor='k', lw=0),
                Line2D([0], [0], c = color_land, marker='o', markersize=10, markeredgecolor='k', lw=0)]
ax3.legend(custom_lines, ['with displacement', 'without displacement', 'ocean point', 'land point'], bbox_to_anchor=(0.,0.6), loc='center left', borderaxespad=0.,framealpha=1)

ax2.quiverkey(q1, 1.3, 0.9, 2, 'displacement field', coordinates='axes')
ax2.quiverkey(q2, 1.3, 0.8, 1e-5, 'particle displacement', coordinates='axes')
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.

In [25]:
d2s_cmap = copy(plt.get_cmap('cmo.deep_r'))
d2s_cmap.set_over('gold')
In [26]:
fig = plt.figure(figsize=(11,6), constrained_layout=True)

ax0 = fig.add_subplot()
ax0.set_title('Figure 6. Distance to shore', fontsize=18)
land = ax0.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')
ax0.scatter(lons_plot, lats_plot, c=landmask[latmin:latmax,lonmin:lonmax],s=50,cmap='Reds_r', edgecolor='k',vmin=-0.05,vmax=0.05)
ax0.plot(ds_SMOC_disp['lon'].T, ds_SMOC_disp['lat'].T,linewidth=3, zorder=1)
d2s = ax0.scatter(ds_SMOC_disp['lon'], ds_SMOC_disp['lat'], c=ds_SMOC_disp['d2s'],cmap=d2s_cmap, s=20,vmax=0.5, zorder=2)
q2 = ax0.quiver(ds_SMOC_disp['lon'], ds_SMOC_disp['lat'],ds_SMOC_disp['dU'], ds_SMOC_disp['dV'], color='k',angles='xy', scale_units='xy', scale=2.3e-4, width=0.003, zorder=3)

ax0.set_xlim(6.9, 8)
ax0.set_ylim(53.4, 53.8)
ax0.set_ylabel('Latitude [degrees]')
ax0.set_xlabel('Longitude [degrees]')
plt.colorbar(d2s,ax=ax0, label='Distance [gridcells]',extend='max')

color_land = copy(plt.get_cmap('Reds'))(0)
color_ocean = copy(plt.get_cmap('Reds'))(128)

custom_lines = [Line2D([0], [0], c = color_ocean, marker='o', markersize=10, markeredgecolor='k', lw=0),
                Line2D([0], [0], c = color_land, marker='o', markersize=10, markeredgecolor='k', lw=0)]
ax0.legend(custom_lines, ['ocean point', 'land point'], bbox_to_anchor=(.01,.95), loc='center left', borderaxespad=0.,framealpha=1)
Out[26]:
<matplotlib.legend.Legend at 0x1fb80d801d0>

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 a staggered C-grid. Each interpolation of the velocity field must then be corrected with a factor depending on the direction of the boundary. To do this, for each cell the direction to shore is mapped and the velocity in the advection kernel is modified according to this direction and the relative position in the cell.

Map boundary directions

In [27]:
def get_boundary_direction(landmask):
    """Function that detects the boundaries for a gridcell [i,j] defined by the nodes [i,j], [i+1,j], [i+1,j+1] and [i,j+1]
    When i and j increase with longitude and latitude, this defines the gridcell [i,j] with node [i,j] in the southwest corner.
    Directions are mapped according to figure 7.
    """
    bound_dir = np.logical_and(landmask, np.roll(landmask, -1, axis=1)).astype(int)
    bound_dir += 2*np.logical_and(np.roll(landmask, -1, axis=0), np.roll(landmask, (-1,-1), axis=(0,1))).astype(int)
    bound_dir += 4*np.logical_and(landmask, np.roll(landmask, -1, axis=0)).astype(int)
    bound_dir += 8*np.logical_and(np.roll(landmask, -1, axis=1), np.roll(landmask, (-1,-1), axis=(0,1))).astype(int)
    bound_dir = bound_dir[:-1,:-1]
    
    return bound_dir
In [28]:
bound_dir = get_boundary_direction(landmask)
In [29]:
# Colormap
dir_cmap = copy(plt.get_cmap('jet',11))
dir_cmap.set_under([1,1,1,1])
dir_cmap.set_over([0,0,0,1])
In [30]:
fig = plt.figure(figsize=(11,5), constrained_layout=True)
fig.suptitle('Figure 7. Mapping boundary directions', fontsize=18, y=1.06)
gs = gridspec.GridSpec(ncols=2, nrows=1, figure=fig)

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

land1 = ax0.pcolormesh(lon_vals[latmin:latmax+1,lonmin:lonmax+1], lat_vals[latmin:latmax+1,lonmin:lonmax+1], lmask.mask[latmin:latmax,lonmin:lonmax],cmap='Reds_r')

lm = ax0.scatter(lons_plot, lats_plot,c=landmask[latmin:latmax,lonmin:lonmax],s=30,cmap='Reds_r', vmax=0.5,edgecolor='k')
bd = ax0.scatter(lon_centers[latmin:latmax,lonmin:lonmax], lat_centers[latmin:latmax,lonmin:lonmax],c=bound_dir[latmin:latmax,lonmin:lonmax],s=80,cmap=dir_cmap, vmin=0.5,vmax=11.5)

ax0.set_title('Map')

plt.colorbar(bd, label = 'Direction', extend='both')
ax0.set_xlim(6.5,7.5)
ax0.set_ylim(53,54)
ax0.set_ylabel('Latitude [degrees]')
ax0.set_xlabel('Longitude [degrees]')


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

ax1.pcolormesh(np.array([[0,0],[1,1]]),np.array([[0,1],[0,1]]),np.array([[1]]), cmap='Reds',edgecolor='k')
ax1.scatter(np.array([[0,0],[1,1]]),np.array([[0,1],[0,1]]), c='w', edgecolor='k')
ax1.quiver(np.array([0,0.5,1,0.5]),np.array([0.5,0,0.5,1]),np.array([-1,0,1,0]),np.array([0,-1,0,1]))
ax1.text(-0.08,0.55,'4',fontsize=14)
ax1.text(0.55,-0.08,'1',fontsize=14)
ax1.text(1.05,0.55,'8',fontsize=14)
ax1.text(0.55,1.05,'2',fontsize=14)
ax1.text(-0.1,-0.1,'[i,j]',fontsize=11)
ax1.text(0.45,0.5,'[i,j]',fontsize=11)
ax1.text(-0.1,1.05,'[i,j+1]',fontsize=11)
ax1.text(0.95,-0.1,'[i+1,j]',fontsize=11)
ax1.text(0.95,1.05,'[i+1,j+1]',fontsize=11)
ax1.set_xlim(-0.3,1.3)
ax1.set_ylim(-0.3,1.3)
ax1.set_title('Binary values for boundary directions')
Out[30]:
Text(0.5, 1.0, 'Binary values for boundary directions')

Modify velocities

There are different ways to modify the velocity to prevent particles from reaching the land. In the previous example we have added a displacement field. Here we modify the velocity in the advection kernel, after Parcels' bilinear interpolation. If you would like to contribute to Parcels, feel free to try to implement an interpolation scheme in Parcels that includes the current corrections.

Here two kernels are presented. The first kernel modifies the tangential velocity as if a free-slip condition were implemented. The second kernel modifies the tangential velocity to represent a partial-slip boundary condition.

In [31]:
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,1)

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-.5*xsi)/(1-xsi)
C:\Users\Gebruiker\anaconda3\envs\py3_parcels\lib\site-packages\ipykernel_launcher.py:14: RuntimeWarning: invalid value encountered in true_divide
  
C:\Users\Gebruiker\anaconda3\envs\py3_parcels\lib\site-packages\ipykernel_launcher.py:17: RuntimeWarning: invalid value encountered in true_divide
In [32]:
fig = plt.figure(figsize=(15,4), constrained_layout=True)
fig.suptitle('Figure 8. Boundary conditions', fontsize=18, y=1.06)
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., 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., 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')
Out[32]:
<matplotlib.legend.Legend at 0x1fb811395c0>