Failed advection past east-west domain boundaries in global HYCOM

[selipot]

I am releasing particles in a global HYCOM simulation but I am encountering an issue where Parcels is not able to find the next position of particles that would be advected east or west past the domain boundaries. I am attaching a figure that shows the longitude values on the 7055 by 9000 grid of my simulation. In my case, particles are not advected across the 1 or 9000 X-positions in the Indian Ocean and I get a warning at execution. I expected Parcels to be able to deal with the cyclic boundary condition. I am defining my fieldset from NetCDF files that contains the velocity fields and the Latitude(Y,X) and Longitude(Y,X) arrays with coordinates Y from 1 to 7055 and X from 1 to 9000.

Is there something else I should specify explicitly in defining my fieldset or kernel? Any advice would be greatly appreciated!

[selipot]

It looks like there's a solution already, thank you @philippemiron for the pointer! https://nbviewer.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_periodic_boundaries.ipynb

[selipot]

The tutoral_periodic_boundaries.ipynb is a start but I don’t think it is readily applicable when (i) the longitude at the end of the domain varies as a function of latitude, and (ii) when the limits of the domain do not correspond to the min and max of the longitude (0,359) or (-180,179) as an example.

With respect to (i), in that same tutorial, the kernel periodicBC conducts a comparison between the particle longitude and a boundary longitude but in my case this fails because the boundary longitude, the halo I believe, is actually a vector ... and with respect to (ii), many particles will have longitude greater or less than the longitude at the boundary but will still be in the domain (because of the way it is organized as shown above). Would a solution be to track the X, Y indices instead of longitude?

[erikvansebille]

Hi @selipot. I think that @VictorOnink has used global HYCOM fields for particle advection (also across the domain boundary) is his paper at https://iopscience.iop.org/article/10.1088/1748-9326/abecbd. Perhaps you can find some hints how he did it in his code at https://github.com/VictorOnink/Modeling-Global-Plastic-Beaching?

[selipot]

@erikvansebille Thank you for the pointer. Looking at the code I can't tell the type of HYCOM files used? @VictorOnink could you elaborate? I am also seeing that the boundary condition is a simple comparison to +/- 180 longitude.

[VictorOnink]

Hi @selipot, I used the HYCOM files from the GOFS 4.1 global reanalysis here: https://www.hycom.org/dataserver/gofs-3pt1/reanalysis. In this case, a simple comparison of +/- 180 degrees was sufficient to fix any BC issues I came across. I'm not entirely sure if I understand the issue you are running into, but if the problem is that your HYCOM files have longitude coordinates [0, 359], then same approach should also work.

[selipot]

Thank you for your reply @VictorOnink. As I tried to explain above and shown in the figure, the boundary longitude of the HYCOM grid I use is a function of latitude, and the east-west boundaries of the domain do not correspond to (-180,180) or (0,360) but instead to approximately 47 degrees. So a kernel that compares the particle position to +/- 180 or (0,360) to wrap around the particles does not work.

I am thinking of a hack that would consist of tweaking the Fieldset to make the longitude constant at the boundaries and also offset all the longitudes to place them between -180 and 180. @erikvansebille do you know if I can modify the fieldset values at the setting up stage or does Parcels re-read the Latitude and Longitude field each time it accesses a new Netcdf file?

Alternatively I could try to recast all my HYCOM files on an offset and regular grid? That would require to know how to go from an irregular to a regular grid ...

[VictorOnink]

Ok, then I get it I think. @selipot, have you tried maybe adding an additional lon/lat field to fieldset, which for each latitude would give you the longitude of the boundary? Your BC fixing kernel could then just sample what the boundary longitude is given the particle's latitude, and then use that to apply the BC.

[selipot]

Something like that could work. I am also thinking of adding an X field varying zonally from 1 to 9000 which would just be the coordinate along the longitude dimension. An additional kernel would test if the particle sample something larger than 9000 or smaller than 1 then wrap it around as needed ...

[erikvansebille]

Have you made progress with solving this issue, @selipot? Anything I can further help with?

[selipot]

Thank you for checking in. I am at the stage of trying to implement my proposed solution above, that is adding a fieldset with the X coordinate to use for the cyclic boundary condition ... getting a bit confused with the halos and how to use two fieldsets but I'll ask again if I get stuck.

[selipot]

@erikvansebille I am hitting a wall trying to implement my proposed solution. I am not even at the stage of implementing the cyclic boundary condition but trying to add a constant field to my fieldset to use for the boundary condition. I have put a notebook with the execution error at

https://github.com/selipot/hycom-unleashed/blob/master/particles_hycom_equator_to_share.ipynb

You'll see that I am trying to add a constant field to a fieldset containing global velocities from one set of files and global ssh from another set of files. Any suggestion or advice would be greatly appreciated!

[selipot]

Thank you @erikvansebille for finding a solution to this issue. The HYCOM output that I am using is such that the last longitude of the grid (to the east) is less than the first longitude of the grid (to the west) and the fieldset is lacking the needed data to wrap around (for my lack of better explanation), i.e lon[-1]<lon[0]. A solution that seems to work is to extend the grid to the east as follows:

lon_c = np.expand_dims(fieldsetH.U.grid.lon[:, 0], axis=1) fieldsetH.U.grid.lon = np.hstack((fieldsetH.U.grid.lon, lon_c))