diff --git a/examples/example_biodegradation.py b/examples/example_biodegradation.py new file mode 100755 index 000000000..b163f72a8 --- /dev/null +++ b/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 diff --git a/opendrift/models/openoil/openoil.py b/opendrift/models/openoil/openoil.py index 836c3c15b..0872ca9bc 100644 --- a/opendrift/models/openoil/openoil.py +++ b/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: