Sea bottom interaction

[Oklyagina]

Dear Knut-Frode, good afternoon! I am wondering, how does OpenDrift handle with sea bottom topology?

The question is here because I met the following issue: when particles are released, for example, at depth of -200m in the region, where bottom is at -100m, they won't move (info about velocities there is filled with NaNs). But if particles are released at -300m in the same region, for some reason they move (see the attached picture, latitudes 43-44).

I tried to enable and disable filling NaNs during interpolation, but it didn't affect the result much. Could you please advice something on this?

[knutfrode]

Hi,

Elements seeded below seafloor, should be lifted to seafloor, as long as config setting general:seafloor_action is lift_to_seafloor, and thereafter they should drift like normal.

Anyway, what you observe seems strange in any case I tried to reproduce with a minimalistic example, but did not see something similar. One thing to notice is that your plot shows land-colors where your particles are, indicating land, and not water, so I guess you are disabling the GSHHG landmask with o.set_config('general:use_auto_landmask', False), as this does not include Caspian Sea? Anyway, I cannot sea how this should make a difference here.

It could be something with the ocean current reader you are using, can you please provide the output of

>>> print(<ocean_reader>)

and perhaps also your full script.

[Oklyagina]

Thank you for the fast answer! Yes, you're right, I had to disable auto landmask.

Here is the output from reader:

and the script:

`import numpy as np
from datetime import datetime, timedelta
from opendrift.models.oceandrift import OceanDrift
from opendrift.readers import reader_shape
from opendrift.readers import reader_netCDF_CF_generic

reader_norkyst =  reader_netCDF_CF_generic.Reader('data_z.nc')
import cartopy.io.shapereader as shpreader
shpfilename = shpreader.natural_earth(resolution='110m',
                                     category='physical',
                                     name='coastline')
reader_natural = reader_shape.Reader.from_shpfiles('/shorelines/GSHHS_shp/h/GSHHS_h_L2.shp', invert = True)

lons = np.linspace(48, 52, 5)
lats = np.linspace(41, 44, 5)

print(lons,lats)
lons, lats = np.meshgrid(lons, lats)
lons = lons.ravel()
lats = lats.ravel()

print(lons,lats)
depth = np.linspace(0,-400,5)
#depth = [-300]
print(depth)
time_step = 5000

for z in depth:
    days = 180
    o = OceanDrift(loglevel=20)  # Set loglevel to 0 for debug information
    o.set_config('seed:ocean_only', True)
    o.set_config('general:use_auto_landmask', False)
    o.set_config('general:coastline_action', 'previous')
    o.set_config('drift:advection_scheme', 'runge-kutta')
    o.add_reader(reader_norkyst)
    o.add_reader(reader_natural)

    print(reader_norkyst)
    print(o)
    #reader_norkyst.buffer = 300
    o.seed_elements(lons, lats, radius=0, number=25, z = z,
                time=reader_norkyst.start_time)

    filename2 = 'z_coord/middle/map_lon =' + str(np.mean(lons)) + '_lat =' + \
                str(np.mean(lats)) +  '_mod_time_' + str(days) + 'depth' + str(z) + '.nc'
    o.run(time_step=time_step, duration=timedelta(days=days), outfile=filename2,
          export_variables=['lon', 'lat', 'z'])  # , outfile=ncfile)

    filename = 'z_coord/middle/map_lon =' + str(np.mean(lons)) + '_lat =' +\
        str(np.mean(lats)) + '_mod_time_' + str(days) + 'depth' + str(z) + '.jpg'
    o.plot(buffer=0.2, show_particles=True, filename=filename, fast=True, hide_landmask=False, linecolor='z')

[knutfrode]

Some things to notice here:

• "norkyst" is a specific ocean model by MET Norway, so you'd probably like to call your reader reader_current or something like this, to avoid confusion
• Since you are not using shpfilename you could delete this line
• Since your landmask seems to be inverted on the plot, you should perhaps change inverted=True to inverted=False, though then the simulation might perhaps not work?
• You are running the simulations in a loop over all depths. This will give a large amount of overlap and performance loss. You could instead simply do only the seeding over the same loop. Or even the seeding might be done in a single call, using meshgrid also over z:

  lons = np.linspace(48, 52, 5)
  lats = np.linspace(41, 44, 5)
  zs = np.linspace(0, -20, 5)
  lon, lat, z = np.meshgrid(lons, lats, zs)
  o.seed_elements(lon=lon.ravel(), lat=lat.ravel(), z=z.ravel(), time=time)

But the main thing here might be the fact that you seem to have no reader for the seafloor depth, thus OpenDrift will here assume a depth of 10.000m, as there is no built-in global seafloor depth dataset. I suspect this is the reason for the strange behaviour. Could you add this to your netCDF input file, or add another reader containing this for the Caspian Sea? It should have standard_name: sea_floor_depth_below_sea_level If adding this to your output file, you might also consider adding a landmask (land_binary_mask), to be used instead of the shapefile.

[Oklyagina]

Thank you for the advice, especially for about meshgrid, I was looking for this solution! Yes, the simulation does not work, if the shapefile is inverted.

Sorry for the stupid question: Is there any way to make a binary mask via OpenDrift only? When experimenting, I tried to pass the land_binary_mask like other variables from reader (like variable['land_binary_mask']). This variable is not in the .nc file, so I just made a masked array from the velocities on z = 0. The model didn't feel it and asked to provide a landmask file.

I provided sea_floor_depth_below_sea_level info, it works correctly now! Thank you!

[knutfrode]

Regarding the question of making binary mask from OpenDrift only: OpenDrift will never modify the reader input files, but there is a mechanism to derive some environment variables from others. E.g. here is a method to derive land_binary_mask from sea_floor_depth_below_sea_level (landmask is 0 where depth is <= 0). https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/basereader/variables.py#L381

This one is not activated by default, so to use it, you have to type

<your_reader>.activate_environment_mapping('land_binary_mask_from_ocean_depth')

It is also possible to modify or add new such methods. E.g. one could say that land_binary_mask is 1 wherever current speed is 0, although this is a bit dangerous.

[Oklyagina]

Sorry I didn't answer earlier, but this method worked. Thank you!