partial slip and diffusion

[cgrimaldi]

Question

Question

Hi
I am doing standard particle tracking using the partial slip option to avoid particles ending up on land, but because I am also using uniform diffusion, the particles still end up on land. Looking for recommendations on how to avoid this happening?

I was wondering if this option could be viable: identify the coastal nodes and turn off the diffusion for the particles in there?

Supporting code/error messages

# Paste your code within this block

#!/usr/bin/env python

coding: utf-8

Import libraries

import random
import numpy as np
import glob
import numpy.ma as ma
import xarray as xr
import scipy
import zarr
from shapely.geometry import Point
from datetime import timedelta as delta
from datetime import datetime as datetime

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
import matplotlib.gridspec as gridspec
from matplotlib.colors import ListedColormap
from matplotlib.lines import Line2D
from copy import copy
import netCDF4 as nc
import pandas
import geopandas as gpd
import os

from parcels import StatusCode,FieldSet,Variable,ParticleSet,JITParticle,ParcelsRandom,ScipyParticle, AdvectionRK4, Variable,DiffusionUniformKh, Field,GeographicPolar,Geographic

Define kernels

def CheckOutOfBounds(particle, fieldset, time):
if particle.state == StatusCode.ErrorOutOfBounds:
particle.delete()

Make the pset: release from shapefiles

coord = gpd.read_file('mng_ex_ph.gpkg')

make the pset

num_particles=100
lats, lons = [], []
for ii in range(num_particles):
coord['Point' + str(ii)] = coord['geometry'].apply(random_point_in_poly)
lats = [*lats,*list(coord['Point'+str(ii)].y)]
lons = [*lons,*list(coord['Point'+str(ii)].x)]

List and open the ncfiles

path='/export/scratch/'
#set working dir
os.chdir(path)
arr =glob.glob('*.nc')

for i in range(0,len(arr)):
ds = xr.open_dataset(arr[i], decode_times=False)
# Make the fieldset
variables = {'U': 'u_eastward', 'V': 'v_northward'}
dimensions = {'U': { 'time': 'ocean_time','depth':'s_rho','lat': 'lat_rho','lon': 'lon_rho'},
'V': { 'time': 'ocean_time','depth':'s_rho','lat': 'lat_rho','lon': 'lon_rho'}}
fieldset = FieldSet.from_xarray_dataset(ds, variables, dimensions, mesh='spherical',allow_time_extrapolation = 'True', interp_method={
"U": "partialslip",
"V": "partialslip"})
fieldset.add_constant_field("Kh_zonal",10, mesh='spherical') # spherical in degree
fieldset.add_constant_field("Kh_meridional",10, mesh='spherical')

#Make the landmask
landmask = ds.variables['u_eastward'][0, 0,:,:]
landmask = (landmask == 0).astype('int')
fieldset.add_field(Field(name="landmask", data=landmask, lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat, mesh="spherical"))
fieldset.add_constant('lon_max',125.1)
fieldset.add_constant('lon_min',107.7)
fieldset.add_constant('lat_max',-14.6)
fieldset.add_constant('lat_min',-27.1)

depths = np.full(len(lons), -0.01666667)
repeatdt = delta(days=3) 

# Run the ptracking - One time release
pathout='/export/scratch/'

# Advection Diff
pset = ParticleSet(fieldset=fieldset, pclass=JITParticle, lon=lons, lat=lats, depth=depths,repeatdt=repeatdt)

output_file= pset.ParticleFile(name=pathout+arr[i],outputdt=delta(days=3))

kernels = [AdvectionRK4,DiffusionUniformKh,CheckOutOfBounds]

pset.execute(kernels, runtime=delta(days=180), dt=delta(days=3),output_file=output_file) 

[erikvansebille]

Dear @cgrimaldi, thanks for your question. Unfortunately, there is no out-of-the-box option to avoid the DiffusionUniformKh from moving particles to land. You will have to write your own kernel that controls that 'behaviour'. You can write a kernel that tests whether a particle is moved to land by diffusion, and if so don't do the diffusion? Or reflect the particle? Or half the tilmestep?

Sorry I can't be of more help now; in my team we typically don't use diffusion anymore. Perhaps others can weigh in how they do this?