Keeping particles off of land for lake simulation

[bwahyatt]

Hi all,

Hope this is the correct place to ask this question. I am currently using parcels to track particles for a lake simulation model. For now, I have been writing the output of the particle trajectories to .csv files which I am then reading into an independent animation script which uses the velocity field data outputted from an MITgcm executable. This U,V,W velocity field data is then used to create the .nc files that are loaded into my parcels script. These animations reveal that particles can temporarily cross the "shore" of the lake and travel on the land.

Is there a good way to prevent this effect from happening? I am looking to have the particles move along the edge of the shore when they encounter it, rather than crossing over and onto the land. In my case the bathymetry of the simulation is very simple, just a rectangular peninsula in a rectangular lake domain. I have been using this function as a quick fix to correct the trajectories:

def stop(particle,fieldset,time):
  if particle.lon > 7960:     # max_lon
    particle.lon = 7960  
  elif particle.lon < 20:     # min_lon
    particle.lon = 20

  if particle.lat > 7960:     # max_lat
    particle.lat = 7960
  elif particle.lat < 20:     # min_lat
    particle.lat = 20
            
  #keep this as is for now
  if particle.depth > 15.7:   # max_depth
    particle.depth = 15.7
  elif particle.depth < 0.1:  # min_depth
    particle.depth = 0.1

  if particle.lon < 3860 and particle.lat > 3040: # east_edge and south_edge
  #if you're coming onto land from the east: it should be just under x~3880
    if (particle.lon + 30) > 3880: 
      particle.lon = 3860                           # east_edge
    #or if it's from the south
    elif (particle.lat - 30) < 3060:
      particle.lat = 3040                           # south_edge
    else:
      print('Particle is too far inland')

Which I use by adding to the kernel in pset.execute, i.e. pset.execute(AdvectionRK4+pset.Kernel(stop), ...). The hand written values in this function are relevant to my specific bathymetry.

Is there a more general way to go about fixing this? I.e. how should I approach this for a bathymetry with an arbitrary/realistic geometry? I know the documentation mentions using a land=True parameter, but in this case I am using a flat mesh/cartesian coordinates with a custom made bathymetry (not modelling real oceans/continental coastlines). The bathymetry data is stored in a .bin file, and I have access to other independent scripts which use it produce cross sections of the lake as numpy arrays.

For additional context, I am a 2nd year physics undergrad and very new to parcels, not a veteran programmer by any means (feel free to dumb it down, I will not be offended). Thanks so much!

[reint-fischer]

Hi @bwahyatt ,

This is indeed a good place to ask these types of questions! Particles crossing the "shore" or other solid-fluid boundaries can be a complex problem, depending on the types of fields you are working with and in many cases we have not found a general 'answer' for all purposes. In fact, many simulations have used similar solutions to the one you have implemented.

I am working on a notebook to explain the different options in different cases. Since you are working with MITgcm, am I correct to assume that the velocities come on a staggered grid: https://mitgcm.readthedocs.io/en/latest/algorithm/horiz-grid.html?

If your velocities are indeed defined on this staggered grid (C-grid), the velocities at the cell edges that form your "shore" should be zero and particles should not be able to cross them. If they still do, could you either try decreasing your computational timestep (dt) or share the way you have set up your fieldset and the resulting animation?

If your velocities are defined on an unstaggered grid, could you confirm this and maybe show how/where your particles cross the "shore"? Then I can recommend some options based on that.

Cheers,

[bwahyatt]

Hi @reint-fischer

Yes I have confirmed that the velocities from MITgcm are on a C-grid.

For context, I am running a simulation for a 24h duration. For the MITgcm executable, the "computation" time step is 1 second and the "output" time step is currently 300s. I.e. the velocity field info is written to file every 300s, so there are 289 files each for U,V and W data. These velocity field files are then condensed into .nc files (289 of these too, with info every 300s). These .nc files are then read into my parcels .py script.

In my parcels script, I am using the "output" time step of 300s for my output_file "outputdt". For a while I was using this same time step as my dt time step in pset.execute, but I have also tried dt=1 to match the "computation" time step above but it didn't seem to make much difference. i.e. my script has:

time_step = int(fieldset.U.grid.time[1])  
nc_fname = "parcels_output_"+d_string+str(time_step)+"s.nc"

