Merge pull request #1247 from knutfrode/dev

Added experimental alternative oil biodegradation with decay_rate, co…

examples/example_biodegradation.py

@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+"""
+Biodegradation of oil
+======================
+"""
+
+import numpy as np
+from datetime import datetime, timedelta
+from opendrift.readers import reader_netCDF_CF_generic
+from opendrift.models.openoil import OpenOil
+
+o = OpenOil(loglevel=0)  # Set loglevel to 0 for debug information
+time = datetime.now()
+
+# No motion is needed for this test
+o.set_config('environment:constant', {k: 0 for k in
+             ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']})
+o.set_config('drift', {'current_uncertainty': 0, 'wind_uncertainty': 0})
+
+#%%
+# Seeding some particles
+o.set_config('drift:vertical_mixing', False)  # Keeping oil at fixed depths
+o.set_config('processes:biodegradation', True)
+o.set_config('biodegradation:method', 'decay_rate')
+
+#%%
+# Fast decay for droplets, and slow decay for slick 
+decay = {'biodegradation_decay_rate_slick': np.log(2)/timedelta(days=100).total_seconds(),
+         'biodegradation_decay_rate_droplet': np.log(2)/timedelta(days=3).total_seconds(),
+         'oil_type': 'GENERIC MEDIUM CRUDE'}
+
+#%%
+# Seed 5000 oil elements at surface, and 5000 elements at 100m depth
+o.seed_elements(lon=4, lat=60.0, z=0, number=5000, time=datetime.now(), **decay)
+o.seed_elements(lon=4, lat=60.0, z=-100, number=5000, time=datetime.now(), **decay)
+
+#%%
+# Running model
+o.run(duration=timedelta(hours=24), time_step=900) 
+
+#%%
+# Print and plot results
+o.plot_oil_budget(show_watercontent_and_viscosity=False, show_wind_and_current=False)
+
+#%%
+# .. image:: /gallery/animations/example_biodegradation_0.gif

opendrift/models/openoil/openoil.py

@@ -167,13 +167,17 @@ class Oil(Lagrangian3DArray):
         }),
         ('mass_biodegraded', {
             'dtype': np.float32,
-            'units': 'kg',
+            'units': 'fraction',
             'seed': False,
             'default': 0
         }),
+        ('biodegradation_decay_rate_droplet',  # TODO: should have valid_min and valid_max for element properties
+                {'dtype': np.float32, 'units': 'kg', 'seed': False, 'default': np.log(2)/(3600*24)}),  # 1 day default
+        ('biodegradation_decay_rate_slick',
+                {'dtype': np.float32, 'units': 'kg', 'seed': False, 'default': np.log(2)/(3600*24*3)}),  # 3 days default
         ('fraction_evaporated', {
             'dtype': np.float32,
-            'units': '%',
+            'units': '%',  # TODO: should be fraction and not percent
             'seed': False,
             'default': 0
         }),
@@ -430,6 +434,12 @@ def __init__(self, weathering_model='noaa', *args, **kwargs):
                 'description': 'Oil mass is biodegraded (eaten by bacteria).',
                 'level': CONFIG_LEVEL_BASIC
             },
+            'biodegradation:method': {
+                'type': 'enum',
+                'enum': ['Adcroft', 'decay_rate'],
+                'default': 'Adcroft', 'level': CONFIG_LEVEL_ADVANCED,
+                'description': 'Alogorithm to be used for biodegradation of oil'
+            },
             'processes:update_oilfilm_thickness': {
                 'type': 'bool',
                 'default': False,
@@ -538,32 +548,57 @@ def update_surface_oilfilm_thickness(self):
 
     def biodegradation(self):
         if self.get_config('processes:biodegradation') is True:
-            '''
-            Oil biodegradation function based on the article:
-            Adcroft et al. (2010), Simulations of underwater plumes of
-            dissolved oil in the Gulf of Mexico.
-            '''
-            logger.debug('Calculating: biodegradation')
-
-            swt = self.environment.sea_water_temperature.copy()
-            swt[swt > 100] -= 273.15  # K to C
-            age0 = self.time_step.total_seconds() / (3600 * 24)
-
-            # Decay rate in days (temperature in Celsius)
-            tau = (12) * (3**((20 - swt) / 10))
-
-            fraction_biodegraded = (1 - np.exp(-age0 / tau))
-            biodegraded_now = self.elements.mass_oil * fraction_biodegraded
-
-            self.elements.mass_biodegraded = \
-                self.elements.mass_biodegraded + biodegraded_now
-            self.elements.mass_oil = \
-                self.elements.mass_oil - biodegraded_now
-            if self.oil_weathering_model == 'noaa':
-                self.noaa_mass_balance['mass_components'][self.elements.ID - 1, :] = \
-                self.noaa_mass_balance['mass_components'][self.elements.ID - 1, :]*(1-fraction_biodegraded[:, np.newaxis])
-            else:
-                pass
+            method = self.get_config('biodegradation:method')
+            logger.debug(f'Calculating: biodegradation ({method})')
+            if method == 'Adcroft':
+                self.biodegradation_adcroft()
+            elif method == 'decay_rate':
+                self.biodegradation_decay_rate()
+
+    def biodegradation_decay_rate(self):
+        '''Oil biodegradation with exponential decay'''
+
+        surface = np.where(self.elements.z == 0)[0]  # of active elements
+        age0 = self.time_step.total_seconds()
+
+        half_time = np.log(2) / self.elements.biodegradation_decay_rate_droplet
+        half_time[surface] = np.log(2) / self.elements.biodegradation_decay_rate_slick[surface]
+
+        fraction_biodegraded = (1 - np.exp(-age0 / half_time))
+        biodegraded_now = self.elements.mass_oil * fraction_biodegraded
+
+        self.elements.mass_biodegraded = \
+            self.elements.mass_biodegraded + biodegraded_now
+        self.elements.mass_oil = \
+            self.elements.mass_oil - biodegraded_now
+        if self.oil_weathering_model == 'noaa':
+            self.noaa_mass_balance['mass_components'][self.elements.ID - 1, :] = \
+            self.noaa_mass_balance['mass_components'][self.elements.ID - 1, :]*(1-fraction_biodegraded[:, np.newaxis])
+
+    def biodegradation_adcroft(self):
+        '''
+        Oil biodegradation function based on the article:
+        Adcroft et al. (2010), Simulations of underwater plumes of
+        dissolved oil in the Gulf of Mexico.
+        '''
+
+        swt = self.environment.sea_water_temperature.copy()
+        swt[swt > 100] -= 273.15  # K to C
+        age0 = self.time_step.total_seconds() / (3600 * 24)
+
+        # Decay rate in days (temperature in Celsius)
+        tau = (12) * (3**((20 - swt) / 10))
+
+        fraction_biodegraded = (1 - np.exp(-age0 / tau))
+        biodegraded_now = self.elements.mass_oil * fraction_biodegraded
+
+        self.elements.mass_biodegraded = \
+            self.elements.mass_biodegraded + biodegraded_now
+        self.elements.mass_oil = \
+            self.elements.mass_oil - biodegraded_now
+        if self.oil_weathering_model == 'noaa':
+            self.noaa_mass_balance['mass_components'][self.elements.ID - 1, :] = \
+            self.noaa_mass_balance['mass_components'][self.elements.ID - 1, :]*(1-fraction_biodegraded[:, np.newaxis])
 
     def disperse(self):
         if self.get_config('processes:dispersion') is True: