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)

20:13:32 INFO    opendrift.models.basemodel:515: OpenDriftSimulation initialised (version 1.12.0 / v1.12.0-26-g390e945)
20:13:32 INFO    opendrift.models.openoil.adios.dirjs:86: Querying ADIOS database for oil: TROLL, STATOIL
20:13:32 INFO    opendrift.models.openoil.openoil:1715: Using density 891.412528 and viscosity 0.0032285361678330597 of oiltype TROLL, STATOIL
20:13:32 INFO    opendrift.models.basemodel.environment:218: Adding a dynamical landmask with max. priority based on assumed maximum speed of 1.3 m/s. Adding a customised landmask may be faster...
20:13:35 INFO    opendrift.models.basemodel.environment:245: Fallback values will be used for the following variables which have no readers:
20:13:35 INFO    opendrift.models.basemodel.environment:248:    x_sea_water_velocity: -0.200000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    y_sea_water_velocity: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    x_wind: 6.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    y_wind: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_surface_height: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    upward_sea_water_velocity: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_surface_wave_significant_height: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_surface_wave_stokes_drift_x_velocity: 0.300000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_surface_wave_stokes_drift_y_velocity: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_surface_wave_period_at_variance_spectral_density_maximum: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_surface_wave_mean_period_from_variance_spectral_density_second_frequency_moment: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_ice_area_fraction: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_ice_x_velocity: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_ice_y_velocity: 0.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_water_temperature: 5.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_water_salinity: 34.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    sea_floor_depth_below_sea_level: 10000.000000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    ocean_vertical_diffusivity: 0.020000
20:13:35 INFO    opendrift.models.basemodel.environment:248:    ocean_mixed_layer_thickness: 50.000000
20:13:36 INFO    opendrift.models.basemodel:935: Using existing reader for land_binary_mask
20:13:36 INFO    opendrift.models.basemodel:946: All points are in ocean
20:13:36 INFO    opendrift.models.openoil.openoil:684: Oil-water surface tension is 0.032037 Nm
20:13:36 INFO    opendrift.models.openoil.openoil:693: Using max water fractions [0.74, 0.75] for temperatures [5.0, 15.0] for oiltype TROLL, STATOIL
20:13:36 INFO    opendrift.models.openoil.openoil:694: Corresponding max water fraction from GNOME is 0.5262214946654028
20:13:36 INFO    opendrift.models.basemodel:2035: 2024-12-12 20:13:32.521591 - step 1 of 2 - 10000 active elements (0 deactivated)
20:13:36 INFO    opendrift.models.basemodel:2035: 2024-12-12 21:13:32.521591 - 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