Note
Go to the end to download the full example code
Droplet distribution
Plotting different droplet size distributions used in Opendrift
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
from opendrift.models.openoil import OpenOil
droplet size interval for plotting
dmin = 1.e-6
dmax = 3.e-3
simulation with Johansen et al. (2015) droplet spectrum
o = OpenOil(loglevel=20, weathering_model='noaa')
o.set_config('environment:fallback:land_binary_mask', 0)
o.set_config('environment:fallback:x_sea_water_velocity', -.2)
o.set_config('environment:fallback:y_sea_water_velocity', 0)
o.set_config('environment:fallback:x_wind', 6.)
o.set_config('environment:fallback:y_wind', 0)
o.set_config('environment:fallback:sea_water_temperature', 5.)
o.set_config('environment:fallback:sea_surface_wave_stokes_drift_x_velocity', .3)
o.set_config('environment:fallback:sea_surface_wave_stokes_drift_y_velocity', 0)
o.set_config('wave_entrainment:droplet_size_distribution', 'Johansen et al. (2015)')
o.set_config('processes:evaporation', False)
o.set_config('processes:dispersion', False)
o.seed_elements(lon=4, lat=60, time=datetime.utcnow(), number=10000, radius=100,
z=0, oil_type='TROLL, STATOIL', oil_film_thickness=0.005)
o.run(duration=timedelta(hours=2), time_step=3600)
droplet_diameters = o.elements.diameter
sd = 0.4
Sd = np.log(10.)*sd
DV_50 = np.median(droplet_diameters)
DN_50 = np.exp( np.log(DV_50) - 3*Sd**2 )
print('DV_50: %f' % DV_50)
print('DN_50: %f' % DN_50)
17:06:27 INFO opendrift.models.basemodel:533: OpenDriftSimulation initialised (version 1.11.2 / v1.11.2-17-g1eadfa7)
17:06:27 INFO opendrift.models.openoil.adios.dirjs:90: Querying ADIOS database for oil: TROLL, STATOIL
17:06:27 INFO opendrift.models.openoil.openoil:1720: Using density 891.412528 and viscosity 0.0032285361678330597 of oiltype TROLL, STATOIL
17:06:27 INFO opendrift.models.basemodel.environment:218: Adding a dynamical landmask with max. priority based on assumed maximum speed of 1.0 m/s. Adding a customised landmask may be faster...
17:06:33 INFO opendrift.models.basemodel.environment:245: Fallback values will be used for the following variables which have no readers:
17:06:33 INFO opendrift.models.basemodel.environment:248: x_sea_water_velocity: -0.200000
17:06:33 INFO opendrift.models.basemodel.environment:248: y_sea_water_velocity: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: x_wind: 6.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: y_wind: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_surface_height: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: upward_sea_water_velocity: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_surface_wave_significant_height: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_surface_wave_stokes_drift_x_velocity: 0.300000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_surface_wave_stokes_drift_y_velocity: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_surface_wave_period_at_variance_spectral_density_maximum: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_surface_wave_mean_period_from_variance_spectral_density_second_frequency_moment: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_ice_area_fraction: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_ice_x_velocity: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_ice_y_velocity: 0.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_water_temperature: 5.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_water_salinity: 34.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: sea_floor_depth_below_sea_level: 10000.000000
17:06:33 INFO opendrift.models.basemodel.environment:248: ocean_vertical_diffusivity: 0.020000
17:06:33 INFO opendrift.models.basemodel.environment:248: ocean_mixed_layer_thickness: 50.000000
17:06:33 INFO opendrift.models.basemodel:909: Using existing reader for land_binary_mask
17:06:33 INFO opendrift.models.basemodel:920: All points are in ocean
17:06:33 INFO opendrift.models.openoil.openoil:688: Oil-water surface tension is 0.032037 Nm
17:06:33 INFO opendrift.models.openoil.openoil:697: Using max water fractions [0.74, 0.75] for temperatures [5.0, 15.0] for oiltype TROLL, STATOIL
17:06:33 INFO opendrift.models.openoil.openoil:698: Corresponding max water fraction from GNOME is 0.5262214946654028
17:06:33 INFO opendrift.models.basemodel:2009: 2024-04-12 17:06:27.779005 - step 1 of 2 - 10000 active elements (0 deactivated)
17:06:33 INFO opendrift.models.basemodel:2009: 2024-04-12 18:06:27.779005 - step 2 of 2 - 10000 active elements (0 deactivated)
/root/project/examples/example_droplet_distribution_compareJohansen2015.py:42: RuntimeWarning: divide by zero encountered in log
DN_50 = np.exp( np.log(DV_50) - 3*Sd**2 )
DV_50: 0.000000
DN_50: 0.000000
Plotting
plt.figure(figsize=[14,14])
plt.subplot(3, 2, 1)
nVpdf, binsV, patches = plt.hist(droplet_diameters, 100, range=(dmin,dmax), align='mid')
plt.xlabel('Droplet diameter d [m]', fontsize=8)
plt.ylabel('V(d)', fontsize=8)
plt.title('volume spectrum\nJohansen et al. (2015) distribution in OpenOil', fontsize=10)
plt.subplot(3,2,2)
nVcum, binsV, patches = plt.hist(droplet_diameters, 100, range=(dmin,dmax), align='mid', cumulative=True)
plt.xlabel('Droplet diameter d [m]', fontsize=8)
plt.ylabel('V(d)', fontsize=8)
plt.title('cumulative volume spectrum', fontsize=10)
Text(0.5, 1.0, 'cumulative volume spectrum')
calculate number spectrum from volume spectrum
d = 0.5* (binsV[1:] + binsV[:-1])
V = 4./3. * np.pi * (d/2.)**3
Npdf = nVpdf / V
Ncum = np.cumsum(Npdf)
calculate theoretical cumulative Number distribution
spectrum = (np.exp(-(np.log(d) - np.log(DN_50))**2 / (2 * Sd**2))) / (d * Sd * np.sqrt(2 * np.pi)) # from OpenOil median diameter
DN_50_johansen = 0.000483
spectrumJ = (np.exp(-(np.log(d) - np.log(DN_50_johansen))**2 / (2 * Sd**2))) / (d * Sd * np.sqrt(2 * np.pi)) # from parameters in Johansen et al. 2015
spectrum_cum = np.cumsum(spectrum)
spectrumJ_cum = np.cumsum(spectrumJ)
# calculate theoretical number distribution pdf
#spectrum_pdf = spectrum / np.sum(spectrum)
#spectrumJ_pdf = spectrumJ / np.sum(spectrumJ)
plt.subplot(3,2,3)
plt.plot(d, Npdf/np.sum(Npdf), lw=2)
#plt.plot(d, spectrum/np.sum(spectrum), label='curve fit from OpenOil', lw=2)
plt.plot(d, spectrumJ/np.sum(spectrumJ), label='curve fit from Johansen et al. 2015', lw=2)
plt.xlabel('Droplet diameter d [m]', fontsize=8)
plt.ylabel('N(d)', fontsize=8)
plt.title('number spectrum', fontsize=10)
plt.subplot(3,2,4)
plt.plot(d, Ncum/np.max(Ncum), label='OpenOil result', lw=2)
#plt.plot(d, spectrum_cum/np.max(spectrum_cum), label='OpenOil par.', lw=2)
plt.plot(d, spectrumJ_cum/np.max(spectrumJ_cum), label='Johansen et al. 2015', lw=2)
plt.xlabel('Droplet diameter d [m]', fontsize=8)
plt.ylabel('N(d)', fontsize=8)
plt.title('cumulative number spectrum', fontsize=10)
plt.legend(loc='lower right')
plt.subplot(3,2,5)
plt.loglog(d, Npdf/np.sum(Npdf), lw=2)
#plt.loglog(d, spectrum/np.sum(spectrum), label='curve fit from OpenOil', lw=2)
plt.loglog(d, spectrumJ/np.sum(spectrumJ), label='curve fit from Johansen et al. 2015', lw=2)
plt.xlabel('Droplet diameter d [m]', fontsize=8)
plt.ylabel('N(d)', fontsize=8)
plt.title('number spectrum', fontsize=10)
plt.subplot(3,2,6)
plt.semilogx(d, Ncum/np.max(Ncum), label='OpenOil result', lw=2)
#plt.semilogx(d, spectrum_cum/np.max(spectrum_cum), label='OpenOil par.', lw=2)
plt.semilogx(d, spectrumJ_cum/np.max(spectrumJ_cum), label='Johansen et al. 2015', lw=2)
plt.xlabel('Droplet diameter d [m]', fontsize=8)
plt.ylabel('N(d)', fontsize=8)
plt.title('cumulative number spectrum', fontsize=10)
#plt.legend(loc='lower right')
plt.show()
/root/project/examples/example_droplet_distribution_compareJohansen2015.py:71: RuntimeWarning: divide by zero encountered in log
spectrum = (np.exp(-(np.log(d) - np.log(DN_50))**2 / (2 * Sd**2))) / (d * Sd * np.sqrt(2 * np.pi)) # from OpenOil median diameter
Total running time of the script: (0 minutes 8.777 seconds)