ThroughSurfaceError

[ajithacyriac]

Dear Ocean Parcels team,

I am new to Ocean Parcels and Python. I have been trying to use parcels on 3D velocity data following the example notebook using NEMO data (https://nbviewer.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_3D.ipynb). When I try to execute pset, I keep getting the same error, 'ThroughSurfaceError'. The code that I use is as follows:

from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4_3D
from glob import glob
from datetime import timedelta as delta
from parcels import logger, XarrayDecodedFilter
logger.addFilter(XarrayDecodedFilter()) # Add a filter for the xarray decoding warning

filenames = {'U':'BRAN_3Dcurrents_2015.nc',
'V':'BRAN_3Dcurrents_2015.nc',
'W':'BRAN_3Dcurrents_2015.nc'}
variables = {'U': 'U',
'V': 'V',
'W': 'W'}
dimensions = {'U': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'},
'V': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'},
'W': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'}}
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
pset = ParticleSet.from_line(fieldset=fieldset, pclass=JITParticle,
size=10,
start=(130, -35), # releasing on a line: the start longitude and latitude
finish=(132, -35), # end longitude and latitude
depth=0)
kernels = pset.Kernel(AdvectionRK4_3D)
pset.execute(kernels, runtime=delta(days=4), dt=delta(hours=1))

I am getting the following error while I run the above codes:

---

ThroughSurfaceError Traceback (most recent call last)
Cell In [7], line 2
1 kernels = pset.Kernel(AdvectionRK4_3D)
----> 2 pset.execute(kernels, runtime=delta(days=4), dt=delta(hours=1))

File ~\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\particleset\baseparticleset.py:450, in BaseParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, moviedt, recovery, output_file, movie_background_field, verbose_progress, postIterationCallbacks, callbackdt)
448 # If we don't perform interaction, only execute the normal kernel efficiently.
449 if self.interaction_kernel is None:
--> 450 self.kernel.execute(self, endtime=next_time, dt=dt, recovery=recovery, output_file=output_file,
451 execute_once=execute_once)
452 # Interaction: interleave the interaction and non-interaction kernel for each time step.
453 # E.g. Inter -> Normal -> Inter -> Normal if endtime-time == 2*dt
454 else:
455 cur_time = time

File ~\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\kernel\kernelsoa.py:227, in KernelSOA.execute(self, pset, endtime, dt, recovery, output_file, execute_once)
225 recovery_kernel = recovery_map[p.state]
226 p.set_state(StateCode.Success)
--> 227 recovery_kernel(p, self.fieldset, p.time)
228 if p.isComputed():
229 p.reset_state()

File ~\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\tools\statuscodes.py:207, in recovery_kernel_through_surface(particle, fieldset, time)
203 """Default sampling error kernel that throws OutOfBoundsError"""
204 if particle.exception is None:
205 # TODO: JIT does not yet provide the context that created
206 # the exception. We need to pass that info back from C.
--> 207 raise ThroughSurfaceError(particle, fieldset)
208 else:
209 error = particle.exception

ThroughSurfaceError: 0
Particle P[10](lon=130.000000, lat=-35.000000, depth=0.000000, time=0.000000)
Time: 2015-01-16T12:00:00.000000000, timestep dt: 3600.000000
Through-surface sampling by particle at (130.000000, -35.000000, 0.000000)

I would be very thankful if you could help me with it.

[erikvansebille]

Dear @ajithacyriac, could it be that the top layer in your velocity fields is not at 0m, but slightly deeper? If so, a simple solution would be not to specify the depth=0 when you create the ParticleSet. If you don't specify depth, then Parcels will by default use the depth of the top-most layer

[ajithacyriac]

Dear @erikvansebille, thanks for the response. I get the same error when I put "depth=2.5", which is the top layer of my velocity field. I tried without the depth option and still got the same error. But when I put "depth=10m", which is a deeper layer in my field, it works. I wonder why.

[erikvansebille]

Hmm, could it be that BRAN is a B-grid? Then you might want to use `FieldSet.from_b_grid_dataset instead? Also note that in some B-grids the upper layer is on the bottom of the cell; see also #573 and #578.

[ajithacyriac]

I tried with `FieldSet.from_b_grid_dataset'. It still gives me error for depth values less than 10m. Could there be a physical explanation for it? Do the particles need to be released at a certain depth which is deep enough for them to be advected? But it is not a problem for the example data in https://nbviewer.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_3D.ipynb.

[erikvansebille]

Hmm, that is strange... Could you send me some of the files (via transfer service/email?) so that I can try it myself?

[erikvansebille]

I checked your files, and think I know what's going on. You are using the AdvectionRK4_3D kernel, and your vertical velocities are so high that your particles actually move through fieldset.U.grid.depth[0] at 5m. So this is because of the velocity fields themselves, not because of Parcels. If you use AdvectionRK4 (so not the 3D-version) and don't specify a depth in your ParticleSet, then the particles move correctly

What you need to do is 'capture' these particles that go through the surface and e.g. push them back. Something like

from parcels import ErrorCode

def pushback(particle, fieldset, time):
    particle.depth = 10.

pset.execute(kernels, runtime=delta(days=4), dt=delta(hours=1), output_file=output_file, 
                       recovery={ErrorCode.ErrorThroughSurface: pushback})

[ajithacyriac]

Dear @erikvansebille, it worked!! Thank you so much!