# This file is part of OpenDrift.
#
# OpenDrift is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2
#
# OpenDrift is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with OpenDrift. If not, see <https://www.gnu.org/licenses/>.
#
# Copyright 2015, Knut-Frode Dagestad, MET Norway
# Copyright 2021, Gaute Hope, MET Norway
"""
OpenOil is a 3D oil drift module bundled within the OpenDrift framework.
The oil weathering calculations is based on the NOAA-ERR-ERD OilLibrary
package, which is installed as a dependency. The code for evaporation and
emulsification in OpenOil is borrowed from the NOAA PyGnome code, and adapted
to the OpenDrift architecture.
Example of ship leaking oil along the coast of Northern Norway
##############################################################
.. image:: https://dl.dropboxusercontent.com/s/ty6dmqf0oohewky/oilspill_tromsoe.gif?dl=0
Simulation of Deepwater Horizon (Macondo) accident, initiated from satellite images
###################################################################################
.. image:: https://dl.dropboxusercontent.com/s/ghi7crtmwpyjgto/macondo_simulation.gif?dl=0
Satellite images provided by Prof. Chuanmin Hu, and ocean model output provided by Prof. Vassiliki Kourafalou
Example oil budget for a simulation
###################################
.. image:: https://dl.dropboxusercontent.com/s/pb0h6tlev9pnoh3/oil_budget_draugen.png?dl=0
Oil properties affecting the drift
***********************************
The vertical (and thus indirectly also the horizontal) motion of oil (droplets) is affected by oil density and droplet diameters.
When using the NOAA oil weathering model (``o = OpenOil(weathering_model='noaa')``), which is the default, the density is obtained from the NOAA database according to the oiltype selected when seeding. This value can not be overridden by the user, and it will also change during the simulation due to oil weathering processes (evaporation and emulsification).
The droplet diameter may be given explicitly when seeding, e.g.:
.. testcode::
o = OpenOil()
o.set_config('environment:constant:x_wind', 0)
o.set_config('environment:constant:y_wind', 0)
o.set_config('environment:constant:x_sea_water_velocity', 0)
o.set_config('environment:constant:y_sea_water_velocity', 0)
o.seed_elements(4, 60, number=100, time=datetime.now(), diameter=1e-5)
In this case, the diameter will not change during the simulation, which is useful e.g. for sensitivity tests. The same diameter will be used for all elements for this example, but an array of the same length as the number of elements may also be provided.
If a constant droplet diameter is not given by the user, it will be chosen randomly within given config limits for a subsea spill ('blowout'), and modified after any later wave breaking event. Oil droplets seeded under sea surface (z<0) will be assigned initial diameters between the following limits, typical for a subsea blowout (Johansen, 2000)::
.. code::
o.set_config('seed:droplet_diameter_min_subsea', 0.0005) # 0.5 mm
o.set_config('seed:droplet_diameter_max_subsea', 0.005) # 5 mm
Alternatively, the user can specify normal or lognormal initial subsea droplet size distributions, which are later modified by wave breaking events. In these cases the user must specify the mean and standard deviation of the distribution::
.. code::
o.set_config('seed:droplet_size_distribution','lognormal')
o.set_config('seed:droplet_diameter_mu',0.001) # 1 mm
o.set_config('seed:droplet_diameter_sigma',0.0008) # 0.8 mm
Note that these config settings must be adjusted before the seeding call.
After each wave breaking event, a new droplet diameter will be chosen based on the config setting for droplet size distribution.
"""
from io import open
import numpy as np
from datetime import datetime
import pyproj
import matplotlib.pyplot as plt
import logging
import json
from importlib import resources
logger = logging.getLogger(__name__)
from opendrift.models.oceandrift import OceanDrift, Lagrangian3DArray
from . import noaa_oil_weathering as noaa
from . import adios
from adios_db.computation.physical_properties import KinematicViscosity, Density
from adios_db.computation import gnome_oil
from opendrift.models.physics_methods import oil_wave_entrainment_rate_li2017
from opendrift.config import CONFIG_LEVEL_ESSENTIAL, CONFIG_LEVEL_BASIC, CONFIG_LEVEL_ADVANCED
# Defining the oil element properties
[docs]
class Oil(Lagrangian3DArray):
"""Extending LagrangianArray with variables relevant for oil particles."""
variables = Lagrangian3DArray.add_variables([
('mass_oil', {
'dtype': np.float32,
'units': 'kg',
'seed': False,
'default': 1
}),
(
'viscosity',
{
'dtype': np.float32,
'units': 'm2/s',
'seed': False, # Taken from NOAA database
'description': 'Kinematic viscosity of oil (emulsion)',
'default': 0.005
}),
(
'density',
{
'dtype': np.float32,
'units': 'kg/m^3',
'seed': False, # Taken from NOAA database
'default': 880
}),
(
'wind_drift_factor',
{
'dtype':
np.float32, # TODO: inherit from
'units':
'%', # OceanDrift
'description':
'Elements at the ocean surface are moved by '
'this fraction of the wind vector, in addition to '
'currents and Stokes drift',
'default':
0.03
}),
('bulltime', {
'dtype': np.float32,
'units': 's',
'seed': False,
'default': 0
}),
('interfacial_area', {
'dtype': np.float32,
'units': 'm2',
'seed': False,
'default': 0
}),
('mass_dispersed', {
'dtype': np.float32,
'units': 'kg',
'seed': False,
'default': 0
}),
('mass_evaporated', {
'dtype': np.float32,
'units': 'kg',
'seed': False,
'default': 0
}),
('mass_biodegraded', {
'dtype': np.float32,
'units': 'kg',
'seed': False,
'default': 0
}),
('biodegradation_half_time_droplet', # Note: if unit is "days", Xarray tries to decode as datetime with error
{'dtype': np.float32, 'units': 'Days', 'seed': False, 'default': 1,
'description': 'Biodegradation half time in days for submerged oil droplets'}),
('biodegradation_half_time_slick',
{'dtype': np.float32, 'units': 'Days', 'seed': False, 'default': 3,
'description': 'Biodegradation half time in days for surface oil slick'}),
('fraction_evaporated', {
'dtype': np.float32,
'units': '%', # TODO: should be fraction and not percent
'seed': False,
'default': 0
}),
('water_fraction', {
'dtype': np.float32,
'units': '%',
'seed': False,
'default': 0
}),
('oil_film_thickness', {
'dtype': np.float32,
'units': 'm',
'default': 0.001
}),
(
'diameter',
{
'dtype': np.float32, # Particle diameter
'units': 'm',
'seed': False,
'default': 0.
})
])
[docs]
class OpenOil(OceanDrift):
"""Open source oil trajectory model based on the OpenDrift framework.
Developed at MET Norway based on oil weathering parameterisations
found in open/published litterature.
Under construction.
"""
ElementType = Oil
required_variables = {
'x_sea_water_velocity': {
'fallback': None
},
'y_sea_water_velocity': {
'fallback': None
},
'x_wind': {
'fallback': None
},
'y_wind': {
'fallback': None
},
'sea_surface_height': {'fallback': 0},
'upward_sea_water_velocity': {
'fallback': 0,
'important': False
},
'sea_surface_wave_significant_height': {
'fallback': 0,
'important': False
},
'sea_surface_wave_stokes_drift_x_velocity': {
'fallback': 0,
'important': False
},
'sea_surface_wave_stokes_drift_y_velocity': {
'fallback': 0,
'important': False
},
'sea_surface_wave_period_at_variance_spectral_density_maximum': {
'fallback': 0,
'important': False
},
'sea_surface_wave_mean_period_from_variance_spectral_density_second_frequency_moment':
{
'fallback': 0,
'important': False
},
'sea_ice_area_fraction': {
'fallback': 0,
'important': False
},
'sea_ice_x_velocity': {
'fallback': 0,
'important': False
},
'sea_ice_y_velocity': {
'fallback': 0,
'important': False
},
'sea_water_temperature': {
'fallback': 10,
'profiles': True
},
'sea_water_salinity': {
'fallback': 34,
'profiles': True
},
'sea_floor_depth_below_sea_level': {
'fallback': 10000
},
'ocean_vertical_diffusivity': {
'fallback': 0.02,
'important': False,
'profiles': True
},
'land_binary_mask': {
'fallback': None
},
'ocean_mixed_layer_thickness': {
'fallback': 50,
'important': False
},
}
# Default colors for plotting
status_colors = {
'initial': 'green',
'active': 'blue',
'missing_data': 'gray',
'stranded': 'red',
'evaporated': 'yellow',
'dispersed': 'magenta'
}
duplicate_oils = ['ALVHEIM BLEND, STATOIL', 'DRAUGEN, STATOIL',
'EKOFISK BLEND 2000', 'EKOFISK BLEND, STATOIL',
'EKOFISK, CITGO', 'EKOFISK, EXXON', 'EKOFISK, PHILLIPS',
'EKOFISK, STATOIL', 'ELDFISK', 'ELDFISK B',
'GLITNE, STATOIL', 'GOLIAT BLEND, STATOIL',
'GRANE BLEND, STATOIL', 'GUDRUN BLEND, STATOIL',
'GULLFAKS A, STATOIL', 'GULLFAKS C, STATOIL',
'GULLFAKS, SHELL OIL', 'GULLFAKS SOR',
'GULLFAKS, STATOIL', 'HEIDRUN, STATOIL',
'NJORD, STATOIL', 'NORNE, STATOIL',
'OSEBERG BLEND, STATOIL', 'OSEBERG EXXON',
'OSEBERG, PHILLIPS', 'OSEBERG, SHELL OIL',
'SLEIPNER CONDENSATE, STATOIL',
'STATFJORD BLEND, STATOIL', 'VARG, STATOIL']
def __init__(self, weathering_model='noaa', *args, **kwargs):
self.oil_weathering_model = weathering_model
if self.oil_weathering_model == 'noaa': # Currently the only option
self.oiltypes = adios.get_oil_names(
location=kwargs.get('location', None))
# Update config with oiltypes
self.oiltypes.extend(adios.oil_name_alias.keys())
# Sort alphabetically, but put GENERIC oils first
generic_oiltypes = [o for o in self.oiltypes if o[0:7] == 'GENERIC']
other_oiltypes = [o for o in self.oiltypes if o[0:7] != 'GENERIC']
self.oiltypes = sorted([o for o in generic_oiltypes]) + sorted([o for o in other_oiltypes])
self.oiltypes = [ot for ot in self.oiltypes if ot not in self.duplicate_oils]
# For Norwegian oils, max water fraction from Sintef is overriding NOAA value
self.max_water_fraction = None
else:
raise ValueError('Weathering model unknown: ' + weathering_model)
kwargs.pop('location', None)
# Calling general constructor of parent class
super(OpenOil, self).__init__(*args, **kwargs)
self._add_config({
'seed:m3_per_hour': {
'type': 'float',
'default': 1,
'min': 0,
'max': 1e10,
'units': 'm3 per hour',
'description':
'The amount (volume) of oil released per hour (or total amount if release is instantaneous)',
'level': CONFIG_LEVEL_ESSENTIAL
},
'seed:droplet_size_distribution': {
'type':
'enum',
'enum': ['uniform', 'normal', 'lognormal'],
'default':
'uniform',
'level':
CONFIG_LEVEL_ADVANCED,
'description':
'Droplet size distribution used for subsea release.'
},
'seed:droplet_diameter_mu': {
'type': 'float',
'default': 0.001,
'min': 1e-8,
'max': 1,
'units': 'meters',
'description':
'The mean diameter of oil droplet for a subsea release, used in normal/lognormal distributions.',
'level': CONFIG_LEVEL_BASIC
},
'seed:droplet_diameter_sigma': {
'type': 'float',
'default': 0.0005,
'min': 1e-8,
'max': 1,
'units': 'meters',
'description':
'The standard deviation in diameter of oil droplet for a subsea release, used in normal/lognormal distributions.',
'level': CONFIG_LEVEL_BASIC
},
'seed:droplet_diameter_min_subsea': {
'type': 'float',
'default': 0.0005,
'min': 1e-8,
'max': 1,
'units': 'meters',
'description':
'The minimum diameter of oil droplet for a subsea release, used in unifrom distribution.',
'level': CONFIG_LEVEL_BASIC
},
'seed:droplet_diameter_max_subsea': {
'type': 'float',
'default': 0.005,
'min': 1e-8,
'max': 1,
'units': 'meters',
'description':
'The maximum diameter of oil droplet for a subsea release, used in uniform distribution.',
'level': CONFIG_LEVEL_BASIC
},
'processes:dispersion': {
'type': 'bool',
'default': True,
'description':
'Oil is removed from simulation (dispersed), if entrained as very small droplets.',
'level': CONFIG_LEVEL_BASIC
},
'processes:evaporation': {
'type': 'bool',
'default': True,
'description': 'Surface oil is evaporated.',
'level': CONFIG_LEVEL_BASIC
},
'processes:emulsification': {
'type': 'bool',
'default': True,
'description':
'Surface oil is emulsified, i.e. water droplets are mixed into oil due to wave mixing, with resulting increas of viscosity.',
'level': CONFIG_LEVEL_BASIC
},
'processes:biodegradation': {
'type': 'bool',
'default': False,
'description': 'Oil mass is biodegraded (eaten by bacteria).',
'level': CONFIG_LEVEL_BASIC
},
'biodegradation:method': {
'type': 'enum',
'enum': ['Adcroft', 'half_time'],
'default': 'Adcroft', 'level': CONFIG_LEVEL_ADVANCED,
'description': 'Alogorithm to be used for biodegradation of oil'
},
'processes:update_oilfilm_thickness': {
'type': 'bool',
'default': False,
'description':
'Oil film thickness is calculated at each time step. The alternative is that oil film thickness is kept constant with value provided at seeding.',
'level': CONFIG_LEVEL_ADVANCED
},
'wave_entrainment:droplet_size_distribution': {
'type':
'enum',
'enum': ['Johansen et al. (2015)', 'Li et al. (2017)'],
'default':
'Johansen et al. (2015)',
'level':
CONFIG_LEVEL_ADVANCED,
'description':
'Algorithm to be used for calculating oil droplet size spectrum after entrainment by breaking waves.'
},
'wave_entrainment:entrainment_rate': {
'type':
'enum',
'enum': ['Li et al. (2017)'],
'default':
'Li et al. (2017)',
'level':
CONFIG_LEVEL_ADVANCED,
'description':
'Algorithm to be used for calculating the entrainment rate of oil due to wave breaking.'
},
'seed:oil_type': {
'type':
'enum',
'enum':
self.oiltypes,
'default':
self.oiltypes[0],
'level':
CONFIG_LEVEL_ESSENTIAL,
'description':
'Oil type to be used for the simulation, from the NOAA ADIOS database.'
},
})
self._set_config_default('drift:vertical_advection', False)
self._set_config_default('drift:vertical_mixing', True)
self._set_config_default('drift:current_uncertainty', 0.05)
self._set_config_default('drift:wind_uncertainty', 0.5)
self._set_config_default('drift:max_speed', 1.3)
[docs]
def update_surface_oilfilm_thickness(self):
'''The mass of oil is summed within a grid of 100x100
cells covering the oil at a given time. Each oil particle
within each cell is given a film thickness as the amount of
oil divided by the cell area.
'''
from scipy.stats import binned_statistic_2d
surface = np.where(self.elements.z == 0)[0]
if len(surface) == 0:
logger.debug('No oil at surface, no film thickness to update')
return
logger.debug(
'Updating oil film thickness for %s of %s elements at surface' %
(len(surface), self.num_elements_active()))
meanlon = self.elements.lon[surface].mean()
meanlat = self.elements.lat[surface].mean()
# Using stereographic coordinates to get regular X and Y
psproj = pyproj.Proj('+proj=stere +lat_0=%s +lat_ts=%s +lon_0=%s' %
(meanlat, meanlat, meanlon))
X, Y = psproj(self.elements.lon[surface], self.elements.lat[surface])
mass_bin, x_edge, y_edge, binnumber = binned_statistic_2d(
X,
Y,
self.elements.mass_oil[surface],
expand_binnumbers=True,
statistic='sum',
bins=100)
bin_area = (x_edge[1] - x_edge[0]) * (y_edge[1] - y_edge[0])
oil_density = 1000 # ok approximation here
film_thickness = (mass_bin / oil_density) / bin_area
# Postulating min and max film thickness
max_thickness = 0.01 # 1 cm
min_thickness = 1e-9 # 1 nanometer
if film_thickness.max() > max_thickness:
logger.debug(
'Warning: decreasing thickness to %sm for %s of %s bins' %
(max_thickness, np.sum(film_thickness > max_thickness),
film_thickness.size))
film_thickness[film_thickness > max_thickness] = max_thickness
num_too_thin = np.sum((film_thickness < min_thickness)
& (film_thickness > 0))
if num_too_thin > 0:
logger.debug(
'Warning: increasing thickness to %sm for %s of %s bins' %
(min_thickness, num_too_thin, film_thickness.size))
film_thickness[film_thickness < min_thickness] = min_thickness
# https://github.com/scipy/scipy/issues/7010
binnumber = binnumber - 1
bx = binnumber[0, :]
by = binnumber[1, :]
# Update thickness
self.elements.oil_film_thickness[
surface] = self.elements.oil_film_thickness[surface] * np.nan
self.elements.oil_film_thickness[surface] = \
film_thickness[bx, by]
[docs]
def biodegradation(self):
if self.get_config('processes:biodegradation') is True:
method = self.get_config('biodegradation:method')
logger.debug(f'Calculating: biodegradation ({method})')
if method == 'Adcroft':
self.biodegradation_adcroft()
elif method == 'half_time':
self.biodegradation_half_time()
[docs]
def biodegradation_half_time(self):
'''Oil biodegradation with exponential decay'''
surface = np.where(self.elements.z == 0)[0] # of active elements
age0 = self.time_step.total_seconds()/(3600*24) # days
half_time = self.elements.biodegradation_half_time_droplet.copy()
half_time[surface] = self.elements.biodegradation_half_time_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])
[docs]
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])
[docs]
def disperse(self):
if self.get_config('processes:dispersion') is True:
logger.debug(' Calculating: dispersion')
windspeed = np.sqrt(self.environment.x_wind**2 +
self.environment.y_wind**2)
# From NOAA PyGnome model:
# https://github.com/NOAA-ORR-ERD/PyGnome/
v_entrain = 3.9E-8
sea_water_density = 1028
fraction_breaking_waves = 0.02
wave_significant_height = \
self.significant_wave_height()
wave_significant_height[wave_significant_height == 0] = \
0.0246*windspeed[wave_significant_height == 0]**2
dissipation_wave_energy = \
(0.0034*sea_water_density*9.81*wave_significant_height**2)
c_disp = np.power(dissipation_wave_energy, 0.57) * \
fraction_breaking_waves
# Roy's constant
C_Roy = 2400.0 * np.exp(-73.682 * np.sqrt(
self.elements.viscosity / self.elements.density))
q_disp = C_Roy * c_disp * v_entrain / self.elements.density
oil_mass_loss = (q_disp * self.time_step.total_seconds() *
self.elements.density) * self.elements.mass_oil
self.elements.mass_oil -= oil_mass_loss
self.elements.mass_dispersed += oil_mass_loss
#self.elements.z = \
# self.environment.sea_surface_wave_significant_height
## Marks R. (1987), Marine aerosols and whitecaps in the
## North Atlantic and Greenland sea regions
## Deutsche Hydrografische Zeitschrift, Vol 40, Issue 2 , pp 71-79
#whitecap_coverage = (2.54E-4)*np.power(windspeed, 3.58)
## Martinsen et al. (1994), The operational
## oil drift system at DNMI
## DNMI Technical report No 125, 51
#wave_period = 3.85*np.sqrt(
# self.environment.sea_surface_wave_significant_height) # sec
#time_between_breaking_events = 3.85/whitecap_coverage # TBC
#rho_w = 1025 # kg/m3
#wave_period[wave_period == 0] = 5 # NB: temporal fix for no waves
#dissipation_energy = 0.0034*rho_w*9.81*(wave_period**2)
#dsize_coeff = 2100 # Wettre et al., Appendix A
#c_oil = dsize_coeff*np.power(self.elements.viscosity, -0.4)
## Random numbers between 0 and 1:
#p = np.random.rand(self.num_elements_active(), )
#oil_per_unit_surface = 1
#droplet_size = np.power((p*oil_per_unit_surface)/(c_oil *
# np.power(dissipation_energy, 0.57)),
# 1/1.17)
#self.deactivate_elements(droplet_size < 5E-7, reason='dispersed')
[docs]
def oil_weathering(self):
if self.time_step.days < 0:
logger.debug('Skipping oil weathering for backwards run')
return
self.timer_start('main loop:updating elements:oil weathering')
if self.oil_weathering_model == 'noaa':
self.oil_weathering_noaa()
self.timer_end('main loop:updating elements:oil weathering')
[docs]
def prepare_run(self):
if self.oil_weathering_model == 'noaa':
self.noaa_mass_balance = {}
# Populate with seeded mass spread on oiltype.mass_fraction
mass_oil = np.atleast_1d(self.elements_scheduled.mass_oil)
if len(mass_oil) == 1:
mass_oil = mass_oil * np.ones(self.num_elements_total())
self.noaa_mass_balance['mass_components'] = \
np.asarray(self.oiltype.mass_fraction)*(mass_oil.reshape(
(self.num_elements_total(), -1)))
self.noaa_mass_balance['mass_evaporated'] = \
self.noaa_mass_balance['mass_components']*0
self.oil_water_interfacial_tension = \
self.oiltype.oil_water_surface_tension()
logger.info('Oil-water surface tension is %f Nm' %
self.oil_water_interfacial_tension)
try:
max_water_fractions = json.loads(
resources.read_text('opendrift.models.openoil.adios', 'max_water_fraction.json'))
if self.oil_name in max_water_fractions:
self.max_water_fraction = max_water_fractions[self.oil_name]
T = self.max_water_fraction['temperatures']
wf = self.max_water_fraction['max_water_fraction']
logger.info(f'Using max water fractions {wf} for temperatures {T} for oiltype {self.oil_name}')
logger.info('Corresponding max water fraction from GNOME is '
f'{self.oiltype.gnome_oil["emulsion_water_fraction_max"]}')
else:
logger.info(f'Max water fraction not available for {self.oil_name}, using default')
except Exception as e:
logger.warning('Could not load max water content file')
print(e)
super(OpenOil, self).prepare_run()
[docs]
def oil_weathering_noaa(self):
'''Oil weathering scheme adopted from NOAA PyGNOME model:
https://github.com/NOAA-ORR-ERD/PyGnome
'''
logger.debug('NOAA oil weathering')
# C to K
self.environment.sea_water_temperature[
self.environment.sea_water_temperature < 100] += 273.15
#########################################################
# Update density and viscosity according to temperature
#########################################################
self.timer_start(
'main loop:updating elements:oil weathering:updating viscosities')
oil_viscosity = self.KinematicViscosity.at_temp(
self.environment.sea_water_temperature)
self.timer_end(
'main loop:updating elements:oil weathering:updating viscosities')
self.timer_start(
'main loop:updating elements:oil weathering:updating densities')
oil_density = self.Density.at_temp(self.environment.sea_water_temperature)
self.timer_end(
'main loop:updating elements:oil weathering:updating densities')
# Calculate emulsion density
self.elements.density = (
self.elements.water_fraction * self.sea_water_density() +
(1 - self.elements.water_fraction) * oil_density)
# Calculate emulsion viscosity
visc_f_ref = 0.84 # From PyGNOME
visc_curvfit_param = 1.5e3 # units are sec^0.5 / m
fw_d_fref = self.elements.water_fraction / visc_f_ref
kv1 = np.sqrt(oil_viscosity) * visc_curvfit_param
kv1[kv1 < 1] = 1
kv1[kv1 > 10] = 10
# NB: neglecting dispersed and biodegraded in calculation of fraction_evaporated
self.elements.fraction_evaporated = self.elements.mass_evaporated / (
self.elements.mass_oil + self.elements.mass_evaporated)
self.elements.viscosity = (
oil_viscosity * np.exp(kv1 * self.elements.fraction_evaporated) *
(1 + (fw_d_fref / (1.187 - fw_d_fref)))**2.49)
if self.get_config('processes:evaporation') is True:
self.timer_start(
'main loop:updating elements:oil weathering:evaporation')
self.evaporation_noaa()
self.timer_end(
'main loop:updating elements:oil weathering:evaporation')
if self.get_config('processes:emulsification') is True:
self.timer_start(
'main loop:updating elements:oil weathering:emulsification')
self.emulsification_noaa()
self.timer_end(
'main loop:updating elements:oil weathering:emulsification')
if self.get_config('processes:dispersion') is True:
self.timer_start(
'main loop:updating elements:oil weathering:dispersion')
self.disperse_noaa()
self.timer_end(
'main loop:updating elements:oil weathering:dispersion')
if self.get_config('processes:biodegradation') is True:
self.timer_start(
'main loop:updating elements:oil weathering:biodegradation')
self.biodegradation()
self.timer_end(
'main loop:updating elements:oil weathering:biodegradation')
if self.elements.mass_oil.min() < 0: # Should not happen
logger.warning('NEGATIVE OIL MASS!')
[docs]
def disperse_noaa(self):
logger.debug(' Calculating: dispersion - NOAA')
# From NOAA PyGnome model:
# https://github.com/NOAA-ORR-ERD/PyGnome/
c_disp = np.power(self.wave_energy_dissipation(), 0.57) * \
self.sea_surface_wave_breaking_fraction()
# Roy's constant
C_Roy = 2400.0 * np.exp(
-73.682 * np.sqrt(self.elements.viscosity / self.elements.density))
v_entrain = 3.9E-8
q_disp = C_Roy * c_disp * v_entrain / self.elements.density
fraction_dispersed = (q_disp * self.time_step.total_seconds() *
self.elements.density)
if fraction_dispersed.max() >= 1:
logger.warning(
'Dispersion fraction larger than 1 -> setting to 0.99')
fraction_dispersed[fraction_dispersed >= 1] = .99
oil_mass_loss = fraction_dispersed * self.elements.mass_oil
self.noaa_mass_balance['mass_components'][self.elements.ID - 1, :] = \
self.noaa_mass_balance['mass_components'][self.elements.ID - 1, :]*(1-fraction_dispersed[:, np.newaxis])
self.elements.mass_oil -= oil_mass_loss
self.elements.mass_dispersed += oil_mass_loss
[docs]
def plot_droplet_spectrum(self):
'''Plotting distribution of droplet radii, for debugging'''
plt.hist(self.elements.diameter / 2.0)
plt.show()
[docs]
def evaporation_noaa(self):
#############################################
# Evaporation, for elements at surface only
#############################################
logger.debug(' Calculating evaporation - NOAA')
surface = np.where(self.elements.z == 0)[0] # of active elements
if len(surface) == 0:
logger.debug('All elements submerged, no evaporation')
return
if self.elements.age_seconds[surface].min() > 3600 * 24:
logger.debug('All surface oil elements older than 24 hours, ' +
'skipping further evaporation.')
return
surfaceID = self.elements.ID[surface] - 1 # of any elements
# Area for each element, repeated for each component
volume = (self.elements.mass_oil[surface] /
self.elements.density[surface])
area = volume / self.elements.oil_film_thickness[surface]
evap_decay_constant = noaa.evap_decay_constant(
self.oiltype,
self.wind_speed()[surface],
self.environment.sea_water_temperature[surface], area,
self.noaa_mass_balance['mass_components'][surfaceID, :])
mass_remain = \
(self.noaa_mass_balance['mass_components'][surfaceID, :] *
np.exp(evap_decay_constant*self.time_step.total_seconds()))
self.elements.mass_evaporated[surface] += \
np.sum(self.noaa_mass_balance[
'mass_components'][surfaceID, :] - mass_remain, 1)
self.noaa_mass_balance['mass_components'][surfaceID, :] = \
mass_remain
self.elements.mass_oil[surface] = np.sum(mass_remain, 1)
[docs]
def emulsification_noaa(self):
#############################################
# Emulsification (surface only)
#############################################
logger.debug(' Calculating emulsification - NOAA')
emul_time = self.oiltype.bulltime
emul_constant = self.oiltype.bullwinkle
# Emulsify...
# f ((le_age >= emul_time && emul_time >= 0.) || frac_evap[i] >= emul_C && emul_C > 0.)
start_emulsion = np.where((
(self.elements.age_seconds >= emul_time) & (emul_time >= 0))
| ((self.elements.fraction_evaporated >= emul_constant)
& (emul_constant > 0)))[0]
if len(start_emulsion) == 0:
logger.debug(' Emulsification not yet started')
return
# max water content fraction - get from database
Y_max = np.atleast_1d(self.oiltype.emulsion_water_fraction_max)
if self.max_water_fraction is not None:
wf = self.max_water_fraction['max_water_fraction']
wft = self.max_water_fraction['temperatures']
if len(wf) == 1:
wf = [wf, wf]
wft = [wft, wft]
swt = self.environment.sea_water_temperature[start_emulsion] - 273.15 # to Celcius
weights = (wft[1] - swt) / (wft[1] - wft[0])
weights[swt>wft[1]] = 0
weights[swt<=wft[0]] = 1
max_water_fraction_sintef = weights*wf[0] + (1-weights)*wf[1]
if (Y_max - max_water_fraction_sintef).min() > 0:
logger.debug(
f'Overriding max water fraction {Y_max} with linear fit to SINTEF max values:'
f' T: {wft}, Fraction: {wf}')
Y_max = np.array(np.minimum(Y_max, max_water_fraction_sintef))
# emulsion
if Y_max.max() <= 0:
logger.debug('Oil does not emulsify, returning.')
return
# Constants for droplets
drop_min = 1.0e-6
drop_max = 1.0e-5
S_max = (6. / drop_min) * (Y_max / (1.0 - Y_max))
S_min = (6. / drop_max) * (Y_max / (1.0 - Y_max))
if self.oiltype.bulltime > 0: # User has set value
start_time = self.oiltype.bulltime * np.ones(len(start_emulsion))
else:
start_time = self.elements.age_seconds[start_emulsion]
# TODO: it should be possible to specify oil age at seeding
start_time[self.elements.age_seconds[start_emulsion] >=
0] = self.elements.bulltime[start_emulsion]
# Update droplet interfacial area
k_emul = noaa.water_uptake_coefficient(self.oiltype, self.wind_speed()[start_emulsion])
self.elements.interfacial_area[start_emulsion] = \
self.elements.interfacial_area[start_emulsion] + \
(k_emul*self.time_step.total_seconds()* np.exp((-k_emul/S_max)*(
self.elements.age_seconds[start_emulsion] - start_time)))
self.elements.interfacial_area[start_emulsion] = np.minimum(self.elements.interfacial_area[start_emulsion], S_max)
# Update water fraction
self.elements.water_fraction[start_emulsion] = (
self.elements.interfacial_area[start_emulsion] * drop_max / (6.0 +
(self.elements.interfacial_area[start_emulsion] * drop_max)))
self.elements.water_fraction[start_emulsion] = np.minimum(self.elements.water_fraction[start_emulsion], Y_max)
[docs]
def update_terminal_velocity(self,
Tprofiles=None,
Sprofiles=None,
z_index=None):
"""Calculate terminal velocity for oil droplets
according to:
* Tkalich et al. (2002): Vertical mixing of oil droplets by breaking waves
* Marine Pollution Bulletin 44, 1219-1229
If profiles of temperature and salt are passed into this function,
they will be interpolated from the profiles.
if not, T,S will be fetched from reader.
"""
g = 9.81 # ms-2
r = self.elements.diameter # NB: r is diameter, not radius
# Prepare interpolation of temp, salt
if not (Tprofiles is None and Sprofiles is None):
if z_index is None:
z_i = range(Tprofiles.shape[0]) # evtl. move out of loop
# evtl. move out of loop
z_index = interp1d(-self.environment_profiles['z'],
z_i,
bounds_error=False)
zi = z_index(-self.elements.z)
upper = np.maximum(np.floor(zi).astype(np.uint8), 0)
lower = np.minimum(upper + 1, Tprofiles.shape[0] - 1)
weight_upper = 1 - (zi - upper)
# Do interpolation of temp, salt if profiles were passed into
# this function, if not, use reader by calling self.environment
if Tprofiles is None:
T0 = self.environment.sea_water_temperature
else:
T0 = Tprofiles[upper, range(Tprofiles.shape[1])] * \
weight_upper + \
Tprofiles[lower, range(Tprofiles.shape[1])] * \
(1-weight_upper)
if Sprofiles is None:
S0 = self.environment.sea_water_salinity
else:
S0 = Sprofiles[upper, range(Sprofiles.shape[1])] * \
weight_upper + \
Sprofiles[lower, range(Sprofiles.shape[1])] * \
(1-weight_upper)
T0 = T0 - 273.15 # convert to Celcius - needed for the calcs in this method
rho_oil = self.elements.density
rho_water = self.sea_water_density(T=T0, S=S0)
# dynamic water viscosity
my_w = 0.001 * (1.7915 - 0.0538 * T0 + 0.007 *
(T0**(2.0)) - 0.0023 * S0)
# ~0.0014 kg m-1 s-1
# kinemativ water viscosity
ny_w = my_w / rho_water
rhopr = rho_oil / rho_water
# terminal velocity for low Reynolds numbers
kw = 2 * g * (1 - rhopr) / (9 * ny_w)
W = kw * (r/2)**2 # r is diameter so divide by 2 for radius
# check if we are in a high Reynolds number regime
Re = r * W / ny_w # r is diameter so no need to multiply by 2
highRe = np.where(Re > 50)
# Terminal velocity in high Reynolds numbers
kw = (16 * g * (1 - rhopr) / 3)**0.5
W2 = kw * (r/2)**0.5
W[highRe] = W2[highRe]
self.elements.terminal_velocity = W
[docs]
def oil_wave_entrainment_rate(self):
er = self.get_config('wave_entrainment:entrainment_rate')
if er == 'Li et al. (2017)':
entrainment_rate = oil_wave_entrainment_rate_li2017(
dynamic_viscosity=self.elements.viscosity *
self.elements.density,
oil_density=self.elements.density,
interfacial_tension=self.oil_water_interfacial_tension,
significant_wave_height=self.significant_wave_height(),
wave_breaking_fraction=self.sea_surface_wave_breaking_fraction(
),
sea_water_density=self.sea_water_density())
return entrainment_rate
else:
logger.error("no entrainment rate mechanism configured")
raise Exception("no entrainment rate mechanism configured")
[docs]
def prepare_vertical_mixing(self):
'''Calculate entrainment probability before main loop'''
self.oil_entrainment_probability = \
1 - np.exp(-self.oil_wave_entrainment_rate()*\
self.get_config('vertical_mixing:timestep'))
# Calculate a random droplet diameter for each particle,
# to be used if this particle gets entrained
self.droplet_diameter_if_entrained = \
self.get_wave_breaking_droplet_diameter()
# Uncomment lines below to plot droplet size distribution at each step
#import matplotlib.pyplot as plt
#plt.hist(self.droplet_diameter_if_entrained, 200)
#plt.gca().set_xscale("log")
#plt.gca().set_yscale("log")
#plt.show()
[docs]
def surface_wave_mixing(self, time_step_seconds):
"""Mix surface oil into water column."""
# Entrain oil into uppermost layer (whitecapping from waves)
# TODO: optimise this by only calculate for surface elements
surface = self.elements.z >= 0
random_number = np.random.uniform(0, 1, len(self.elements.z))
entrained = np.logical_and(
surface, random_number < self.oil_entrainment_probability)
# Intrusion depth for wave entrainment from
# Delvigne and Sweeney (1988), Li et al. (2017):
if entrained.sum() > 0:
logger.debug('Entraining %i of %i surface elements' %
(entrained.sum(), surface.sum()))
zb = 1.5 * self.significant_wave_height() # between 0 and zb
intrusion_depth = np.random.uniform(0, np.mean(zb),
entrained.sum())
self.elements.z[entrained] = -intrusion_depth
if self.keep_droplet_diameter is False:
# Give entrained elements a random diameter
self.elements.diameter[entrained] = \
self.droplet_diameter_if_entrained[entrained]
[docs]
def surface_stick(self):
"""set surfaced particles to exactly zero depth to let them form a slick """
surface = np.where(self.elements.z >= 0)
if len(surface[0]) > 0:
self.elements.z[surface] = 0.
[docs]
def get_wave_breaking_droplet_diameter(self):
dm = self.get_config('wave_entrainment:droplet_size_distribution')
if dm == 'Johansen et al. (2015)':
return self.get_wave_breaking_droplet_diameter_johansen2015()
elif dm == 'Li et al. (2017)':
return self.get_wave_breaking_droplet_diameter_liz2017()
else:
raise Exception("no wave entrainment droplet size distribution specified")
[docs]
def get_wave_breaking_droplet_diameter_liz2017(self):
# Li,Zhengkai, M. Spaulding, D. French-McCay, D. Crowley, J.R. Payne: "Development of a unified oil droplet size distribution model
# with application to surface breaking waves and subsea blowout releases considering dispersant effects" Mar. Pol. Bul.
# DOI: 10.1016/j.marpolbul.2016.09.008
# Should be prefered when the oil film thickness is unknown.
if not hasattr(self, 'droplet_spectrum_pdf'):
# Generate droplet spectrum as in Li (Zhengkai) et al. (2017)
# Bounds are hardcoded to 1 micron and 3mm
logger.debug('Generating wave breaking droplet size spectrum')
self.droplet_spectrum_diameter = np.linspace(1e-6, 3e-3, 1000000)
g = 9.81
interfacial_tension = self.oil_water_interfacial_tension
delta_rho = self.sea_water_density() - self.elements.density
d_o = 4 * (interfacial_tension / (delta_rho * g))**0.5
we = (self.sea_water_density() * g *
self.significant_wave_height() * d_o) / interfacial_tension
oh = self.elements.viscosity * self.elements.density * (
self.elements.density * interfacial_tension *
d_o)**-0.5 # From kin. to dyn. viscosity by * density
r = 1.791
p = 0.460
q = -0.518
dV_50 = d_o * r * (
1 + 10 * oh
)**p * we**q # median droplet diameter in volume distribution
sd = 0.4 # log standard deviation in log10 units
Sd = np.log(10) * sd # log standard deviation in natural log units
# TODO: calculation below with scalars, but we have arrays, with varying oil properties
# treat all particle in one go:
dV_50 = np.mean(dV_50) # mean log diameter
dN_50 = np.exp(
np.log(dV_50) - 3 *
Sd**2) # convert number distribution to volume distribution
logger.debug(
'Droplet distribution median diameter dV_50: %f, dN_50: %f ' %
(dV_50, np.mean(dN_50)))
spectrum = (np.exp(
-(np.log(self.droplet_spectrum_diameter) - np.log(dV_50))**2 /
(2 * Sd**2))) / (self.droplet_spectrum_diameter * Sd *
np.sqrt(2 * np.pi))
self.droplet_spectrum_pdf = spectrum / np.sum(spectrum)
if ~np.isfinite(np.sum(self.droplet_spectrum_pdf)) or \
np.abs(np.sum(self.droplet_spectrum_pdf) - 1) > 1e-6:
logger.warning('Could not update droplet diameters.')
return self.elements.diameter
else:
return np.random.choice(self.droplet_spectrum_diameter,
size=self.num_elements_active(),
p=self.droplet_spectrum_pdf)
[docs]
def get_wave_breaking_droplet_diameter_johansen2015(self):
# Johansen O, Reed M, Bodsberg NR, Natural dispersion revisited
# DOI: 10.1016/j.marpolbul.2015.02.026
# requires oil film thickness
if not hasattr(self, 'droplet_spectrum_pdf') or self.get_config(
'processes:update_oilfilm_thickness') is True:
# Generate droplet spectrum as in Johansen et al. (2015)
# Bounds are hardcoded to 1micron and 3mm
logger.debug('Generating wave breaking droplet size spectrum')
self.droplet_spectrum_diameter = np.linspace(1e-6, 3e-3, 1000000)
g = 9.81
interfacial_tension = self.oil_water_interfacial_tension
H = self.significant_wave_height(
) # fall height = 2 * wave amplitude
# Reyolds number (Eq. 7a from Johansen et al. 2015)
re = (self.elements.density * self.elements.oil_film_thickness *
(g * H)**0.5) / (self.elements.viscosity *
self.elements.density)
# Weber number (Eq. 7b from Johansen et al.2015)
we = (self.elements.density * self.elements.oil_film_thickness *
g * H) / interfacial_tension # Weber number
A = 2.251 # parameters from Johansen et al. 2015
Bp = 0.027
B = A * Bp
dN_50 = (A * self.elements.oil_film_thickness * we**-0.6) + (
B * self.elements.oil_film_thickness * re**-0.6)
# median droplet diameter in number distribution
sd = 0.4 # log standard deviation in log10 units
Sd = np.log(10) * sd # log standard deviation in natural log units
# Convert number distribution to volume distribution
dV_50 = np.exp(np.log(dN_50) + 3 * Sd**2)
# TODO: calculation below with scalars, but we have
# arrays, with varying oil properties
# treat all particle in one go:
dV_50 = np.mean(dV_50) # mean log diameter
logger.debug(
'Droplet distribution median diameter dV_50: %f, dN_50: %f ' %
(dV_50, np.mean(dN_50)))
spectrum = (np.exp(
-(np.log(self.droplet_spectrum_diameter) - np.log(dV_50))**2 /
(2 * Sd**2))) / (self.droplet_spectrum_diameter * Sd *
np.sqrt(2 * np.pi))
self.droplet_spectrum_pdf = spectrum / np.sum(spectrum)
if ~np.isfinite(np.sum(self.droplet_spectrum_pdf)) or \
np.abs(np.sum(self.droplet_spectrum_pdf) - 1) > 1e-6:
logger.warning('Could not update droplet diameters.')
return self.elements.diameter
else:
return np.random.choice(self.droplet_spectrum_diameter,
size=self.num_elements_active(),
p=self.droplet_spectrum_pdf)
[docs]
def resurface_elements(self, minimum_depth=None):
"""Oil elements reaching surface (or above) form slick, not droplet"""
surface = np.where(self.elements.z >= 0)[0]
self.elements.z[surface] = 0
[docs]
def advect_oil(self):
# Calculating various drift factors according to ice concentration
if hasattr(self.environment, 'sea_ice_area_fraction'):
A = self.environment.sea_ice_area_fraction
# According to
# Nordam T, Beegle-Krause CJ, Skancke J, Nepstad R, Reed M.
# Improving oil spill trajectory modelling in the Arctic.
# Mar Pollut Bull. 2019;140:65-74.
# doi:10.1016/j.marpolbul.2019.01.019
k_ice = (A - 0.3) / (0.8 - 0.3)
k_ice[A < 0.3] = 0
k_ice[A > 0.8] = 1
if k_ice.max() > 0:
logger.info(
'Ice concentration above 30%, using Nordam scheme for advection in ice'
)
# Using decreased Stokes drift according to
# Arneborg, L. (2017). Oil drift modellling in pack ice
# - Sensitivity of oil-in-ice parameters.
# Ocean Engineering 144 (2017) 340-350
factor_stokes = (0.7 - A) / 0.7
factor_stokes[A > 0.7] = 0
else:
k_ice = 0
factor_stokes = 1
# Simply move particles with ambient current
self.advect_ocean_current(factor=1 - k_ice)
# Wind drag for elements at ocean surface
self.advect_wind(factor=1 - k_ice)
# Stokes drift
self.stokes_drift(factor_stokes)
# Advect with ice
self.advect_with_sea_ice(factor=k_ice)
[docs]
def update(self):
"""Update positions and properties of oil particles."""
# TODO: move all config-checking inside respective methods
if self.get_config('processes:update_oilfilm_thickness') is True:
self.update_surface_oilfilm_thickness()
# Oil weathering (inherited from OpenOil)
self.oil_weathering()
# Turbulent Mixing
if self.get_config('drift:vertical_mixing') is True:
self.update_terminal_velocity()
self.vertical_mixing()
del self.droplet_spectrum_pdf
# Vertical advection
if self.get_config('drift:vertical_advection') is True:
self.vertical_advection()
# Horizontal advection (inherited from OpenOil)
self.advect_oil()
[docs]
def get_oil_budget(self):
"""Get oil budget for the current simulation
The oil budget consists of the following categories:
* surface: the sum of variable mass_oil for all active elements where z = 0
* submerged: the sum of variable mass_oil for all active elements where z < 0
* stranded: the sum of variable mass_oil for all elements which are stranded
* evaporated: the sum of variable mass_evaporated for all elements
* dispersed: the sum of variable mass_dispersed for all elements
The sum (total mass) shall equal the mass released. Note that the mass of oil
is conserved, whereas the volume may change with changes in density and
water uptake (emulsification). Therefore mass should be used for budgets,
eventually converted to volume (by dividing on density) in the final step
before presentation.
Note that mass_oil is the mass of pure oil. The mass of oil emulsion
(oil containing entrained water droplets) can be calculated as:
.. code::
mass_emulsion = mass_oil / (1 - water_fraction)
I.e. water_fraction = 0 means pure oil, water_fraction = 0.5 means mixture of
50% oil and 50% water, and water_fraction = 0.9 (which is maximum)
means 10% oil and 90% water.
"""
if self.time_step.days < 0: # Backwards simulation
return None
z, dummy = self.get_property('z')
mass_oil, status = self.get_property('mass_oil')
density = self.get_property('density')[0][0, 0]
if 'stranded' not in self.status_categories:
self.status_categories.append('stranded')
mass_submerged = np.ma.masked_where(
((status == self.status_categories.index('stranded')) |
(z == 0.0)), mass_oil)
mass_submerged = np.ma.sum(mass_submerged, axis=1).filled(0)
mass_surface = np.ma.masked_where(
((status == self.status_categories.index('stranded')) | (z < 0.0)),
mass_oil)
mass_surface = np.ma.sum(mass_surface, axis=1).filled(0)
mass_stranded = np.ma.sum(np.ma.masked_where(
status != self.status_categories.index('stranded'), mass_oil),
axis=1).filled(0)
mass_evaporated, status = self.get_property('mass_evaporated')
mass_evaporated = np.sum(mass_evaporated, axis=1).filled(0)
mass_dispersed, status = self.get_property('mass_dispersed')
mass_dispersed = np.sum(mass_dispersed, axis=1).filled(0)
mass_biodegraded, status = self.get_property('mass_biodegraded')
mass_biodegraded = np.sum(mass_biodegraded, axis=1).filled(0)
oil_budget = {
'oil_density':
density,
'mass_dispersed':
mass_dispersed,
'mass_submerged':
mass_submerged,
'mass_surface':
mass_surface,
'mass_stranded':
mass_stranded,
'mass_evaporated':
mass_evaporated,
'mass_biodegraded':
mass_biodegraded,
'mass_total': (mass_dispersed + mass_submerged + mass_surface +
mass_stranded + mass_evaporated + mass_biodegraded)
}
return oil_budget
[docs]
def plot_oil_budget(self,
filename=None,
ax=None,
show_watercontent_and_viscosity=True,
show_wind_and_current=True):
if self.time_step.days < 0: # Backwards simulation
fig = plt.figure(figsize=(10, 6.))
plt.text(0.1, 0.5, 'Oil weathering deactivated for '
'backwards simulations')
plt.axis('off')
if filename is not None:
plt.savefig(filename)
plt.close()
else:
plt.show()
return
b = self.get_oil_budget()
oil_budget = np.row_stack(
(b['mass_dispersed'], b['mass_submerged'], b['mass_surface'],
b['mass_stranded'], b['mass_evaporated'], b['mass_biodegraded']))
oil_density = b['oil_density']
budget = np.cumsum(oil_budget, axis=0)
time, time_relative = self.get_time_array()
time = np.array([t.total_seconds() / 3600. for t in time_relative])
if ax is None:
# Left axis showing oil mass
nrows = 1
if show_watercontent_and_viscosity is True:
nrows = nrows + 1
if show_wind_and_current is True:
nrows = nrows + 1
fig, axs = plt.subplots(
nrows=nrows, ncols=1,
figsize=(10, 6. + (nrows - 1) * 3)) # Suitable aspect ratio
#ax1 = fig.add_subplot(nrows=nrows, 1, 1)
if nrows == 1:
ax1 = axs
elif nrows >= 2:
ax1 = axs[0]
if show_watercontent_and_viscosity is True:
self.plot_oil_watercontent_and_viscosity(ax=axs[1], show=False)
if show_wind_and_current is True:
self.plot_environment(ax=axs[nrows - 1], show=False)
else:
ax1 = ax
# Hack: make some emply plots since fill_between does not support label
if np.sum(b['mass_dispersed']) > 0:
ax1.add_patch(
plt.Rectangle((0, 0),
0,
0,
color='darkslategrey',
label='dispersed'))
ax1.fill_between(time, 0, budget[0, :], facecolor='darkslategrey')
if np.sum(b['mass_submerged']) > 0:
ax1.add_patch(
plt.Rectangle((0, 0),
0,
0,
color='darkblue',
label='submerged'))
ax1.fill_between(time,
budget[0, :],
budget[1, :],
facecolor='darkblue')
if np.sum(b['mass_surface']) > 0:
ax1.add_patch(
plt.Rectangle((0, 0), 0, 0, color='royalblue',
label='surface'))
ax1.fill_between(time,
budget[1, :],
budget[2, :],
facecolor='royalblue')
if np.sum(b['mass_stranded']) > 0:
ax1.add_patch(
plt.Rectangle((0, 0), 0, 0, color='black', label='stranded'))
ax1.fill_between(time,
budget[2, :],
budget[3, :],
facecolor='black')
if np.sum(b['mass_evaporated']) > 0:
ax1.add_patch(
plt.Rectangle((0, 0),
0,
0,
color='skyblue',
label='evaporated'))
ax1.fill_between(time,
budget[3, :],
budget[4, :],
facecolor='skyblue')
if np.sum(b['mass_biodegraded']) > 0:
ax1.add_patch(
plt.Rectangle((0, 0),
0,
0,
color='indigo',
label='biodegraded'))
ax1.fill_between(time,
budget[4, :],
budget[5, :],
facecolor='indigo')
ax1.set_ylim([0, budget.max()])
ax1.set_xlim([0, time.max()])
ax1.set_ylabel('Mass oil [%s]' %
self.elements.variables['mass_oil']['units'])
ax1.set_xlabel('Time [hours]')
# Right axis showing volume
ax2 = ax1.twinx()
mass_total = b['mass_total'][-1]
ax2.set_ylim([0, mass_total / oil_density])
ax2.set_ylabel('Volume oil [m3]')
plt.title('%s (%.1f kg/m3) - %s to %s' %
(self.get_oil_name(), oil_density,
self.start_time.strftime('%Y-%m-%d %H:%M'),
self.time.strftime('%Y-%m-%d %H:%M')))
# Shrink current axis's height by 10% on the bottom
#box = ax1.get_position()
#ax1.set_position(
# [box.x0, box.y0 + box.height * 0.15, box.width, box.height * 0.85])
#ax2.set_position(
# [box.x0, box.y0 + box.height * 0.5, box.width, box.height * 0.6])
ax1.legend(bbox_to_anchor=(0., -0.1, 1., -0.04),
loc=1,
ncol=6,
mode="expand",
borderaxespad=0.,
fontsize=10)
if filename is not None:
plt.savefig(filename)
plt.close()
else:
plt.show()
[docs]
def get_oil_name(self):
try:
return self.get_config('seed:oil_type')
except: # fallback if importing old files
return 'unknown oiltype'
[docs]
def cumulative_oil_entrainment_fraction(self):
'''Returns the fraction of oil elements which has been entrained vs time'''
z = self.get_property('z')[0].copy()
z = np.ma.masked_where(z == 0, z)
me = np.ma.notmasked_edges(z, axis=0)
maskfirst = me[0][0]
maskrow = me[0][1]
z = z * 0
for mf, mr in zip(maskfirst, maskrow):
z[mf:z.shape[0], mr] = 1 # has been entrained
totentrained = np.sum(z, 1)
cumulative_fraction_entrained = np.sum(z, 1) / z.shape[1]
return cumulative_fraction_entrained
[docs]
def plot_oil_watercontent_and_viscosity(self, ax=None, show=True):
if ax is None:
fig, ax = plt.subplots()
import matplotlib.dates as mdates
time, time_relative = self.get_time_array()
time = np.array([t.total_seconds() / 3600. for t in time_relative])
kin_viscosity = self.history['viscosity']
dyn_viscosity = kin_viscosity * self.history['density'] * 1000 # unit of mPas
dyn_viscosity_mean = dyn_viscosity.mean(axis=0)
dyn_viscosity_std = dyn_viscosity.std(axis=0)
density = self.history['density'].mean(axis=0)
density_std = self.history['density'].std(axis=0)
watercontent = self.history['water_fraction'].mean(axis=0)*100
watercontent_std = self.history['water_fraction'].std(axis=0)*100
ax.plot(time,
dyn_viscosity_mean,
'g',
lw=2,
label='Emulsion viscosity')
ax.fill_between(time,
dyn_viscosity_mean - dyn_viscosity_std,
dyn_viscosity_mean + dyn_viscosity_std,
color='g',
alpha=0.5)
ax.set_ylim([0, max(dyn_viscosity_mean + dyn_viscosity_std)])
ax.set_ylabel(r'Emulsion viscosity [cPoise] / [mPas]', color='g')
ax.tick_params(axis='y', colors='g')
axb = ax.twinx()
axb.plot(time, watercontent, 'b', lw=2, label='Water content')
axb.fill_between(time,
watercontent - watercontent_std,
watercontent + watercontent_std,
color='b', alpha=0.5)
ax.set_xlim([0, time.max()])
ax.set_xlabel('Time [hours]')
axb.set_ylim([0, 100])
axb.set_ylabel(r'Water content [%]', color='b')
axb.tick_params(axis='y', colors='b')
ax.legend(loc='upper left')
axb.legend(loc='upper center')
if show is True:
plt.show()
[docs]
def set_oiltype(self, oiltype):
"""
Sets the oil type by specifying the name, the first match will be chosen. See the `ADIOS database <https://adios.orr.noaa.gov/oils>`_ for a list. OpenDrift provides a small set of extra oils.
"""
if self.get_config('seed:oil_type') != oiltype:
self.__set_seed_config__('seed:oil_type', oiltype)
logger.info(f'setting oil_type to: {oiltype}')
oiltype = adios.oil_name_alias.get(oiltype, oiltype)
self.oil_name = oiltype
if self.oil_weathering_model == 'noaa':
self.oiltype = adios.find_full_oil_from_name(self.oil_name)
if not self.oiltype.valid():
logger.error(
f"{self.oiltype} is not a valid oil for Opendrift simulations"
)
raise ValueError()
else:
raise ValueError("unsupported oil weathering model")
[docs]
def set_oiltype_by_id(self, oiltypeid):
"""
Sets the oil type by specifying the ADIOS ID. See the `ADIOS database <https://adios.orr.noaa.gov/oils>`_ for a list. OpenDrift provides a small set of extra oils.
"""
if self.oil_weathering_model == 'noaa':
self.oiltype = adios.get_full_oil_from_id(oiltypeid)
self.oil_name = self.oiltype.name
if not self.oiltype.valid():
logger.error(
f"{self.oiltype} is not a valid oil for Opendrift simulations"
)
raise ValueError()
else:
raise ValueError("unsupported oil weathering model")
[docs]
def set_oiltype_by_json(self, json):
"""
Sets the oil type by specifing a JSON dict. The format should be the same as the ADIOS database. See the `ADIOS database <https://adios.orr.noaa.gov/oils>`_ for a list.
"""
if self.oil_weathering_model == 'noaa':
#o = { 'data': { 'attributes' : json } }
#o['data']['_id'] = o['data']['attributes']['oil_id']
#o['data']['attributes']['metadata']['location'] = 'Norway'
self.oiltype = adios.oil.OpendriftOil(json)
self.oil_name = self.oiltype.name
if not self.oiltype.valid():
logger.error(
f"{self.oiltype} is not a valid oil for Opendrift simulations"
)
raise ValueError()
else:
raise ValueError("unsupported oil weathering model")
[docs]
def set_oiltype_from_file(self, path):
"""
Sets the oil type by specifing a JSON file. The format should be the same as the ADIOS database. See the `ADIOS database <https://adios.orr.noaa.gov/oils>`_ for a list.
>>> o = OpenOil()
>>> o.set_oiltype_from_file('opendrift/models/openoil/adios/extra_oils/AD04001.json')
"""
if self.oil_weathering_model == 'noaa':
import json
with open(path, 'r') as fd:
j = json.load(fd)
self.set_oiltype_by_json(j)
else:
raise ValueError("unsupported oil weathering model")
[docs]
def seed_elements(self, *args, **kwargs):
if len(args) == 2:
kwargs['lon'] = args[0]
kwargs['lat'] = args[1]
args = {}
if 'number' not in kwargs:
number = self.get_config('seed:number')
else:
number = kwargs['number']
if 'diameter' in kwargs:
logger.info('Droplet diameter is provided, and will '
'be kept constant during simulation')
self.keep_droplet_diameter = True
else:
self.keep_droplet_diameter = False
if 'z' not in kwargs or kwargs['z'] is None:
if self.get_config('seed:seafloor') is True:
kwargs['z'] = 'seafloor'
else:
kwargs['z'] = self.get_config('seed:z')
if isinstance(kwargs['z'], str) and \
kwargs['z'][0:8] == 'seafloor':
z = -np.ones(number)
else:
z = np.atleast_1d(kwargs['z'])
if len(z) == 1:
z = z * np.ones(number) # Convert scalar z to array
subsea = z < 0
if np.sum(subsea) > 0 and 'diameter' not in kwargs:
dsd = self.get_config('seed:droplet_size_distribution')
if dsd == 'uniform':
# Droplet min and max for particles seeded below sea surface
sub_dmin = self.get_config('seed:droplet_diameter_min_subsea')
sub_dmax = self.get_config('seed:droplet_diameter_max_subsea')
logger.info('Using uniform droplet size distribution between %s and %s m for '
'elements seeded below sea surface.' %
(sub_dmin, sub_dmax))
kwargs['diameter'] = np.random.uniform(sub_dmin, sub_dmax, number)
elif dsd == 'normal':
# Droplet mu and sigma for particles seeded below sea surface
sub_mu = self.get_config('seed:droplet_diameter_mu')
sub_sigma = self.get_config('seed:droplet_diameter_sigma')
logger.info('Using normal droplet size distribution with '
'mu = %s and sigma = %s m for elements seeded below sea surface.' %
(sub_mu, sub_sigma))
kwargs['diameter'] = np.random.normal(sub_mu, sub_sigma, number)
elif dsd == 'lognormal':
# Droplet mu and sigma for particles seeded below sea surface
sub_mu = self.get_config('seed:droplet_diameter_mu')
sub_sigma = self.get_config('seed:droplet_diameter_sigma')
logger.info('Using lognormal droplet size distribution with '
'mu = %s and sigma = %s m for elements seeded below sea surface.' %
(sub_mu, sub_sigma))
# From numpy.random.lognormal:
# "Note that the mean and standard deviation are not the values for the
# distribution itself, but of the underlying normal distribution it is derived from."
# So we need to compute the input to the function from the mean and
# standard deviation of the data we want to generate (assumed as input)
sub_sigma2 = sub_sigma**2
sub_sigma2_lognormal = np.log(sub_sigma2/sub_mu**2 + 1)
sub_mu_lognormal = np.log(sub_mu) - sub_sigma2_lognormal/2
sub_sigma_lognormal = sub_sigma2_lognormal**0.5
kwargs['diameter'] = np.random.lognormal(sub_mu_lognormal, sub_sigma_lognormal, number)
# check it worked
# print('median = '+str(np.median(kwargs['diameter'])))
# print('mean = '+str(np.mean(kwargs['diameter'])))
# print('sigma = '+str(np.std(kwargs['diameter'])))
else:
raise Exception("no valid initial subsea droplet size distribution specified")
if 'oiltype' in kwargs:
logger.warning(
'Seed argument *oiltype* is deprecated, use *oil_type* instead'
)
kwargs['oil_type'] = kwargs['oiltype']
del kwargs['oiltype']
if 'oil_type' in kwargs:
if self.get_config('seed:oil_type') != kwargs['oil_type']:
self.__set_seed_config__('seed:oil_type', kwargs['oil_type'])
del kwargs['oil_type']
else:
logger.info('Oil type not specified, using default: ' +
self.get_config('seed:oil_type'))
self.set_oiltype(self.get_config('seed:oil_type'))
if self.oil_weathering_model == 'noaa':
self.Density = Density(self.oiltype.oil)
self.KinematicViscosity = KinematicViscosity(self.oiltype.oil)
oil_density = self.Density.at_temp(285)
oil_viscosity = self.KinematicViscosity.at_temp(285)
logger.info(
'Using density %s and viscosity %s of oiltype %s' %
(oil_density, oil_viscosity, self.get_config('seed:oil_type')))
kwargs['density'] = oil_density
kwargs['viscosity'] = oil_viscosity
if 'm3_per_hour' in kwargs:
m3_per_hour = kwargs['m3_per_hour']
del kwargs['m3_per_hour']
else:
m3_per_hour = self.get_config('seed:m3_per_hour')
if 'number' in kwargs:
num_elements = kwargs['number']
else:
num_elements = self.get_config('seed:number')
time = kwargs['time']
if type(time) is list:
duration_hours = ((time[1] - time[0]).total_seconds()) / 3600
if duration_hours == 0:
duration_hours = 1.
else:
duration_hours = 1. # For instantaneous spill, we use 1h
kwargs['mass_oil'] = (m3_per_hour * duration_hours / num_elements *
kwargs['density'])
self.store_oil_seed_metadata(**kwargs)
super(OpenOil, self).seed_elements(*args, **kwargs)
[docs]
def seed_cone(self, *args, **kwargs):
if 'oiltype' in kwargs:
logger.warning(
'Seed argument *oiltype* is deprecated, use *oil_type* instead'
)
kwargs['oil_type'] = kwargs['oiltype']
del kwargs['oiltype']
self.store_oil_seed_metadata(**kwargs)
super(OpenOil, self).seed_cone(*args, **kwargs)
[docs]
def seed_from_gml(self, gmlfile, num_elements=1000, *args, **kwargs):
"""Read oil slick contours from GML file, and seed particles within."""
# Specific imports
import datetime
from matplotlib.path import Path
from xml.etree import ElementTree
from matplotlib.patches import Polygon
import pyproj
namespaces = {
'od': 'http://cweb.ksat.no/cweb/schema/geoweb/oil',
'gml': 'http://www.opengis.net/gml'
}
slicks = []
with open(gmlfile, 'rt') as e:
tree = ElementTree.parse(e)
pos1 = 'od:oilDetectionMember/od:oilDetection/od:oilSpill/gml:Polygon'
pos2 = 'gml:exterior/gml:LinearRing/gml:posList'
# This retrieves some other types of patches, found in some files only
# Should be combines with the above, to get all patches
#pos1 = 'od:oilDetectionMember/od:oilDetection/od:oilSpill/
# gml:Surface/gml:polygonPatches'
#pos2 = 'gml:PolygonPatch/gml:exterior/gml:LinearRing/gml:posList'
# Find detection time
time_pos = 'od:oilDetectionMember/od:oilDetection/od:detectionTime'
oil_time = datetime.datetime.strptime(
tree.find(time_pos, namespaces).text, '%Y-%m-%dT%H:%M:%S.%fZ')
for patch in tree.findall(pos1, namespaces):
pos = patch.find(pos2, namespaces).text
c = np.array(pos.split()).astype(float)
lon = c[0::2]
lat = c[1::2]
slicks.append(Polygon(list(zip(lon, lat))))
# Find boundary and area of all patches
lons = np.array([])
lats = lons.copy()
for slick in slicks:
ext = slick.get_extents()
lons = np.append(lons, [ext.xmin, ext.xmax])
lats = np.append(lats, [ext.ymin, ext.ymax])
# Make a stereographic projection centred on the polygon
lonmin = lons.min()
lonmax = lons.max()
latmin = lats.min()
latmax = lats.max()
# Place n points within the polygons
proj = pyproj.Proj(
'+proj=aea +lat_1=%f +lat_2=%f +lat_0=%f +lon_0=%f' %
(latmin, latmax, (latmin + latmax) / 2, (lonmin + lonmax) / 2))
slickarea = np.array([])
for slick in slicks:
lonlat = slick.get_xy()
lon = lonlat[:, 0]
lat = lonlat[:, 1]
x, y = proj(lon, lat)
area_of_polygon = 0.0
for i in range(-1, len(x) - 1):
area_of_polygon += x[i] * (y[i + 1] - y[i - 1])
area_of_polygon = abs(area_of_polygon) / 2.0
slickarea = np.append(slickarea, area_of_polygon) # in m2
# Make points
deltax = np.sqrt(np.sum(slickarea) / num_elements)
lonpoints = np.array([])
latpoints = np.array([])
for i, slick in enumerate(slicks):
lonlat = slick.get_xy()
lon = lonlat[:, 0]
lat = lonlat[:, 1]
x, y = proj(lon, lat)
xvec = np.arange(x.min(), x.max(), deltax)
yvec = np.arange(y.min(), y.max(), deltax)
x, y = np.meshgrid(xvec, yvec)
lon, lat = proj(x, y, inverse=True)
lon = lon.ravel()
lat = lat.ravel()
points = np.c_[lon, lat]
ind = Path(slick.xy).contains_points(points)
lonpoints = np.append(lonpoints, lon[ind])
latpoints = np.append(latpoints, lat[ind])
# Finally seed at found positions
kwargs['lon'] = lonpoints
kwargs['lat'] = latpoints
self.seed_elements(time=oil_time, **kwargs)
[docs]
def seed_from_geotiff_thickness(self,
filename,
number=50000,
*args,
**kwargs):
'''Seed from files as provided by Prof. Chuanmin Hu'''
from osgeo import gdal, ogr, osr
if not 'time' in kwargs:
try: # get time from filename
time = datetime.strptime(filename[-28:-13], '%Y%m%d.%H%M%S')
logger.info('Parsed time from filename: %s' % time)
except:
time = datetime.now()
logger.warning('Could not pase time from filename, '
'using present time: %s' % time)
ds = gdal.Open(filename)
srcband = ds.GetRasterBand(1)
data = srcband.ReadAsArray()
thickness_microns = [0.04, 0.44, 4.4, 16] # NB: approximate
categories = [1, 2, 3, 4] # categories
# Make memory raster bands for each category
memrastername = filename + '.mem'
memdriver = gdal.GetDriverByName('MEM')
mem_ds = memdriver.CreateCopy(memrastername, ds)
for cat in categories:
mem_ds.AddBand(gdal.GDT_Byte)
mem_band = mem_ds.GetRasterBand(cat)
mem_band.WriteArray(data == cat)
# Make memory polygons for each category
drv = ogr.GetDriverByName('MEMORY')
mem_vector_ds = [0] * len(categories)
mem_vector_layers = [0] * len(categories)
for cat in categories:
memshapename = filename + '%i.mem' % cat
mem_vector_ds[cat - 1] = drv.CreateDataSource(memshapename)
mem_vector_layers[cat-1] = \
mem_vector_ds[cat-1].CreateLayer(
'thickness%i' % cat, srs=None)
gdal.Polygonize(mem_ds.GetRasterBand(cat),
None,
mem_vector_layers[cat - 1],
-1, [],
callback=None)
total_area = np.zeros(len(categories))
layers = [0] * len(categories)
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(4269)
tgt_srs = osr.SpatialReference()
tgt_srs.ImportFromEPSG(3857)
transform = osr.CoordinateTransformation(src_srs, tgt_srs)
for cat in categories:
memshapename = filename + '%i.shp' % cat
layers[cat - 1] = mem_vector_layers[cat - 1]
areas = np.zeros(layers[cat - 1].GetFeatureCount())
for i, feature in enumerate(layers[cat - 1]):
geom = feature.GetGeometryRef()
geom.Transform(transform) # To get area in m2
areas[i] = geom.GetArea()
# Delete largest polygon, which is outer border
outer = np.where(areas == max(areas))[0]
areas[outer] = 0
total_area[cat - 1] = np.sum(areas)
layers[cat - 1].DeleteFeature(outer)
layers[cat - 1].ResetReading()
# Calculate how many elements to be seeded for each category
areas_weighted = total_area * thickness_microns
numbers = number * areas_weighted / np.sum(areas_weighted)
numbers = np.round(numbers).astype(int)
oil_density = 1000
mass_oil = (total_area * thickness_microns / 1e6) * oil_density
for i, num in enumerate(numbers):
self.seed_from_shapefile([mem_vector_layers[i]],
oil_film_thickness=thickness_microns[i] /
1000000.,
mass_oil=mass_oil[i] / num,
number=num,
time=time,
*args,
**kwargs)
[docs]
def _substance_name(self):
return self.get_oil_name()