output_file = pset.ParticleFile(name=nc_fname, outputdt=time_step)

pset.execute(AdvectionRK4_3D,
             runtime=duration,                     
             dt=1,                                  
             output_file=output_file)

But should I have dt=1 or any dt < 300? For instance If the .nc files I am loading into the script only have data every 300s, will having dt=1 compute the particle trajectories on a "static" velocity field for the first 300 steps, then the field changes, and the trajectories are computed for another 300 time steps on the updated field?

This is what my fieldset looks like:

filenames = {"U":{"lon": "output_0000000000.nc", "lat": "output_0000000000.nc", "depth": "output_0000000000.nc", "data": "output*.nc"},
             "V":{"lon": "output_0000000000.nc", "lat": "output_0000000000.nc", "depth": "output_0000000000.nc", "data": "output*.nc"},
             "W":{"lon": "output_0000000000.nc", "lat": "output_0000000000.nc", "depth": "output_0000000000.nc", "data": "output*.nc"}}
             
variables = {"U":"U", "V":"V", "W":"W"}
dimensions = {"U": {"lon": "x", "lat": "y", "depth": "z", "time": "time"}, 
             "V": {"lon": "x", "lat": "y", "depth": "z", "time": "time"},
             "W": {"lon": "x", "lat": "y", "depth": "z", "time": "time"}} 
              
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True, mesh="flat")

All of the files I am loading in are have a numeric suffix which is the time stamp of their data (e.g. I have output_0000000000.nc, output_0000000300.nc, output_0000000600.nc, .... output_0000086400.nc).

The animation is done by writing the particle trajectories to csv files (with their x,y,z,t coordinates at every time step) which are read into a separate animation script (independent of parcels). I have attached a png of a snapshot of a particle at its most inland point.

The particle (red dot) on the left is at x=3835.7m. The domain is a 400x400 grid, with dx and dy = 20m (8km x 8km in physical units) with a roughly rectangular peninsula in the NW corner. I know from the bathymetry file that in the x direction the last grid point with land is indexed at 192 (i.e. x=3840m in my case) and the water begins at 193 (x=3860, the value in my "stop" function in my initial post).

Any advice on how to proceed from here? Thanks!

[reint-fischer]

hey @bwahyatt,

To start with your first question: you can indeed have dt=1s or any dt<300s, as parcels will linearly interpolate in time between the two files. In many models it is required for numerical stability to comply with a Courant-Friedrichs-Lewy condition dt <= dx / u_max, which bounds the fluid, or in our case the particles to move less one grid cell every timestep. With dx = 20m and u_mx around 0.3 m/s (?), it would be good to keep dt<<300s.

Then for the main problem, parcels requires a specific set up for C grids: This tutorial shows how you can define the dimensions dictionary for C grids. There are two options:

1. You can keep using fieldset.from_netcdf and specify the arguments interp_method = 'cgrid_velocity' and gridindexingtype='mitgcm'.
2. You can use fieldset.from_c_grid_dataset and specify gridindexingtype='mitgcm'.

Let me know if that works and good luck with your simulations!

[bwahyatt]

Hi @reint-fischer ,

Thank you for your response. I was hoping to clarify some things about the tutorial you linked above on initializing dimensions.

In/above the In [4] block, it says that "we need to provide the coordinates of the f-points. In NEMO, these are called glamf, gphif and depthw (in MITgcm, these would be called XG, YG and Zl)." The code then initializes dimensions with glamf, gphif and depthw for lat, lon and depth. However, parcels' documentation says that dimensions is a dictionary "mapping data dimensions (lon, lat, depth, time, data) to dimensions in the netCF file(s)"

I am currently using a script called datnetcdf.py (not written by me) to create .nc files from my MITgcm velocity field data to read into parcels. I have attached the whole script to this post as a txt file if you wish to look, but in particular there it contains these lines which read in the XG/YG files (and RC for the z direction):

XC = mgu.rdmds('XG')
YC = mgu.rdmds('YG')
RC = mgu.rdmds('RC')

depth = mgu.rdmds('Depth')

