Merge pull request #1641 from knutfrode/dev

Implemented correction activated by config drift:vertical_velocity_co…

Author: knutfrode

examples/example_vertical_velocity_correction.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+"""
+Vertical advection correction
+=============================
+"""
+
+#%%
+# In OpenDrift, element depth (z) is relative to actual surface, and not mean sea level.
+# Vertical velocity from ocean models contain a contribution which is due to
+# change in sea surface elevation, however, this should in OpenDrift be subtracted
+# due to the choice of defining z relative to actual surface elevation
+# This correction is activated with config setting 'drift:vertical_advection_correction'
+# The effect is illustrated with the simulations below.
+
+
+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
+
+time = datetime.now()
+reader_tidal = reader_oscillating.Reader('sea_surface_height', amplitude=-1,
+                                         period=timedelta(hours=6), zero_time=time)
+lat=59.8113; lon=10.5517  # Oslo fjord
+z = np.arange(0, -60, -5)  # Seeding one particle every 5 meter from surface to 60m depth
+
+# Without stetching
+o = OceanDrift(loglevel=0)
+o.add_reader(reader_tidal)
+o.add_readers_from_list(['https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be'])
+o.set_config('drift:water_column_stretching', False)
+o.set_config('drift:vertical_advection', True)
+o.seed_elements(lon=lon, lat=lat, time=time, z=z, number=len(z))
+o.run(duration=timedelta(hours=24), time_step=1800)
+
+# With stetching
+o2 = OceanDrift(loglevel=0)
+o2.add_reader(reader_tidal)
+o2.add_readers_from_list(['https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be'])
+o2.set_config('drift:water_column_stretching', True)
+o2.set_config('drift:vertical_advection', True)
+o2.seed_elements(lon=lon, lat=lat, time=time, z=z, number=len(z))
+o2.run(duration=timedelta(hours=24), time_step=1800)
+
+# With w correction
+o3 = OceanDrift(loglevel=0)
+o3.add_reader(reader_tidal)
+o3.add_readers_from_list(['https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be'])
+o3.set_config('drift:water_column_stretching', False)
+o3.set_config('drift:vertical_advection', True)
+o3.set_config('drift:vertical_advection_correction', True)
+o3.seed_elements(lon=lon, lat=lat, time=time, z=z, number=len(z))
+o3.run(duration=timedelta(hours=24), time_step=1800)
+
+
+plt.subplot(3,1,1)
+o.result.z.plot.line(x='time', add_legend=False, color='k')
+o2.result.z.plot.line(x='time', add_legend=False, color='r')
+o3.result.z.plot.line(x='time', add_legend=False, color='g')
+plt.plot([], [], color='k', label='No stretching or correction')
+plt.plot([], [], color='r', label='Column stretching')
+plt.plot([], [], color='g', label='Vertical velocity correction')
+plt.legend()
+plt.subplot(3,1,2)
+o.result.sea_floor_depth_below_sea_level.plot.line(x='time', add_legend=False, color='k')
+plt.subplot(3,1,3)
+o.result.upward_sea_water_velocity.plot.line(x='time', add_legend=False, color='k')
+plt.show()

examples/example_water_column_stretching.py

@@ -45,6 +45,7 @@
 # This lifting is in this case unphysical.
 
 #%%
+# We then make a simulation wih correction for the stretching/contraction of the water column.
 o = OceanDrift(loglevel=20)
 o.add_reader(reader_tidal)
 o.set_config('drift:water_column_stretching', True)
@@ -54,6 +55,7 @@
 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

@@ -124,6 +124,9 @@ def __init__(self, *args, **kwargs):
             '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_advection_correction': {'type': 'bool', 'default': False, 'description':
+                'The vertical velocity from ocean model is partly due to sea level changes. But since element depth in OpenDrift (elements.z) is relative to actual surface, this part should be subtracted.',
+                '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':
@@ -301,6 +304,7 @@ def water_column_stretching(self):
         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)
+        self.elements.z = self.elements.z.data
 
     def vertical_advection(self):
         """Move particles vertically according to vertical ocean current
@@ -321,8 +325,19 @@ def vertical_advection(self):
 
         if len(applicable) > 0:
             w = self.environment.upward_sea_water_velocity[applicable]
+
+            if self.get_config('drift:vertical_advection_correction', None) is True:
+                logger.debug('Subtracting motion due to elevation change from vertical water velocity')
+                delta_zeta = self.environment.sea_surface_height[applicable] - self.environment_previous.sea_surface_height[applicable]
+                w_surface = delta_zeta / self.time_step.total_seconds()
+                total_depth = self.environment.sea_surface_height[applicable] + \
+                     self.environment.sea_floor_depth_below_sea_level[applicable]
+                w_elevation = w_surface * (self.elements.z[applicable] + total_depth) / total_depth
+                w = w - w_elevation
+
             self.elements.z[applicable] = np.minimum(0,
-                self.elements.z[applicable] + self.elements.moving[applicable] * w * self.time_step.total_seconds())
+                self.elements.z[applicable] + self.elements.moving[applicable] * w *
+                self.time_step.total_seconds())
 
     def vertical_buoyancy(self):
         """Move particles vertically according to their buoyancy"""