Problem with FieldOutOfBoundSurfaceError when sampling 2D field

[soizicg]

Hello,

I use a 3D fieldset with time varying vertical coordinates and wanted to develop a kernel called "CheckBoundaries" to handle particle interactions with model boundaries (surface, sea bed, land...), but the status code doesn't seem to be working as expected.

If I sample a 2D field in a kernel, the particle will be able to escape through the surface without raising an error, however if I sample a velocity field, the out of bound at the surface error will be raised when the particle reaches negative vertical coordinates (as expected).

Here is the code to define the fieldset:

from parcels import(FieldSet,
JITParticle,
ParticleSet,
AdvectionRK4_3D,
Field,
Variable,
StatusCode
)

scenario = "may2019"
forcing_filename = "../forcing/" + scenario + "_dvm_8weeks_test.txt" 
netcdf_files = open(forcing_filename,'r').read().split('\n')
f = netCDF4.Dataset('F:/NEATL_w_Kz/coordinates.nc')
landmask = f.variables['mask_psi'][:].data

variables = {"U": "u", "V": "v", "W": "w", "depth_w": "z_w", "Kz": "Kz"}

grid_dimensions = {
    "lon": "lon_psi", "lat": "lat_psi", "depth": "not_yet_set", "time": "time"}
dimensions = {
    "U": grid_dimensions,
    "V": grid_dimensions,
    "W": grid_dimensions,
    "depth_w": grid_dimensions,
    "Kz": {"lon": "lon_psi", "lat": "lat_psi","time": "time"}
}     

fieldset = FieldSet.from_netcdf(netcdf_files, variables, dimensions, gridindexingtype='nemo', mesh="spherical", allow_time_extrapolation=True)
fieldset.U.set_depth_from_field(fieldset.depth_w)
fieldset.V.set_depth_from_field(fieldset.depth_w)
fieldset.W.set_depth_from_field(fieldset.depth_w)
fieldset.add_field(Field(name="landmask", data=landmask, lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat, mesh="spherical"))

Then to define the particle set:

LobsterParticle = JITParticle.add_variables(  # It is a JIT particle JITParticle
    [
        Variable("is_land", initial=0, to_write=False),
    ]
)
    
pset = ParticleSet(
    fieldset=fieldset, # the fields that the particleset uses
    pclass=LobsterParticle, # define the type of particle
    lon=-5.6703,  # release longitude
    lat=51.7214, # release latitude
    depth=0,
    time=0 # release time
)        

output_file = pset.ParticleFile(
    name=("test.zarr"),
    outputdt=3600 #delta(minutes=10)    
)

Finally the kernel definition (sampling of the landmask) and execution of the particle simulation:

def CheckBoundaries(particle, fieldset, time):
    is_land = fieldset.landmask[0,0,particle.lat, particle.lon]

pset.execute(  
    [AdvectionRK4_3D, CheckBoundaries],
    runtime=2 *60*60*24,  #length of model run
    dt=600,  #timestep of model
    output_file=output_file,
)

Now when I run this code, the particle is able to pierce the surface, here is the output of the z variable of the particle:

0
0.600053250789642
1.18962800502777
1.61126649379730
1.75607454776764
1.55384206771851
1.09453940391541
0.509002625942230
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266
-0.0293525625020266

However if I change the kernel and sample u instead of the landmask that way:

def CheckBoundaries(particle, fieldset, time):
    particle.is_land = fieldset.U[particle]

pset.execute(  
    [AdvectionRK4_3D, CheckBoundaries],
    runtime=2 *60*60*24,  #length of model run
    dt=600,  #timestep of model
    output_file=output_file,
)

Oceanparcels now raises a FieldOutOfBoundSurfaceError as the particle go through the surface, then the simulation is stopped.
The particle z variable doesn't become negative due to the status code being raised:

0
0.600053250789642
1.18962800502777
1.61126649379730
1.75607454776764
1.55384206771851
1.09453940391541
0.509002625942230

This problem with the status code seems to occur only when I sample a 2D field (I have also tried with the Kz field that I defined earlier). Should I sample 2D field in a different way?
Thanks for the help.

[erikvansebille]

Hi @soizicg, thanks for raising this question. However, I'm not sure I understand the problem. In Parcels, 2D fields don't consider depth at all, so there is no concept of a FieldOutOfBoundSurfaceError for a 2D field, because there is no depth dimension. It's as if these 2D-fields stretch infinitely in the z-dimension

Of course, you could raise such an error yourself, with a simple if-statement, see the tutorial on Status Codes. Something like

if particle.depth < 0:
    particle.state = StatusCode.ErrorThroughSurface

[soizicg]

Hi @erikvansebille, sorry for the late answer, I was away.
In my kernel list I have "AdvectionRK4_3D" which samples the UVW fields to get the particle position at next step, so I was expecting to get an error even if I was sampling a 2D field in my "CheckBoundaries" custom kernel. But from your answer I guess the status code only raise errors on the custom kernels and not the kernels already implemented within parcels?

[erikvansebille]

Hi @soizicg, sorry if my previous reply was not clear; but this is not a matter of built-in versus custom kernels (under the hood, parcels doesn’t care at all), but about 2D vs 3D Fields. If you have a 2D Field, parcels will never throw an FieldOutOfBoundSurfaceError. You will have to create this error yourself, if you want. Hope that clarifies?

[soizicg]

Ok so I am not sure to understand, I am using 3D fields (which is why I use the built-in kernel "AdvectionRK4_3D").
However the landmask is a 2D field because the model grid is using terrain-following sigma coordinates, therefore I was sampling a 2D field in my custom kernel.
Can the status code also work in the built-in kernel if the particle is escaping through the surface?

[erikvansebille]

The interpolation is per Field, so particles can't be deleted based on your 2D landmarks. Sorry, but I can't make it clearer than that... You really have to come up with another solution than the FieldOutOfBoundSurfaceError for your problem, it seems