# 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:

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)


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

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


In [7]:
fig = plt.figure(figsize=(12,5))
gs = gridspec.GridSpec(ncols=2, nrows=1, figure=fig)

ax0.set_title('A) lazy use of pcolormesh', fontsize=11)
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)

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

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

return coastal

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

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

return coastal

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

return shore

In [10]:
coastal = get_coastal_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)

coa = ax0.scatter(lons_plot,lats_plot, c=coastal[latmin:latmax,lonmin:lonmax], cmap='Reds_r', s=50)

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')

sho = ax1.scatter(lons_plot,lats_plot, c=shore[latmin:latmax,lonmin:lonmax], cmap='Reds_r', s=50)

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

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

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

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)

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.set_title('Figure 4. Distance to shore', fontsize=18)
ax0.set_ylabel('Latitude [degrees]')
ax0.set_xlabel('Longitude [degrees]')

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
lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat,
mesh='spherical'))
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'))
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)

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

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

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

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.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.set_title('Figure 6. Distance to shore', fontsize=18)
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 = 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)

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