Particles stuck when small velocity

[mathildejutras]

I use velocities in a 1 degree grid to advect particles. Particles stop moving away from land. By outputting the velocities, I could see that this occurs when the velocities are very small (of the order of 10^-8 m/s), and this occurs in the transitions from positive to negative velocities, which, when interpolated, provides small velocities (see figure). It causes the particles to spiral into a smaller and smaller spiral until they stop (see figure).

I thought increasing the time step would solve the problem, by preventing this spiraling. However, for any time step above 10 minutes, I do not see any significant difference in the generated track. Is there some a kind of threshold that controls the time step, changing it for a default value in some cases? I don't see why the tracks would look the same for different time steps.

Here's my script:
def periodicBC(particle, fieldset, time):
if particle.lon < fieldset.halo_west:
particle_dlon += fieldset.halo_east - fieldset.halo_west
elif particle.lon > fieldset.halo_east:
particle_dlon -= fieldset.halo_east - fieldset.halo_west

class ocean_particle(JITParticle):
#add some variables
gamman = Variable('gamman', dtype=np.float32, to_write='once', initial=0.)
ugeo = Variable('ugeo', dtype=np.float32)
vgeo = Variable('vgeo', dtype=np.float32)

def SampleU(particle, fieldset, time):
(u, v) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
particle.ugeo = u
particle.vgeo = v

if (particle.ugeo==0) or (particle.vgeo==0) :
particle.delete()

#Setting up the fieldset we will create using our data
file_names ={'U':file,
'V':file}
variables = {'U':'ugeo',
'V':'vgeo'}
dimensions = {'lat':'lat',
'lon':'lon'}

print('Create particle set')
fset = FieldSet.from_netcdf(file,variables,dimensions)

# for period ic boundaries
fset.add_constant("halo_west", fset.U.grid.lon[0])
fset.add_constant("halo_east", fset.U.grid.lon[-1])
fset.add_periodic_halo(zonal=True)

# set initialization
print('Set up initialization')

pset = ParticleSet.from_list(fieldset = fset,
pclass = ocean_particle,
lat = Lats,
lon = Lons)

print('Execute run')

output_file = pset.ParticleFile(name = 'outputs_lagrangian/run_MLbase_%s.zarr'%name, outputdt=timedelta(days=30))

kernels = [AdvectionRK4, periodicBC, SampleU]

pset.execute(kernels,
runtime = timedelta(days = 365*30),
dt = timedelta(minutes = 20),
output_file=output_file)

[erikvansebille]

Hi @mathildejutras, thanks for raising this. In our group, we typically see this kind of spiralling into trapped points when the flow field is not evolving in time. If the flow-field is stationary, then particles will stay on streamlines.

Do you have only one velocity snapshot in your netcdf file?

[mathildejutras]

Hi, thank you for the fast answer. Yes, I do, the flow-field is stationary.