myX = np.tile(XC.reshape(1, Ny, Nx), (Nz, 1, 1))
myY = np.tile(YC.reshape(1, Ny, Nx), (Nz, 1, 1))
myZ = np.tile(RC.reshape(Nz, 1, 1), [1, Ny, Nx])

And later on when the .nc files are being created, these lines appear:

fp = Dataset(outfile, 'w', format='NETCDF4')

# Create dimension objects
x_dim = fp.createDimension('x',Nx)
y_dim = fp.createDimension('y',Ny)
z_dim = fp.createDimension('z',Nz)
t_dim = fp.createDimension('time',1)

# Create and write grids
x_grid  = fp.createVariable('x',    np.float64, ('x',))
y_grid  = fp.createVariable('y',    np.float64, ('y',))
z_grid  = fp.createVariable('z',    np.float64, ('z',))
zc_grid = fp.createVariable('zc',   np.float64, ('z','y','x'))
t_grid  = fp.createVariable('time', np.float64, ('time',))

x_grid[:]      = myX[0,0,:]
y_grid[:]      = myY[0,:,0]
z_grid[:]      = -1*np.array(RC)

I suppose my question is, should my dimensions dictionary use the "dimension" objects of these .nc files (like the documentation suggests) or the "grids" (as the tutorial you linked suggests)?

Or is this datnetcdf.py script insufficiently handling the grids/XG etc files? It seems like that information is included in the .nc files it is creating (but it is worth reiterating this is all pretty new to me!)

Thanks Reint

datnetcdf.txt

[reint-fischer]

Hi @bwahyatt,

To read C grids correctly, parcels needs the "grid" coordinates of the f-points of your cells. Since you have a 'flat' domain, these coordinates, which are originally 3D (for possible non-rectilinear grids):

XC = mgu.rdmds('XG')
YC = mgu.rdmds('YG')
RC = mgu.rdmds('RC')

depth = mgu.rdmds('Depth')

myX = np.tile(XC.reshape(1, Ny, Nx), (Nz, 1, 1))
myY = np.tile(YC.reshape(1, Ny, Nx), (Nz, 1, 1))
myZ = np.tile(RC.reshape(Nz, 1, 1), [1, Ny, Nx])

can be collapsed to 1D:

x_grid[:]      = myX[0,0,:]
y_grid[:]      = myY[0,:,0]
z_grid[:]      = -1*np.array(RC)

It seems to me that your postprocessing in 'datnetcdf.txt' stores these values (XG, YG) in the netCDF variables 'x' and 'y':

# Create and write grids
x_grid  = fp.createVariable('x',    np.float64, ('x',))
y_grid  = fp.createVariable('y',    np.float64, ('y',))
z_grid  = fp.createVariable('z',    np.float64, ('z',))
zc_grid = fp.createVariable('zc',   np.float64, ('z','y','x'))
t_grid  = fp.createVariable('time', np.float64, ('time',))

x_grid[:]      = myX[0,0,:]
y_grid[:]      = myY[0,:,0]
z_grid[:]      = -1*np.array(RC)

[bwahyatt]

Hi @reint-fischer ,

Firstly I edited my datnetcdf.py script so that the grid variable names are 'xg', 'yg' etc instead of just 'x', 'y' etc to distinguish them from the "dimensions" variables defined above. Having said that, I reran this script to update my netcdf files and load them back into parcels. However, I am still having particles travelling pretty far inland in my animations (reaching x=3812m, instead of the land boundary x=3840m).

So my dimensions dictionary has:

dimensions = {"U": {"lon": "xg", "lat": "yg", "depth": "zg", "time": "timeg"}, 
             "V": {"lon": "xg", "lat": "yg", "depth": "zg", "time": "timeg"},
             "W": {"lon": "xg", "lat": "yg", "depth": "zg", "time": "timeg"}} 

I tried using the 2 fieldset options you mentioned above:

fieldset = FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True, mesh="flat", interp_method='cgrid_velocity', gridindexingtype='mitgcm')

fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions, allow_time_extrapolation=True, mesh='flat', gridindexingtype='mitgcm')

I also tried using from_mitgcm which I happened to come across in the documentation:

fieldset = FieldSet.from_mitgcm(filenames, variables, dimensions, allow_time_extrapolation=True, mesh="flat")

But all three of these choices are producing the same trajectories. Any suggestions? Thanks