Merge pull request #1640 from knutfrode/dev

[run-ex] New config setting 'drift:water_column_stretching' for Ocean…

Author: knutfrode

examples/example_water_column_stretching.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+"""
+Water column stretching
+========================
+"""
+
+from datetime import datetime, timedelta
+import numpy as np
+import matplotlib.pyplot as plt
+from opendrift.readers import reader_oscillating
+from opendrift.models.oceandrift import OceanDrift
+
+
+#%%
+# In OpenDrift, the vertical position of elements ("z") is defined as relative to actual surface,
+# and not an absolute reference level (e.g. mean sea surface height).
+# Thus if sea surface elevation changes with time (e.g. tides),
+# we need to add a "correction / perturbation" to z, otherwise elements at/near seafloor
+# will be lifted if surface elevation increases and z (relative to surface) remains unchanged.
+# This correction is presently only implemented for OceanDrift, and must be switched on
+# with config setting "drift:water:column_stretching"
+
+# To illustrate, we add a reader with oscillating sea surface elevation (tidal)
+# with amplitude of 1m and peroid of 6 hours
+time = datetime.now()
+reader_tidal = reader_oscillating.Reader('sea_surface_height', amplitude=-1,
+                                         period=timedelta(hours=6), zero_time=time)
+
+#%%
+# First an illustration withouth this correction.
+o = OceanDrift(loglevel=20)
+o.add_reader(reader_tidal)
+o.set_config('drift:water_column_stretching', False)
+o.set_config('environment:constant:sea_floor_depth_below_sea_level', 10)
+z = np.arange(0, -11, -1)  # Seeding one particle every meter from surface to 10m depth
+o.seed_elements(lon=0, lat=0, time=time, z=z, number=11)
+o.run(duration=timedelta(hours=24), time_step=1800)
+o.result.z.plot.line(x='time', add_legend=False)
+plt.show()
+
+#%%
+# We see that the particles remain at their initial depths (since we have no vertical advection or mixing),
+# except for the element starting at seafloor, which is lifted up when sea level rises,
+# since the config setting `drift:seafloor_action` is `lift_to_seafloor` by default.
+# This lifting is in this case unphysical.
+
+#%%
+o = OceanDrift(loglevel=20)
+o.add_reader(reader_tidal)
+o.set_config('drift:water_column_stretching', True)
+o.set_config('environment:constant:sea_floor_depth_below_sea_level', 10)
+o.seed_elements(lon=0, lat=0, time=time, z=z, number=11)
+o.run(duration=timedelta(hours=24), time_step=1800)
+o.result.z.plot.line(x='time', add_legend=False)
+plt.show()
+
+# Here we see that element depth (z, relative to surface) is changed so that
+# elements at surface and seafloor remain at resp surface (z=0) and
+# seafloor (z = sea_floor_depth + sea_surface_elevation)

opendrift/models/oceandrift.py

@@ -121,6 +121,9 @@ def __init__(self, *args, **kwargs):
             'drift:vertical_advection_at_surface': {'type': 'bool', 'default': False, 'description':
                 'If vertical advection is activated, surface elements (z=0) can only be advected (downwards) if this setting it True.',
                 'level': CONFIG_LEVEL_ADVANCED},
+            'drift:water_column_stretching': {'type': 'bool', 'default': False, 'description':
+                'If sea surface elevation changes, vertical particle position (element property "z" is relative to surface, and not absolute zero-level) is adjusted so that elements at surface and seafloor remains at resp surface and seafloor.',
+                'level': CONFIG_LEVEL_ADVANCED},
             'drift:vertical_mixing': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_BASIC,
                 'description': 'Activate vertical mixing scheme with inner loop'},
             'drift:vertical_mixing_at_surface': {'type': 'bool', 'default': False, 'description':
@@ -292,8 +295,12 @@ def prepare_run(self):
     def water_column_stretching(self):
         """If sea_water_level changes, adjust z for continuity"""
 
-        return  # To be implemented
+        if self.get_config('drift:water_column_stretching') is False:
+            return
+
         delta_zeta = self.environment.sea_surface_height - self.environment_previous.sea_surface_height
+        logger.info('Compensating for change in surface elevation')
+        self.elements.z = self.elements.z + delta_zeta*(self.elements.z/self.environment.sea_floor_depth_below_sea_level)
 
     def vertical_advection(self):
         """Move particles vertically according to vertical ocean current