Introduce new category of seafloor interaction

[nordam]

For modelling things like mineral particles or microplastics that sink, I would want to use a new type of seafloor interaction, where we reflect for diffusion, and deactivate when particles sink to the floor. One way to implement this is to use the following sequence of steps in the vertical transport:

• Displace randomly due to diffusion
• Reflect particles about both the sea surface and the sea floor (and handle the special case where it is so shallow that particles who reflect about one boundary end up outside the other boundary)
• Let particles rise/sink due to buoyancy
• Deactivate those particles that are below the floor (and set those that are above the surface to be exactly at the surface)

Just as a test, I tried implementing this by purely changing the interact_with_seafloor method as shown below (it's only the part below elif i == 'deactivate' that is changed), but I don't think this is the optimal way, since this effectively reverses some of the steps that have already been taken in order to tell apart the diffusion and the buoyancy.

def _interact_with_seafloor(self):
    import logging; logger = logging.getLogger(__name__)
    """Seafloor interaction according to configuration setting"""
    if self.num_elements_active() == 0:
        return
    if 'sea_floor_depth_below_sea_level' not in self.env.priority_list:
        return

    if not hasattr(self, 'environment'):
        logger.warning('Seafloor check not being run because environment is missing. '
                       'This will happen the first time the function is run but if it happens '
                       'subsequently there is probably a problem.')
        return

    if not hasattr(self.environment, 'sea_surface_height'):
        logger.warning('Seafloor check not being run because sea_surface_height is missing. ')
        return

    # the shape of these is different than the original arrays
    # because it is for active drifters
    sea_floor_depth = self.sea_floor_depth()
    sea_surface_height = self.sea_surface_height()

    # Check if any elements are below sea floor
    # But remember that the water column is the sea floor depth + sea surface height
    ibelow = self.elements.z < -(sea_floor_depth + sea_surface_height)
    below = np.where(ibelow)[0]

    if len(below) == 0:
        logger.debug('No elements hit seafloor.')
        return

    i = self.get_config('general:seafloor_action')
    if i == 'lift_to_seafloor':
        logger.debug('Lifting %s elements to seafloor.' % len(below))
        self.elements.z[below] = -sea_floor_depth[below]

    elif i == 'deactivate':
        # Only deactivate particles that were placed below the sea floor by sinking

        # First, calculate the position the particle would have after the diffusion step
        # by subtracting the vertical motion due to buoyancy
        z_diffusion = self.elements.z - self.time_step.total_seconds() * self.elements.terminal_velocity
        # For those particles that are below after diffusion,
        # reflect about sea floor
        mask = z_diffusion < -(sea_floor_depth + sea_surface_height)
        inds = np.where(mask)[0]
        self.elements.z[inds] = -2*(sea_floor_depth + sea_surface_height)[inds] - z_diffusion[inds]
        # Re-apply vertical motion due to buoyancy to those that were reflected
        self.elements.z[inds] = self.elements.z[inds] + self.time_step.total_seconds() * self.elements.terminal_velocity[inds]
        # Deactivate particles if the are below sea floor after buoyancy step
        mask = self.elements.z < -(sea_floor_depth + sea_surface_height)
        inds = np.where(mask)[0]
        self.deactivate_elements(mask, reason='seafloor')
        self.elements.z[inds] = -sea_floor_depth[inds]

    elif i == 'previous':  # Go back to previous position (in water)
        logger.warning('%s elements hit seafloor, '
                       'moving back ' % len(below))
        below_ID = self.elements.ID[below]
        self.elements.lon[below] = \
            np.copy(self.previous_lon[below_ID - 1])
        self.elements.lat[below] = \
            np.copy(self.previous_lat[below_ID - 1])

[nordam]

See sections 2.4, 3.2.1, and 4.2 in this paper, for further discussion of this type of boundary treatment: https://gmd.copernicus.org/articles/16/5339/2023/

[nordam]

On further searching, it looks to me like interact_with_seafloor is called twice in each step, after rising/sinking, and after diffusion:

• Line 427 of oceandrift.py:

  opendrift/opendrift/models/oceandrift.py

  Line 427 in 5bea84b

  self.interact_with_seafloor()

• Line 637 of oceandrift.py:

  opendrift/opendrift/models/oceandrift.py

  Line 637 in 5bea84b

  self.interact_with_seafloor()

Then I think a reasonable solution could be to add an argument to the interact_with_seafloor method, to indicate if it is called after rising/sinking, or after diffusion, and then allow it to handle the two cases differently. Say for example that the method took an argument stage, which could be either 'buoyancy' or 'diffusion'.

Edit 1:

Digging a bit deeper, I see that actually, interact_with_seafloor is called both in vertical_mixing and in vertical_buoyancy, but it looks like only one of those are called for each step in oceandrift:

opendrift/opendrift/models/oceandrift.py

Lines 191 to 196 in 5bea84b

	# Turbulent Mixing
	if self.get_config('drift:vertical_mixing') is True:
	self.update_terminal_velocity()
	self.vertical_mixing()
	else: # Buoyancy
	self.vertical_buoyancy()

I don't understand this, isn't it common to use both buoyancy and diffusion?

Edit 2: Ok, I see now that vertical_mixing contains both buoyancy and diffusion.

Edit 3: and now I see that actually, oceandrift handles the sea floor exactly as I would like, dealing with diffusion by reflection, and sinking by deactivation.

[nordam]

And on an even closer look, I see that interact_with_seafloor is also called twice in the main loop in basemodel/__init__.py. This makes it a bit tricky to get the desired behaviour out of the seafloor interaction in OceanDrift. After debugging, I found that when interact_with_seafloor is called from the main loop in basemodel/__init__.py, sometimes particles with positive (upwards) terminal velocity will for some reason have ended up below the sea floor (possibly due to sideways advection?), and will then be deactivated. I'm trying to model a case where only particles with negative velocities will be deactivated at the sea floor.

A behaviour matching the boundary condition I'm trying to implement would be to lift to the floor any particles that end up under the floor for any reason other than sinking. However, interact_with_seafloor reads only one global setting, and does the same thing both when called in the main loop, and when called as part of the vertical transport routine from the OceanDrift class.

Some more fine grained control would be convenient. Maybe defining a few different seafloor functions, such as interact_with_seafloor_after_advection and interact_with_seafloor_after_buoyancy, etc.?