Algorithm for detecting grounding/beaching?

[selipot]

I am currently dealing with some numerical drifters released with Parcels in a model HYCOM run and I want to assess a posteriori when my drifters actually become “grounded” or "beached" in the model. Because I did not specify any specific boundary condition in my experiment (sort of ignored that aspect), I think my particles become grounded by infinitesimally slowing down towards coasts. Has anyone devised an algorithm to flag such grounding? I am currently looking at an algorithm where I detect grounding if the particle speed is continuously less than a threshold value (0.003 m/s) for a number of days (1). Any suggestion?

[erikvansebille]

Hi @selipot, have you seen the tutorial that @reint-fischer has written, on grounding/beaching? See https://docs.oceanparcels.org/en/latest/examples/documentation_stuck_particles.html. That could provide some useful background?

As for which condition to use to flag beaching: ideally you'd use a mask file from the model that differentiates between land and ocean points (and use the right type of interpolation). If you don't have such a mask feel or can't make it, then a condition on the velocity could work, indeed. But that's always a bit hacky...

[selipot]

Ok so my velocity fields are equivalent to an A grid so I can simply specify the boundary condition in the definition of my fieldset as in the tutorial example and as follows?

interp_method={
        "U": "partialslip",
        "V": "partialslip",
    },  # Setting the interpolation for U and V

Somehow I never saw that in the API reference ...

[erikvansebille]

Yep, correct. Let us know whether this works!

[JamiePringle]

@selipot You mention that you are using HYCOM, which is a c-grid model. If your date is, as you state above "effectively an a-grid" that means that there has been averaging of velocity onto the rho-grid. In that case, it is very problematic to compute beaching.

The underlying reason this is true is that the no-normal-flow-through-land boundary is implemented in a c-grid by setting the velocity on u and v grids to zero at the land edge. When you average to an A grid, you will either loose this boundary condition or, if hacking it in, loosing the non-divergence of the flow field. Either leads to a bunch of spurious grounding and/or non-volume-density-of-floats conservation near coasts. With the C-grid grounding should be very rare, and due to numerical truncation issues.

Even if your particles are fixed depth or otherwise non-Lagrangian, so you could expect grounding, you will tend to get strange results and greatly increased grounding when using C-grid averaged to A-grid data.

Be very careful with this.

Jamie

[LauraBachmaier]

I asked about the problem of particles travelling on land a couple months ago and was advised to use C grid data. However, C grid data is extremely hard to come by and I have continued working with A grid data as, even when using C grid data, particles seemingly are still ending up on land due to a horizontal diffusivity field that I have added. To deal with this I have added a Boolean landmask. But it seems due to the way the A grid is structured the simulation isn't incorporating the landmask with my fieldset properly or not sampling from it properly.

To deal with this I have tried using

fieldset.add_field(Field("landmask", data=landmask, lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat, mesh="spherical", interp_method='nearest'))

which seems to make it work properly. I just wanted to ask if anyone could explain to me what the differences are between interp_method = 'nearest' vs. linear vs. linear_invdist_land_tracer vs. cgrid_tracer and when to use which one in what instance. I'm adding a number of other additional fields to my fieldset (e.g. a unbeaching field). Would all additionally added fields need to use an interpolation method when working with an A grid?

I have also tried using the above

interp_method={ "U": "partialslip", "V": "partialslip", }, # Setting the interpolation for U and V

to avoid having to add a landmask altogether but this does not seem to work for me for some reason. Whether I add a slip boundary condition or not, particles end up on land.

Could someone please elaborate what the differences are between these interpolation methods. And if I use a partial or full slip, would I still need to use one of the other interpolation methods for adding new fields? I would really appreciate any help as I am getting more confused the more I read about interpolation methods. Many thanks :)

[erikvansebille]

Hi @LauraBachmaier, thanks for your question. The differences between the interpolation methods is explained in this tutorial.

Not sure why the slid boundary condition does not work for you. Could it be because of one of the other kernels? Perhaps try not using these, and see if the slip boundary condition works as expected when you're only using AdvectionRK4?