Source code for opendrift.models.chemicaldrift

# 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 2020, Manuel Aghito, MET Norway

"""
ChemicalDrift is an OpenDrift module for drift and fate of chemicals.
The module is under development within the scope of the Horizon2020 project EMERGE
Manuel Aghito. Norwegian Meteorological Institute. 2021.
The initial version is based on Radionuclides module by Magne Simonsen
"""

import numpy as np
import logging; logger = logging.getLogger(__name__)

from opendrift.models.oceandrift import OceanDrift, Lagrangian3DArray
from opendrift.config import CONFIG_LEVEL_ESSENTIAL, CONFIG_LEVEL_BASIC, CONFIG_LEVEL_ADVANCED
import pyproj
from datetime import datetime

# Defining the Chemical element properties
[docs] class Chemical(Lagrangian3DArray): """Extending Lagrangian3DArray with specific properties for chemicals """ variables = Lagrangian3DArray.add_variables([ ('diameter', {'dtype': np.float32, 'units': 'm', 'default': 0.}), #('neutral_buoyancy_salinity', {'dtype': np.float32, # 'units': '[]', # 'default': 31.25}), # for NEA Cod ('density', {'dtype': np.float32, 'units': 'kg/m^3', 'default': 2650.}), # Mineral particles ('specie', {'dtype': np.int32, 'units': '', 'default': 0}), ('mass', {'dtype': np.float32, 'units': 'ug', 'seed': True, 'default': 1e3}), ('mass_degraded', {'dtype': np.float32, 'units': 'ug', 'seed': True, 'default': 0}), ('mass_degraded_water', {'dtype': np.float32, 'units': 'ug', 'seed': True, 'default': 0}), ('mass_degraded_sediment', {'dtype': np.float32, 'units': 'ug', 'seed': True, 'default': 0}), ('mass_volatilized', {'dtype': np.float32, 'units': 'ug', 'seed': True, 'default': 0}) ])
[docs] class ChemicalDrift(OceanDrift): """Chemical particle trajectory model based on the OpenDrift framework. Developed at MET Norway Generic module for particles that are subject to vertical turbulent mixing with the possibility for positive or negative buoyancy Particles could be e.g. oil droplets, plankton, or sediments Chemical functionality include interactions with solid matter (particles and sediments) through transformation processes, implemented with stochastic approach for dynamic partitioning. Under construction. """ ElementType = Chemical required_variables = { 'x_sea_water_velocity': {'fallback': None}, 'y_sea_water_velocity': {'fallback': None}, 'sea_surface_height': {'fallback': 0}, 'x_wind': {'fallback': 0}, 'y_wind': {'fallback': 0}, 'land_binary_mask': {'fallback': None}, 'sea_floor_depth_below_sea_level': {'fallback': 10000}, 'ocean_vertical_diffusivity': {'fallback': 0.0001, 'profiles': True}, 'sea_water_temperature': {'fallback': 10, 'profiles': True}, 'sea_water_salinity': {'fallback': 34, 'profiles': True}, 'upward_sea_water_velocity': {'fallback': 0}, #'conc3': {'fallback': 1.e-3}, 'spm': {'fallback': 1}, 'ocean_mixed_layer_thickness': {'fallback': 50}, 'active_sediment_layer_thickness': {'fallback': 0.03}, # TODO - currently not used, redundant with 'chemical:sediment:mixing_depth' 'doc': {'fallback': 0.0}, # Variables for dissociation 'sea_water_ph_reported_on_total_scale':{'fallback': 8.1, 'profiles': True}, # water_pH from CMENS with standard name # 'pH_sediment':{'fallback': 6.9, 'profiles': False}, # supplied by the user, with pH_sediment as standard name # }
[docs] def specie_num2name(self,num): return self.name_species[num]
[docs] def specie_name2num(self,name): num = self.name_species.index(name) return num
def __init__(self, *args, **kwargs): # Calling general constructor of parent class super(ChemicalDrift, self).__init__(*args, **kwargs) # TODO: descriptions and units must be added in config setting below self._add_config({ 'chemical:transfer_setup': {'type': 'enum', 'enum': ['Sandnesfj_Al','metals', '137Cs_rev', 'custom', 'organics'], 'default': 'custom', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:dynamic_partitioning': {'type': 'bool', 'default': True, 'level': CONFIG_LEVEL_BASIC, 'description': 'Toggle dynamic partitioning'}, 'chemical:slowly_fraction': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:irreversible_fraction': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:dissolved_diameter': {'type': 'float', 'default': 0, 'min': 0, 'max': 100e-6, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:particle_diameter': {'type': 'float', 'default': 5e-6, 'min': 0, 'max': 100e-6, 'units': 'm', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:particle_concentration_half_depth': {'type': 'float', 'default': 20, 'min': 0, 'max': 100, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:doc_concentration_half_depth': {'type': 'float', 'default': 1000, # TODO: check better 'min': 0, 'max': 1000, 'units': 'm', # Vertical conc drops more slowly slower than for SPM 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, # example: 10.3389/fmars.2017.00436. lower limit around 40 umol/L 'chemical:particle_diameter_uncertainty': {'type': 'float', 'default': 1e-7, 'min': 0, 'max': 100e-6, 'units': 'm', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'seed:LMM_fraction': {'type': 'float','default': .1, 'min': 0, 'max': 1, 'units': '', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'seed:particle_fraction': {'type': 'float','default': 0.9, 'min': 0, 'max': 1, 'units': '', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, # Species 'chemical:species:LMM': {'type': 'bool', 'default': True, 'level': CONFIG_LEVEL_BASIC, 'description': 'Toggle LMM species'}, 'chemical:species:LMMcation': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_BASIC, 'description': ''}, 'chemical:species:LMManion': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_BASIC, 'description': ''}, 'chemical:species:Colloid': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_BASIC, 'description': ''}, 'chemical:species:Humic_colloid': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:species:Polymer': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:species:Particle_reversible': {'type': 'bool', 'default': True, 'level': CONFIG_LEVEL_BASIC, 'description': ''}, 'chemical:species:Particle_slowly_reversible': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_BASIC, 'description': ''}, 'chemical:species:Particle_irreversible': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:species:Sediment_reversible': {'type': 'bool', 'default': True, 'level': CONFIG_LEVEL_BASIC, 'description': ''}, 'chemical:species:Sediment_slowly_reversible': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_BASIC, 'description': ''}, 'chemical:species:Sediment_irreversible': {'type': 'bool', 'default': False, 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, # Transformations 'chemical:transformations:Kd': {'type': 'float', 'default': 2.0, 'min': 0, 'max': 1e9, 'units': 'm3/kg', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:S0': {'type': 'float', 'default': 0.0, 'min': 0, 'max': 100, 'units': 'PSU', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'parameter controlling salinity dependency of Kd for metals'}, 'chemical:transformations:Dc': {'type': 'float', 'default': 1.16e-5, # Simonsen 2019 'min': 0, 'max': 1e6, 'units': '', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:transformations:slow_coeff': {'type': 'float', 'default': 0, #1.2e-7, # Simonsen 2019 'min': 0, 'max': 1e6, 'units': '', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:transformations:volatilization': {'type': 'bool', 'default': False, 'description': 'Chemical is evaporated.', 'level': CONFIG_LEVEL_BASIC}, 'chemical:transformations:degradation': {'type': 'bool', 'default': False, 'description': 'Chemical mass is degraded.', 'level': CONFIG_LEVEL_BASIC}, 'chemical:transformations:degradation_mode': {'type': 'enum', 'enum': ['OverallRateConstants'], 'default': 'OverallRateConstants', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, # sorption/desorption 'chemical:transformations:dissociation': {'type': 'enum', 'enum': ['nondiss','acid', 'base', 'amphoter'], 'default': 'nondiss', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:LogKOW': {'type': 'float', 'default': 3.361, # Naphthalene 'min': -3, 'max': 10, 'units': 'Log L/Kg', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:TrefKOW': {'type': 'float', 'default': 25., # Naphthalene 'min': -3, 'max': 30, 'units': 'C', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:DeltaH_KOC_Sed': {'type': 'float', 'default': -21036., # Naphthalene 'min': -100000., 'max': 100000., 'units': 'J/mol', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:DeltaH_KOC_DOM': {'type': 'float', 'default': -25900., # Naphthalene 'min': -100000., 'max': 100000., 'units': 'J/mol', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:Setchenow': {'type': 'float', 'default': 0.2503, # Naphthalene 'min': 0, 'max': 1, 'units': 'L/mol', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:pKa_acid': {'type': 'float', 'default': -1, 'min': -1, 'max': 14, 'units': '', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:transformations:pKa_base': {'type': 'float', 'default': -1, 'min': -1, 'max': 14, 'units': '', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:transformations:KOC_DOM': {'type': 'float', 'default': -1, 'min': -1, 'max': 10000000000, 'units': 'L/KgOC', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:transformations:KOC_sed': {'type': 'float', 'default': -1, 'min': -1, 'max': 10000000000, 'units': 'L/KgOC', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:transformations:fOC_SPM': {'type': 'float', 'default': 0.05, 'min': 0.01, 'max': 0.1, 'units': 'gOC/g', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:transformations:fOC_sed': {'type': 'float', 'default': 0.05, 'min': 0.01, 'max': 0.1, 'units': 'gOC/g', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:transformations:aggregation_rate': {'type': 'float', 'default': 0, 'min': 0, 'max': 1, 'units': 's-1', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, # Degradation in water column 'chemical:transformations:t12_W_tot': {'type': 'float', 'default': 224.08, # Naphthalene 'min': 1, 'max': None, 'units': 'hours', 'level': CONFIG_LEVEL_ADVANCED, 'description': 'half life in water, total'}, 'chemical:transformations:Tref_kWt': {'type': 'float', 'default': 25., # Naphthalene 'min': -3, 'max': 30, 'units': 'C', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:DeltaH_kWt': {'type': 'float', 'default': 50000., # generic 'min': -100000., 'max': 100000., 'units': 'J/mol', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, # Degradation in sediment layer 'chemical:transformations:t12_S_tot': {'type': 'float', 'default': 5012.4, # Naphthalene 'min': 1, 'max': None, 'units': 'hours', 'level': CONFIG_LEVEL_ADVANCED, 'description': 'half life in sediments, total'}, 'chemical:transformations:Tref_kSt': {'type': 'float', 'default': 25., # Naphthalene 'min': -3, 'max': 30, 'units': 'C', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, 'chemical:transformations:DeltaH_kSt': {'type': 'float', 'default': 50000., # generic 'min': -100000., 'max': 100000., 'units': 'J/mol', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, # Volatilization 'chemical:transformations:MolWt': {'type': 'float', 'default': 128.1705, # Naphthalene 'min': 50, 'max': 1000, 'units': 'amu', 'level': CONFIG_LEVEL_ADVANCED, 'description': 'molecular weight'}, 'chemical:transformations:Henry': {'type': 'float', 'default': 4.551e-4, # Napththalene 'min': None, 'max': None, 'units': 'atm m3 mol-1', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Henry constant'}, # vapour pressure 'chemical:transformations:Vpress': {'type': 'float', 'default': 11.2, # Naphthalene 'min': None, 'max': None, 'units': 'Pa', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Vapour pressure'}, 'chemical:transformations:Tref_Vpress': {'type': 'float', 'default': 25., # Naphthalene 'min': None, 'max': None, 'units': 'C', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Vapour pressure ref temp'}, 'chemical:transformations:DeltaH_Vpress': {'type': 'float', 'default': 55925., # Naphthalene 'min': -100000., 'max': 115000., 'units': 'J/mol', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Enthalpy of volatilization'}, # solubility 'chemical:transformations:Solub': {'type': 'float', 'default': 31.4, # Naphthalene 'min': None, 'max': None, 'units': 'g/m3', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Solubility'}, 'chemical:transformations:Tref_Solub': {'type': 'float', 'default': 25., # Naphthalene 'min': None, 'max': None, 'units': 'C', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Solubility ref temp'}, 'chemical:transformations:DeltaH_Solub': {'type': 'float', 'default': 25300., # Naphthalene 'min': -100000., 'max': 100000., 'units': 'J/mol', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Enthalpy of solubilization'}, # Sedimentation/Resuspension 'chemical:sediment:mixing_depth': {'type': 'float', 'default': 0.03, 'min': 0, 'max': 100, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:density': {'type': 'float', 'default': 2600, 'min': 0, 'max': 10000, 'units': 'kg/m3', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:effective_fraction': {'type': 'float', 'default': 0.9, 'min': 0, 'max': 1, 'units': '', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:corr_factor': {'type': 'float', 'default': 0.1, 'min': 0, 'max': 10, 'units': '', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:porosity': {'type': 'float', 'default': 0.6, 'min': 0, 'max': 1, 'units': '', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:layer_thickness': {'type': 'float', 'default': 1, 'min': 0, 'max': 100, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:desorption_depth': {'type': 'float', 'default': 1, 'min': 0, 'max': 100, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:desorption_depth_uncert': {'type': 'float', 'default': .5, 'min': 0, 'max': 100, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:resuspension_depth': {'type': 'float', 'default': 1, 'min': 0, 'max': 100, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:resuspension_depth_uncert': {'type': 'float', 'default': .5, 'min': 0, 'max': 100, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:resuspension_critvel': {'type': 'float', 'default': .01, 'min': 0, 'max': 1, 'units': 'm/s', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:burial_rate': {'type': 'float', 'default': .00003, # MacKay 'min': 0, 'max': 10, 'units': 'm/year', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'chemical:sediment:buried_leaking_rate': {'type': 'float', 'default': 0, 'min': 0, 'max': 10, 'units': 's-1', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, # 'chemical:compound': {'type': 'enum', 'enum': ['Naphthalene','Phenanthrene','Fluoranthene', 'Benzo-a-anthracene','Benzo-a-pyrene','Dibenzo-ah-anthracene', 'C1-Naphthalene','Acenaphthene','Acenaphthylene','Fluorene', 'Dibenzothiophene','C2-Naphthalene','Anthracene','C3-Naphthalene','C1-Dibenzothiophene', 'Pyrene','C1-Phenanthrene','C2-Dibenzothiophene', 'C2-Phenanthrene','Benzo-b-fluoranthene','Chrysene', 'C3-Dibenzothiophene','C3-Phenanthrene', 'Benzo-k-fluoranthene','Benzo-ghi-perylene','Indeno-123cd-pyrene', 'Copper','Cadmium','Chromium','Lead','Vanadium','Zinc','Nickel',None], 'default': None, 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, })
[docs] def prepare_run(self): logger.info( 'Number of species: {}'.format(self.nspecies) ) for i,sp in enumerate(self.name_species): logger.info( '{:>3} {}'.format( i, sp ) ) logger.info( 'transfer setup: %s' % self.get_config('chemical:transfer_setup')) logger.info('nspecies: %s' % self.nspecies) logger.info('Transfer rates:\n %s' % self.transfer_rates) self.SPM_vertical_levels_given = False for key, value in self.env.readers.items(): if 'spm' in value.variables: if (hasattr(value,'sigma') or hasattr(value,'z') ): self.SPM_vertical_levels_given = True self.DOC_vertical_levels_given = False for key, value in self.env.readers.items(): if 'doc' in value.variables: if (hasattr(value,'sigma') or hasattr(value,'z') ): self.DOC_vertical_levels_given = True super(ChemicalDrift, self).prepare_run()
[docs] def init_species(self): # Initialize specie types if self.get_config('chemical:transfer_setup')=='metals': self.set_config('chemical:species:LMM',True) self.set_config('chemical:species:Particle_reversible', True) self.set_config('chemical:species:Particle_slowly_reversible', True) self.set_config('chemical:species:Sediment_reversible', True) self.set_config('chemical:species:Sediment_slowly_reversible', True) elif self.get_config('chemical:transfer_setup')=='137Cs_rev': self.set_config('chemical:species:LMM',True) self.set_config('chemical:species:Particle_reversible', True) self.set_config('chemical:species:Sediment_reversible', True) elif self.get_config('chemical:transfer_setup')=='Sandnesfj_Al': self.set_config('chemical:species:LMM', False) self.set_config('chemical:species:LMMcation', True) self.set_config('chemical:species:LMManion', True) self.set_config('chemical:species:Humic_colloid', True) self.set_config('chemical:species:Polymer', True) self.set_config('chemical:species:Particle_reversible', True) self.set_config('chemical:species:Sediment_reversible', True) elif self.get_config('chemical:transfer_setup')=='organics': self.set_config('chemical:species:LMM',True) self.set_config('chemical:species:Particle_reversible', True) self.set_config('chemical:species:Particle_slowly_reversible', False) self.set_config('chemical:species:Sediment_reversible', True) self.set_config('chemical:species:Sediment_slowly_reversible', True) self.set_config('chemical:species:Humic_colloid', True) elif self.get_config('chemical:transfer_setup')=='custom': # Do nothing, species must be set manually pass else: logger.error('No valid transfer_setup {}'.format(self.get_config('chemical:transfer_setup'))) self.name_species=[] if self.get_config('chemical:species:LMM'): self.name_species.append('LMM') if self.get_config('chemical:species:LMMcation'): self.name_species.append('LMMcation') if self.get_config('chemical:species:LMManion'): self.name_species.append('LMManion') if self.get_config('chemical:species:Colloid'): self.name_species.append('Colloid') if self.get_config('chemical:species:Humic_colloid'): self.name_species.append('Humic colloid') if self.get_config('chemical:species:Polymer'): self.name_species.append('Polymer') if self.get_config('chemical:species:Particle_reversible'): self.name_species.append('Particle reversible') if self.get_config('chemical:species:Particle_slowly_reversible'): self.name_species.append('Particle slowly reversible') if self.get_config('chemical:species:Particle_irreversible'): self.name_species.append('Particle irreversible') if self.get_config('chemical:species:Sediment_reversible'): self.name_species.append('Sediment reversible') if self.get_config('chemical:species:Sediment_slowly_reversible'): self.name_species.append('Sediment slowly reversible') if self.get_config('chemical:species:Sediment_irreversible'): self.name_species.append('Sediment irreversible') if self.get_config('chemical:species:Sediment_slowly_reversible') and \ self.get_config('chemical:species:Particle_slowly_reversible'): self.set_config('chemical:slowly_fraction', True) if self.get_config('chemical:species:Sediment_irreversible') and \ self.get_config('chemical:species:Particle_irreversible'): self.set_config('chemical:irreversible_fraction', True) self.nspecies = len(self.name_species)
# logger.info( 'Number of species: {}'.format(self.nspecies) ) # for i,sp in enumerate(self.name_species): # logger.info( '{:>3} {}'.format( i, sp ) )
[docs] def seed_elements(self, *args, **kwargs): if hasattr(self,'name_species') == False: self.init_species() self.init_transfer_rates() if 'number' in kwargs: num_elements = kwargs['number'] else: num_elements = self.get_config('seed:number') if 'specie' in kwargs: print('num_elements', num_elements) try: print('len specie:',len(kwargs['specie'])) except: print('specie:',kwargs['specie']) init_specie = np.ones(num_elements,dtype=int) init_specie[:] = kwargs['specie'] else: # Set initial partitioning if 'particle_fraction' in kwargs: particle_frac = kwargs['particle_fraction'] else: particle_frac = self.get_config('seed:particle_fraction') if 'LMM_fraction' in kwargs: lmm_frac = kwargs['LMM_fraction'] else: lmm_frac = self.get_config('seed:LMM_fraction') if not lmm_frac + particle_frac == 1.: logger.error('Fraction does not sum up to 1: %s' % str(lmm_frac+particle_frac) ) logger.error('LMM fraction: %s ' % str(lmm_frac)) logger.error( 'Particle fraction %s '% str(particle_frac) ) raise ValueError('Illegal specie fraction combination : ' + str(lmm_frac) + ' '+ str(particle_frac) ) init_specie = np.ones(num_elements, int) dissolved=np.random.rand(num_elements)<lmm_frac if self.get_config('chemical:transfer_setup')=='Sandnesfj_Al': init_specie[dissolved]=self.num_lmmcation else: init_specie[dissolved]=self.num_lmm init_specie[~dissolved]=self.num_prev kwargs['specie'] = init_specie logger.debug('Initial partitioning:') for i,sp in enumerate(self.name_species): logger.debug( '{:>9} {:>3} {:24} '.format( np.sum(init_specie==i), i, sp ) ) # Set initial particle size if 'diameter' in kwargs: diameter = kwargs['diameter'] else: diameter = self.get_config('chemical:particle_diameter') std = self.get_config('chemical:particle_diameter_uncertainty') init_diam = np.zeros(num_elements,float) init_diam[init_specie==self.num_prev] = diameter + np.random.normal(0, std, sum(init_specie==self.num_prev)) kwargs['diameter'] = init_diam super(ChemicalDrift, self).seed_elements(*args, **kwargs)
[docs] def tempcorr(self,mode,DeltaH,T_C,Tref_C): ''' Temperature correction using Arrhenius or Q10 method ''' if mode == 'Arrhenius': R = 8.3145 # J/(mol*K) T_K = T_C + 273.15 Tref_K = Tref_C + 273.15 corr = np.e**(-(DeltaH/R)*(1/T_K - 1/Tref_K)) elif mode =='Q10': corr = 2**((T_C - Tref_C)/10) return corr
[docs] def salinitycorr(self,Setschenow,Temperature,Salinity): ''' Salinity correction ''' # Setschenow constant for the given chemical (L/mol) # Salinity (PSU =g/Kg) # Temperature (Celsius) MWsalt = 68.35 # average mass of sea water salt (g/mol) Schwarzenbach Gschwend Imboden Environmental Organic Chemistry Dens_sw = self.sea_water_density(T=Temperature, S=Salinity)*1e-3 # (Kg/L) # ConcSalt= (Salinitypsu/MWsalt)∙Dens_sw # = ( g/Kg / g/mol )∙ Kg/L # = mol/Kg ∙ Kg/L = mol/L ConcSalt = (Salinity/MWsalt)*Dens_sw # Log(Kd_fin)=(Setschenow ∙ ConcSalt)+Log(Kd_T) # Kd_fin = 10^(Setschenow ∙ ConcSalt) * Kd_T corr = 10**(Setschenow*ConcSalt) return corr
### Functions to update partitioning coefficients
[docs] def calc_KOC_sedcorr(self, KOC_sed_initial, KOC_sed_n, pKa_acid, pKa_base, KOW, pH_sed, diss): ''' Calculate correction of KOC due to pH of sediments ''' # Estimate KOC for dissociated forms from KOW KOC_sed_diss_acid = (10**(0.11*np.log10(KOW)+1.54)) # KOC for dissociated acid species (L/kg_OC), from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202. KOC_sed_diss_base = 10**(pKa_acid**(0.65*((KOW/(KOW+1))**0.14))) # KOC for ionized form of base species (L/kg_OC) # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 # Updated values of KOC to calculate correction factor KOC_sed_updated = np.empty_like(pH_sed) KOC_sedcorr= np.empty_like(pH_sed) for i in (range(len(pH_sed))): if diss=='acid': Phi_n_sed = 1/(1 + 10**(pH_sed[i]-pKa_acid)) Phi_diss_sed = 1-Phi_n_sed KOC_sed_updated[i] = (KOC_sed_n * Phi_n_sed) + (Phi_diss_sed * KOC_sed_diss_acid) KOC_sedcorr[i] = KOC_sed_updated[i]/KOC_sed_initial elif diss=='base': Phi_n_sed = 1/(1 + 10**(pH_sed[i]-pKa_base)) Phi_diss_sed = 1-Phi_n_sed KOC_sed_updated[i] = (KOC_sed_n * Phi_n_sed) + (Phi_diss_sed * KOC_sed_diss_acid) KOC_sedcorr[i] = KOC_sed_updated[i]/KOC_sed_initial elif diss=='amphoter': Phi_n_sed = 1/(1 + 10**(pH_sed[i]-pKa_acid) + 10**(pKa_base)) Phi_anion_sed = Phi_n_sed * 10**(pH_sed[i]-pKa_acid) Phi_cation_sed = Phi_n_sed * 10**(pKa_base-pH_sed[i]) KOC_sed_updated[i] = (KOC_sed_n * Phi_n_sed) + (Phi_anion_sed * KOC_sed_diss_acid) + (Phi_cation_sed * KOC_sed_diss_base) KOC_sedcorr[i]=KOC_sed_updated[i]/KOC_sed_initial elif diss=='undiss': KOC_sedcorr[i]=1 return KOC_sedcorr
[docs] def calc_KOC_watcorrSPM(self, KOC_SPM_initial, KOC_sed_n, pKa_acid, pKa_base, KOW, pH_water_SPM, diss): ''' Calculate correction of KOC due to pH of water for SPM ''' # Estimate KOC for dissociated forms from KOW KOC_sed_diss_acid = (10**(0.11*np.log10(KOW)+1.54)) # KOC for dissociated acid species (L/kg_OC), from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202. KOC_sed_diss_base = 10**(pKa_acid**(0.65*((KOW/(KOW+1))**0.14))) # KOC for ionized form of base species (L/kg_OC) # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 KOC_SPM_updated = np.empty_like(pH_water_SPM) KOC_SPMcorr = np.empty_like(pH_water_SPM) for i in (range(len(pH_water_SPM))): if diss=='acid': Phi_n_SPM = 1/(1 + 10**(pH_water_SPM[i]-pKa_acid)) Phi_diss_SPM = 1-Phi_n_SPM KOC_SPM_updated[i] = (KOC_sed_n * Phi_n_SPM) + (Phi_diss_SPM * KOC_sed_diss_acid) KOC_SPMcorr[i]=KOC_SPM_updated[i]/KOC_SPM_initial elif diss=='base': Phi_n_SPM = 1/(1 + 10**(pH_water_SPM[i]-pKa_base)) Phi_diss_SPM = 1-Phi_n_SPM KOC_SPM_updated[i] = (KOC_sed_n * Phi_n_SPM) + (Phi_diss_SPM * KOC_sed_diss_acid) KOC_SPMcorr[i]=KOC_SPM_updated[i]/KOC_SPM_initial elif diss=='amphoter': Phi_n_SPM = 1/(1 + 10**(pH_water_SPM[i]-pKa_acid) + 10**(pKa_base)) Phi_anion_SPM = Phi_n_SPM * 10**(pH_water_SPM[i]-pKa_acid) Phi_cation_SPM = Phi_n_SPM * 10**(pKa_base-pH_water_SPM[i]) KOC_SPM_updated[i] = (KOC_sed_n * Phi_n_SPM) + (Phi_anion_SPM * KOC_sed_diss_acid) + (Phi_cation_SPM * KOC_sed_diss_base) KOC_SPMcorr[i]=KOC_SPM_updated[i]/KOC_SPM_initial elif diss=='undiss': KOC_SPMcorr[i]=1 return KOC_SPMcorr
[docs] def calc_KOC_watcorrDOM(self, KOC_DOM_initial, KOC_DOM_n, pKa_acid, pKa_base, KOW, pH_water_DOM, diss): ''' Calculate correction of KOC due to pH of water for DOM ''' Phi_n_DOM = np.empty_like(pH_water_DOM) Phi_diss_DOM = np.empty_like(pH_water_DOM) KOC_DOM_updated = np.empty_like(pH_water_DOM) KOC_DOMcorr = np.empty_like(pH_water_DOM) for i in (range(len(pH_water_DOM))): if diss=='acid': Phi_n_DOM = 1/(1 + 10**(pH_water_DOM[i]-pKa_acid)) Phi_diss_DOM = 1-Phi_n_DOM KOC_DOM_updated[i] = (0.08 * ((Phi_n_DOM*(KOC_DOM_n)) + ((1 - Phi_diss_DOM)*10**(np.log10(KOW)-3.5))))/0.526 # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 KOC_DOMcorr[i]=KOC_DOM_updated[i]/KOC_DOM_initial elif diss=='base': Phi_n_DOM = 1/(1 + 10**(pH_water_DOM[i]-pKa_base)) Phi_diss_DOM = 1-Phi_n_DOM KOC_DOM_updated[i] = (0.08 * ((Phi_n_DOM*(KOC_DOM_n)) + ((1 - Phi_diss_DOM)*10**(np.log10(KOW)-3.5))))/0.526 # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 KOC_DOMcorr[i]=KOC_DOM_updated[i]/KOC_DOM_initial elif diss=='amphoter': Phi_n_DOM = 1/(1 + 10**(pH_water_DOM[i]-pKa_acid) + 10**(pKa_base)) Phi_diss_DOM = 1-Phi_n_DOM KOC_DOM_updated[i] = (0.08 * ((Phi_n_DOM*(KOC_DOM_n)) + ((1 - Phi_diss_DOM)*10**(np.log10(KOW)-3.5))))/0.526 # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 KOC_DOMcorr[i]=KOC_DOM_updated[i]/KOC_DOM_initial elif diss=='undiss': KOC_DOMcorr[i]=1 return KOC_DOMcorr
[docs] def init_transfer_rates(self): ''' Initialization of background values in the transfer rates 2D array. ''' transfer_setup=self.get_config('chemical:transfer_setup') # logger.info( 'transfer setup: %s' % transfer_setup) self.transfer_rates = np.zeros([self.nspecies,self.nspecies]) self.ntransformations = np.zeros([self.nspecies,self.nspecies]) if transfer_setup == 'organics': self.num_lmm = self.specie_name2num('LMM') self.num_humcol = self.specie_name2num('Humic colloid') self.num_prev = self.specie_name2num('Particle reversible') self.num_srev = self.specie_name2num('Sediment reversible') #self.num_psrev = self.specie_name2num('Particle slowly reversible') self.num_ssrev = self.specie_name2num('Sediment slowly reversible') # Values from EMERGE-Aquatox Org2C = 0.526 # kgOC/KgOM #Kd = self.get_config('chemical:transformations:Kd') KOW = 10**self.get_config('chemical:transformations:LogKOW') KOWTref = self.get_config('chemical:transformations:TrefKOW') DH_KOC_Sed = self.get_config('chemical:transformations:DeltaH_KOC_Sed') DH_KOC_DOM = self.get_config('chemical:transformations:DeltaH_KOC_DOM') Setchenow = self.get_config('chemical:transformations:Setchenow') diss = self.get_config('chemical:transformations:dissociation') pKa_acid = self.get_config('chemical:transformations:pKa_acid') if pKa_acid < 0 and diss!='nondiss': raise ValueError("pKa_acid must be positive") # print("pKa_acid must be positive") # UserWarning(("pKa_acid must be positive")) else: pass pKa_base = self.get_config('chemical:transformations:pKa_base') if pKa_base < 0 and diss!='nondiss': raise ValueError("pKa_base must be positive") # print("pKa_base must be positive") # UserWarning(("pKa_base must be positive")) else: pass if diss == 'amphoter' and abs(pKa_acid - pKa_base) < 2: raise ValueError("pKa_base and pKa_acid must differ of at least two units") else: pass # Read water pH to calculate dissociation # pH_water = self.environment.sea_water_ph_reported_on_total_scale pH_water = 8.1 # 8.1 pH_sed = 6.9 # 6.9 fOC_SPM = self.get_config('chemical:transformations:fOC_SPM') # typical values from 0.01 to 0.1 gOC/g fOC_sed = self.get_config('chemical:transformations:fOC_sed') # typical values from 0.01 to 0.1 gOC/g concDOM = 1.e-3 / Org2C # concentration of available dissolved organic matter (kg/m3) # rough initial estimate for coastal waters, doi: 10.1002/lom3.10118 #concDOM = 50.e-3 # HIGHER VALUE FOR TESTING!!!!!!!!!!!! # Values from Simonsen et al (2019a) slow_coeff = self.get_config('chemical:transformations:slow_coeff') concSPM = 50.e-3 # available SPM (kg/m3) sed_L = self.get_config('chemical:sediment:mixing_depth') # sediment mixing depth (m) sed_dens = self.get_config('chemical:sediment:density') # default particle density (kg/m3) sed_phi = self.get_config('chemical:sediment:corr_factor') # sediment correction factor sed_poro = self.get_config('chemical:sediment:porosity') # sediment porosity sed_H = self.get_config('chemical:sediment:layer_thickness') # thickness of seabed interaction layer (m) sed_burial = self.get_config('chemical:sediment:burial_rate') # sediment burial rate (m/y) sed_leaking_rate = self.get_config( 'chemical:sediment:buried_leaking_rate') if diss=='nondiss': KOC_DOM = self.get_config('chemical:transformations:KOC_DOM') if KOC_DOM < 0: KOC_DOM = 2.88 * KOW**0.67 # (L/KgOC), Park and Clough, 2014 KOC_sed = self.get_config('chemical:transformations:KOC_sed') if KOC_sed < 0: KOC_sed = 2.62 * KOW**0.82 # (L/KgOC), Park and Clough, 2014 (334)/Org2C #KOC_Sed = 1.26 * kOW**0.81 # (L/KgOC),Ragas et al., 2019 KOC_SPM = KOC_sed else: if diss=='acid': # Dissociation in water Phi_n_water = 1/(1 + 10**(pH_water-pKa_acid)) Phi_diss_water = 1-Phi_n_water KOC_sed_n = self.get_config('chemical:transformations:KOC_sed') if KOC_sed_n < 0: KOC_sed_n = 2.62 * KOW**0.82 # (L/KgOC), Park and Clough, 2014 (334)/Org2C TO DO Add if choice between input and estimation else: pass KOC_sed_diss_acid = (10**(0.11*np.log10(KOW)+1.54)) # KOC for dissociated acid species (L/kg_OC), from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202. KOC_SPM = (KOC_sed_n * Phi_n_water) + (Phi_diss_water * KOC_sed_diss_acid) KOC_DOM_n = self.get_config('chemical:transformations:KOC_DOM') if KOC_DOM_n <0: KOC_DOM_n = 2.88 * KOW**0.67 # (L/KgOC), Park and Clough, 2014 TO DO Add if choice between input and estimation else: pass KOC_DOM = (0.08 * ((Phi_n_water*(KOC_DOM_n)) + ((1 - Phi_diss_water)*10**(np.log10(KOW)-3.5))))/0.526 # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 # Dissociation in sediments Phi_n_sed = 1/(1 + 10**(pH_sed-pKa_acid)) Phi_diss_sed = 1-Phi_n_sed KOC_sed = (KOC_sed_n * Phi_n_sed) + (Phi_diss_sed * KOC_sed_diss_acid) elif diss=='base': # Dissociation in water Phi_n_water = 1/(1 + 10**(pH_water-pKa_base)) Phi_diss_water = 1-Phi_n_water KOC_sed_n = self.get_config('chemical:transformations:KOC_sed') if KOC_sed_n <0: KOC_sed_n = 2.62 * KOW**0.82 # (L/KgOC), Park and Clough, 2014 (334)/Org2C TO DO Add if choice between input and estimation else: pass KOC_sed_diss_base = 10**(pKa_acid**(0.65*((KOW/(KOW+1))**0.14))) # KOC for ionized form of base species (L/kg_OC) # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 KOC_SPM = (KOC_sed_n * Phi_n_water) + (Phi_diss_water * KOC_sed_diss_base) KOC_DOM_n = self.get_config('chemical:transformations:KOC_DOM') if KOC_DOM_n <0: KOC_DOM_n = 2.88 * KOW**0.67 # (L/KgOC), Park and Clough, 2014 TO DO Add if choice between input and estimation else: pass KOC_DOM = (0.08 * ((Phi_n_water*(KOC_DOM_n)) + ((1 - Phi_diss_water)*10**(np.log10(KOW)-3.5))))/0.526 # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 # Dissociation in sediments Phi_n_sed = 1/(1 + 10**(pH_sed-pKa_base)) Phi_diss_sed = 1-Phi_n_sed KOC_sed = (KOC_sed_n * Phi_n_sed) + (Phi_diss_sed * KOC_sed_diss_base) elif diss=='amphoter': # Dissociation in water # This approach ignores the zwitterionic fraction. 10.1002/etc.115 Phi_n_water = 1/(1 + 10**(pH_water-pKa_acid) + 10**(pKa_base)) Phi_anion_water = Phi_n_water * 10**(pH_water-pKa_acid) Phi_cation_water = Phi_n_water * 10**(pKa_base-pH_water) Phi_diss_water = 1 - Phi_n_water KOC_sed_n = self.get_config('chemical:transformations:KOC_sed') if KOC_sed_n < 0: KOC_sed_n = KOC_sed_n = 2.62 * KOW**0.82 # (L/KgOC), Park and Clough, 2014 (334)/Org2C TO DO Add if choice between input and estimation else: pass KOC_sed_diss_base = 10**(pKa_acid**(0.65*((KOW/(KOW+1))**0.14))) # KOC for ionized form of base species (L/kg_OC) # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 KOC_sed_diss_acid = (10**(0.11*np.log10(KOW)+1.54)) # KOC for dissociated acid species (L/kg_OC), from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202. KOC_SPM = (KOC_sed_n * Phi_n_water) + (Phi_anion_water * KOC_sed_diss_acid) + (Phi_cation_water * KOC_sed_diss_base) KOC_DOM_n = self.get_config('chemical:transformations:KOC_DOM') if KOC_DOM_n <0: KOC_DOM_n = 2.88 * KOW**0.67 # (L/KgOC), Park and Clough, 2014 TO DO Add if choice between input and estimation else: pass KOC_DOM = (0.08 * ((Phi_n_water*(KOC_DOM_n)) + ((1 - Phi_diss_water)*10**(np.log10(KOW)-3.5))))/0.526 # from http://i-pie.org/wp-content/uploads/2019/12/ePiE_Technical_Manual-Final_Version_20191202 # Dissociation in sediments Phi_n_sed = 1/(1 + 10**(pH_sed-pKa_acid) + 10**(pKa_base)) Phi_anion_sed = Phi_n_sed * 10**(pH_sed-pKa_acid) Phi_cation_sed = Phi_n_sed * 10**(pKa_base-pH_sed) KOC_sed = (KOC_sed_n * Phi_n_sed) + (Phi_anion_sed * KOC_sed_diss_acid) + (Phi_cation_sed * KOC_sed_diss_base) logger.debug('Partitioning coefficients (Tref,freshwater)') logger.debug('KOC_sed: %s L/KgOC' % KOC_sed) logger.debug('KOC_SPM: %s L/KgOC' % KOC_SPM) logger.debug('KOC_DOM: %s L/KgOC' % KOC_DOM) #KOM_sed = KOC_sed * Org2C # L/KgOC * KgOC/KgOM = L/KgOM #KOM_SPM = KOC_sed * Org2C # L/KgOC * KgOC/KgOM = L/KgOM #KOM_DOM = KOC_DOM * Org2C # L/KgOC * KgOC/KgOM = L/KgOM # to be calculated separately for sed, SPM, dom (different KOC, pH, fOC) self.Kd_sed = Kd_sed = KOC_sed * fOC_sed # L/KgOC * KgOC/KG = L/Kg self.Kd_SPM = Kd_SPM = KOC_SPM * fOC_SPM # L/KgOC * KgOC/KG = L/Kg self.Kd_DOM = Kd_DOM = KOC_DOM * Org2C # L/KgOC * KgOC/KgOM * 1KgOM/Kg = L/Kg (=KOM_DOM) # TODO Use setconfig() to store these? logger.debug('Kd_sed: %s L/Kg' % Kd_sed) logger.debug('Kd_SPM: %s L/Kg' % Kd_SPM) logger.debug('Kd_DOM: %s L/Kg' % Kd_DOM) # From Karickhoff and Morris 1985 k_ads = 33.3 / (60*60) # L/(Kg*s) = 33 L/(kgOM*h) k_des_sed = k_ads / Kd_sed # 1/s k_des_SPM = k_ads / Kd_SPM # 1/s k_des_DOM = k_ads / Kd_DOM # 1/s # Default corrections, assuming temperature 25 salinity 35 TcorrSed = self.tempcorr("Arrhenius",DH_KOC_Sed,25,KOWTref) TcorrDOM = self.tempcorr("Arrhenius",DH_KOC_DOM,25,KOWTref) Scorr = self.salinitycorr(Setchenow,KOWTref,35) concSPM = concSPM * 1e-3 # (Kg/L) concDOM = concDOM * 1e-3 # (Kg/L) self.k_ads = k_ads self.k21_0 = k_des_DOM self.k31_0 = k_des_SPM self.k41_0 = k_des_sed * sed_phi # TODO Use setconfig() to store these? self.transfer_rates[self.num_lmm,self.num_humcol] = k_ads * concDOM # k12 self.transfer_rates[self.num_humcol,self.num_lmm] = k_des_DOM / TcorrDOM / Scorr# k21 self.transfer_rates[self.num_lmm,self.num_prev] = k_ads * concSPM # k13 self.transfer_rates[self.num_prev,self.num_lmm] = k_des_SPM / TcorrSed / Scorr # k31 self.transfer_rates[self.num_lmm,self.num_srev] = \ k_ads * sed_L * sed_dens * (1.-sed_poro) * sed_phi / sed_H # k14 # TODO CHECK DIMENSIONS!!!!! L-m3 !!!! self.transfer_rates[self.num_srev,self.num_lmm] = \ k_des_sed * sed_phi / TcorrSed / Scorr # k41 #self.transfer_rates[self.num_srev,self.num_ssrev] = slow_coeff # k46 #self.transfer_rates[self.num_ssrev,self.num_srev] = slow_coeff*.1 # k64 # Using slowly reversible specie for burial - TODO buried sediment should be a new specie self.transfer_rates[self.num_srev,self.num_ssrev] = sed_burial / sed_L / 31556926 # k46 (m/y) / m / (s/y) = s-1 self.transfer_rates[self.num_ssrev,self.num_srev] = sed_leaking_rate # k64 self.transfer_rates[self.num_humcol,self.num_prev] = self.get_config('chemical:transformations:aggregation_rate') self.transfer_rates[self.num_prev,self.num_humcol] = 0 # TODO check if valid for organics elif transfer_setup == 'metals': # renamed from radionuclides Bokna_137Cs self.num_lmm = self.specie_name2num('LMM') self.num_prev = self.specie_name2num('Particle reversible') self.num_srev = self.specie_name2num('Sediment reversible') self.num_psrev = self.specie_name2num('Particle slowly reversible') self.num_ssrev = self.specie_name2num('Sediment slowly reversible') # Values from Simonsen et al (2019a) Kd = self.get_config('chemical:transformations:Kd') # (m3/Kg) Dc = self.get_config('chemical:transformations:Dc') # (1/s) slow_coeff = self.get_config('chemical:transformations:slow_coeff') concSPM = 1.e-3 # concentration of available suspended particulate matter (kg/m3) sed_L = self.get_config('chemical:sediment:mixing_depth') # sediment mixing depth (m) sed_dens = self.get_config('chemical:sediment:density') # default particle density (kg/m3) sed_f = self.get_config('chemical:sediment:effective_fraction') # fraction of effective sorbents sed_phi = self.get_config('chemical:sediment:corr_factor') # sediment correction factor sed_poro = self.get_config('chemical:sediment:porosity') # sediment porosity sed_H = self.get_config('chemical:sediment:layer_thickness') # thickness of seabed interaction layer (m) #self.k_ads = Dc * Kd * 1e3 # L/(Kg*s) self.transfer_rates[self.num_lmm,self.num_prev] = Dc * Kd * concSPM self.transfer_rates[self.num_prev,self.num_lmm] = Dc self.transfer_rates[self.num_lmm,self.num_srev] = \ Dc * Kd * sed_L * sed_dens * (1.-sed_poro) * sed_f * sed_phi / sed_H self.transfer_rates[self.num_srev,self.num_lmm] = Dc * sed_phi self.transfer_rates[self.num_srev,self.num_ssrev] = slow_coeff self.transfer_rates[self.num_prev,self.num_psrev] = slow_coeff self.transfer_rates[self.num_ssrev,self.num_srev] = slow_coeff*.1 self.transfer_rates[self.num_psrev,self.num_prev] = slow_coeff*.1 elif transfer_setup == '137Cs_rev': self.num_lmm = self.specie_name2num('LMM') self.num_prev = self.specie_name2num('Particle reversible') self.num_srev = self.specie_name2num('Sediment reversible') # Simpler version of Values from Simonsen et al (2019a) # Only consider the reversible fraction Kd = self.get_config('chemical:transformations:Kd') Dc = self.get_config('chemical:transformations:Dc') concSPM = 1.e-3 # concentration of available suspended particulate matter (kg/m3) sed_L = self.get_config('chemical:sediment:mixing_depth') # sediment mixing depth (m) sed_dens = self.get_config('chemical:sediment:density') # default particle density (kg/m3) sed_f = self.get_config('chemical:sediment:effective_fraction') # fraction of effective sorbents sed_phi = self.get_config('chemical:sediment:corr_factor') # sediment correction factor sed_poro = self.get_config('chemical:sediment:porosity') # sediment porosity sed_H = self.get_config('chemical:sediment:layer_thickness') # thickness of seabed interaction layer (m) self.transfer_rates[self.num_lmm,self.num_prev] = Dc * Kd * concSPM self.transfer_rates[self.num_prev,self.num_lmm] = Dc self.transfer_rates[self.num_lmm,self.num_srev] = \ Dc * Kd * sed_L * sed_dens * (1.-sed_poro) * sed_f * sed_phi / sed_H self.transfer_rates[self.num_srev,self.num_lmm] = Dc * sed_phi elif transfer_setup=='custom': # Set of custom values for testing/development self.num_lmm = self.specie_name2num('LMM') if self.get_config('chemical:species:Colloid'): self.num_col = self.specie_name2num('Colloid') if self.get_config('chemical:species:Particle_reversible'): self.num_prev = self.specie_name2num('Particle reversible') if self.get_config('chemical:species:Sediment_reversible'): self.num_srev = self.specie_name2num('Sediment reversible') if self.get_config('chemical:slowly_fraction'): self.num_psrev = self.specie_name2num('Particle slowly reversible') self.num_ssrev = self.specie_name2num('Sediment slowly reversible') if self.get_config('chemical:irreversible_fraction'): self.num_pirrev = self.specie_name2num('Particle irreversible') self.num_sirrev = self.specie_name2num('Sediment irreversible') if self.get_config('chemical:species:Particle_reversible'): self.transfer_rates[self.num_lmm,self.num_prev] = 5.e-6 #*0. self.transfer_rates[self.num_prev,self.num_lmm] = \ self.get_config('chemical:transformations:Dc') if self.get_config('chemical:species:Sediment_reversible'): self.transfer_rates[self.num_lmm,self.num_srev] = 1.e-5 #*0. self.transfer_rates[self.num_srev,self.num_lmm] = \ self.get_config('chemical:transformations:Dc') * self.get_config('chemical:sediment:corr_factor') # self.transfer_rates[self.num_srev,self.num_lmm] = 5.e-6 if self.get_config('chemical:slowly_fraction'): self.transfer_rates[self.num_prev,self.num_psrev] = 2.e-6 self.transfer_rates[self.num_srev,self.num_ssrev] = 2.e-6 self.transfer_rates[self.num_psrev,self.num_prev] = 2.e-7 self.transfer_rates[self.num_ssrev,self.num_srev] = 2.e-7 elif transfer_setup=='Sandnesfj_Al': # Use values from Simonsen et al (2019b) self.num_lmmanion = self.specie_name2num('LMManion') self.num_lmmcation = self.specie_name2num('LMMcation') self.num_humcol = self.specie_name2num('Humic colloid') self.num_polymer = self.specie_name2num('Polymer') self.num_prev = self.specie_name2num('Particle reversible') self.num_srev = self.specie_name2num('Sediment reversible') Dc = self.get_config('chemical:transformations:Dc') self.salinity_intervals = [0,1,10,20] # Resize transfer rates array self.transfer_rates = np.zeros([len(self.salinity_intervals),self.transfer_rates.shape[0],self.transfer_rates.shape[1]]) # Salinity interval 0-1 psu self.transfer_rates[0,self.num_lmmcation, self.num_humcol] = 1.2e-5 self.transfer_rates[0,self.num_lmmcation, self.num_prev] = 4.e-6 self.transfer_rates[0,self.num_humcol, self.num_lmmcation] = .3*Dc self.transfer_rates[0,self.num_humcol, self.num_prev] = 2.e-6 self.transfer_rates[0,self.num_prev, self.num_lmmcation] = .3*Dc self.transfer_rates[0,self.num_srev, self.num_lmmcation] = .03*Dc # Salinity interval 1-10 psu self.transfer_rates[1,self.num_lmmcation, self.num_humcol] = 1.e-5 self.transfer_rates[1,self.num_lmmcation, self.num_prev] = 3.e-6 self.transfer_rates[1,self.num_lmmcation, self.num_polymer] = 1.2e-4 self.transfer_rates[1,self.num_humcol, self.num_lmmcation] = 7.*Dc self.transfer_rates[1,self.num_humcol, self.num_prev] = 4.e-6 self.transfer_rates[1,self.num_prev, self.num_lmmcation] = .5*Dc self.transfer_rates[1,self.num_srev, self.num_lmmcation] = .05*Dc self.transfer_rates[1,self.num_lmmanion, self.num_polymer] = 5.e-6 self.transfer_rates[1,self.num_polymer, self.num_lmmanion] = 12.*Dc self.transfer_rates[1,self.num_polymer, self.num_prev] = 2.4e-5 # Salinity interval 10-20 psu self.transfer_rates[2,self.num_lmmcation, self.num_humcol] = 8.e-6 self.transfer_rates[2,self.num_lmmcation, self.num_prev] = 2.e-6 self.transfer_rates[2,self.num_lmmcation, self.num_polymer] = 1.4e-4 self.transfer_rates[2,self.num_humcol, self.num_lmmcation] = 7.*Dc self.transfer_rates[2,self.num_humcol, self.num_prev] = 6.e-6 self.transfer_rates[2,self.num_prev, self.num_lmmcation] = .6*Dc self.transfer_rates[2,self.num_srev, self.num_lmmcation] = .06*Dc self.transfer_rates[2,self.num_lmmanion, self.num_polymer] = 5.e-6 self.transfer_rates[2,self.num_polymer, self.num_lmmanion] = 12.*Dc self.transfer_rates[2,self.num_polymer, self.num_prev] = 6.e-5 # Salinity interval >20 psu self.transfer_rates[3,self.num_lmmcation, self.num_humcol] = 6.e-6 self.transfer_rates[3,self.num_lmmcation, self.num_prev] = 1.8e-6 self.transfer_rates[3,self.num_lmmcation, self.num_polymer] = 1.5e-4 self.transfer_rates[3,self.num_humcol, self.num_lmmcation] = 7.*Dc self.transfer_rates[3,self.num_humcol, self.num_prev] = 1.e-5 self.transfer_rates[3,self.num_prev, self.num_lmmcation] = .8*Dc self.transfer_rates[3,self.num_srev, self.num_lmmcation] = .08*Dc self.transfer_rates[3,self.num_lmmanion, self.num_polymer] = 5.e-6 self.transfer_rates[3,self.num_polymer, self.num_lmmanion] = 12.*Dc self.transfer_rates[3,self.num_polymer, self.num_prev] = 8.e-5 else: logger.ERROR('No transfer setup available') # Set diagonal to 0. (not possible to transform to present specie) if len(self.transfer_rates.shape) == 3: for ii in range(self.transfer_rates.shape[0]): np.fill_diagonal(self.transfer_rates[ii,:,:],0.) else: np.fill_diagonal(self.transfer_rates,0.) # # HACK : # self.transfer_rates[:] = 0. # print ('\n ###### \n IMPORTANT:: \n transfer rates have been hacked! \n#### \n ') logger.debug('nspecies: %s' % self.nspecies) logger.debug('Transfer rates:\n %s' % self.transfer_rates)
[docs] def update_terminal_velocity(self, Tprofiles=None, Sprofiles=None, z_index=None): """Calculate terminal velocity for Pelagic Egg according to S. Sundby (1983): A one-dimensional model for the vertical distribution of pelagic fish eggs in the mixed layer Deep Sea Research (30) pp. 645-661 Method copied from ibm.f90 module of LADIM: Vikebo, F., S. Sundby, B. Aadlandsvik and O. Otteraa (2007), Fish. Oceanogr. (16) pp. 216-228 """ g = 9.81 # ms-2 # Particle properties that determine settling velocity partsize = self.elements.diameter # 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) DENSw = self.sea_water_density(T=T0, S=S0) DENSpart = self.elements.density dr = DENSw-DENSpart # density difference # 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 # terminal velocity for low Reynolds numbers W = (1.0/my_w)*(1.0/18.0)*g*partsize**2 * dr #W=np.zeros_like(W) #Setting to zero for debugging self.elements.terminal_velocity = W self.elements.terminal_velocity = W * self.elements.moving
[docs] def update_transfer_rates(self): '''Pick out the correct row from transfer_rates for each element. Modify the transfer rates according to local environmental conditions ''' transfer_setup=self.get_config('chemical:transfer_setup') if transfer_setup == 'metals' or \ transfer_setup=='custom' or \ transfer_setup=='137Cs_rev'or \ transfer_setup=='organics': self.elements.transfer_rates1D = self.transfer_rates[self.elements.specie,:] diss = self.get_config('chemical:transformations:dissociation') # Updating desorption rates according to local temperature, salinity, pH if transfer_setup=='organics' and diss=='nondiss': # filtering out zero values from temperature and salinity # TODO: Find out if problem is in the reader or in the data temperature=self.environment.sea_water_temperature temperature[temperature==0]=np.median(temperature) salinity=self.environment.sea_water_salinity salinity[salinity==0]=np.median(salinity) KOWTref = self.get_config('chemical:transformations:TrefKOW') DH_KOC_Sed = self.get_config('chemical:transformations:DeltaH_KOC_Sed') DH_KOC_DOM = self.get_config('chemical:transformations:DeltaH_KOC_DOM') Setchenow = self.get_config('chemical:transformations:Setchenow') tempcorrSed = self.tempcorr("Arrhenius",DH_KOC_Sed,temperature,KOWTref) tempcorrDOM = self.tempcorr("Arrhenius",DH_KOC_DOM,temperature,KOWTref) salinitycorr = self.salinitycorr(Setchenow,temperature,salinity) # Temperature and salinity correction for desorption rates (inversely proportional to Kd) self.elements.transfer_rates1D[self.elements.specie==self.num_humcol,self.num_lmm] = \ self.k21_0 / tempcorrDOM[self.elements.specie==self.num_humcol] / salinitycorr[self.elements.specie==self.num_humcol] self.elements.transfer_rates1D[self.elements.specie==self.num_prev,self.num_lmm] = \ self.k31_0 / tempcorrSed[self.elements.specie==self.num_prev] / salinitycorr[self.elements.specie==self.num_prev] self.elements.transfer_rates1D[self.elements.specie==self.num_srev,self.num_lmm] = \ self.k41_0 / tempcorrSed[self.elements.specie==self.num_srev] / salinitycorr[self.elements.specie==self.num_srev] elif transfer_setup=='organics' and diss!='nondiss': # Select elements for updating trasfer rates in sediments, SPM, and DOM #Sediments S = (self.elements.specie == self.num_srev) \ + (self.elements.specie == self.num_ssrev) SPM = (self.elements.specie == self.num_prev) DOM = (self.elements.specie == self.num_humcol) pH_sed = self.environment.pH_sediment[S] # pH_sed[pH_sed==0]=np.median(pH_sed) pH_water_SPM=self.environment.sea_water_ph_reported_on_total_scale[SPM] # pH_water_SPM[pH_water_SPM==0]=np.median(TW) pH_water_DOM=self.environment.sea_water_ph_reported_on_total_scale[DOM] # pH_water_DOM[pH_water_DOM==0]=np.median(pH_water_DOM) pKa_acid = self.get_config('chemical:transformations:pKa_acid') if pKa_acid < 0: raise ValueError("pKa_acid must be positive") # print("pKa_acid must be positive") # UserWarning(("pKa_acid must be positive")) else: pass pKa_base = self.get_config('chemical:transformations:pKa_base') if pKa_base < 0: raise ValueError("pKa_base must be positive") # print("pKa_base must be positive") # UserWarning(("pKa_base must be positive")) else: pass KOW = 10**self.get_config('chemical:transformations:LogKOW') KOC_sed_n = self.get_config('chemical:transformations:KOC_sed') if KOC_sed_n < 0: KOC_sed_n = 2.62 * KOW**0.82 # (L/KgOC), Park and Clough, 2014 (334)/Org2C TO DO Add if choice between input and estimation else: pass KOC_DOM_n = self.get_config('chemical:transformations:KOC_DOM') if KOC_DOM_n < 0: KOC_DOM_n = 2.88 * KOW**0.67 # (L/KgOC), Park and Clough, 2014 TO DO Add if choice between input and estimation else: pass fOC_SPM = self.get_config('chemical:transformations:fOC_SPM') # typical values from 0.01 to 0.1 gOC/g fOC_sed = self.get_config('chemical:transformations:fOC_sed') Org2C = 0.526 # kgOC/KgOM # Calculate original KOC_Values # TO DO: Store directly KOC values KOC_sed_initial = (self.Kd_sed)/fOC_sed # L/Kg / KgOC/Kg = L/KgOC KOC_SPM_initial = (self.Kd_SPM)/fOC_SPM # L/Kg / KgOC/Kg = L/KgOC KOC_DOM_initial = (self.Kd_DOM)/Org2C # filtering out zero values from temperature and salinity # TODO: Find out if problem is in the reader or in the data temperature=self.environment.sea_water_temperature temperature[temperature==0]=np.median(temperature) salinity=self.environment.sea_water_salinity salinity[salinity==0]=np.median(salinity) KOWTref = self.get_config('chemical:transformations:TrefKOW') DH_KOC_Sed = self.get_config('chemical:transformations:DeltaH_KOC_Sed') DH_KOC_DOM = self.get_config('chemical:transformations:DeltaH_KOC_DOM') Setchenow = self.get_config('chemical:transformations:Setchenow') tempcorrSed = self.tempcorr("Arrhenius",DH_KOC_Sed,temperature,KOWTref) tempcorrDOM = self.tempcorr("Arrhenius",DH_KOC_DOM,temperature,KOWTref) salinitycorr = self.salinitycorr(Setchenow,temperature,salinity) KOC_sedcorr = self.calc_KOC_sedcorr(KOC_sed_initial, KOC_sed_n, pKa_acid, pKa_base, KOW, pH_sed, diss) KOC_watcorrSPM = self.calc_KOC_watcorrSPM(KOC_SPM_initial, KOC_sed_n, pKa_acid, pKa_base, KOW, pH_water_SPM, diss) KOC_watcorrDOM = self.calc_KOC_watcorrDOM(KOC_DOM_initial, KOC_DOM_n, pKa_acid, pKa_base, KOW, pH_water_DOM, diss) # Temperature and salinity correction for desorption rates (inversely proportional to Kd) #### self.elements.transfer_rates1D[self.elements.specie==self.num_humcol,self.num_lmm] = \ self.k21_0 * KOC_watcorrDOM / tempcorrDOM[self.elements.specie==self.num_humcol] / salinitycorr[self.elements.specie==self.num_humcol] self.elements.transfer_rates1D[self.elements.specie==self.num_prev,self.num_lmm] = \ self.k31_0 * KOC_watcorrSPM / tempcorrSed[self.elements.specie==self.num_prev] / salinitycorr[self.elements.specie==self.num_prev] self.elements.transfer_rates1D[self.elements.specie==self.num_srev,self.num_lmm] = \ self.k41_0 * KOC_sedcorr / tempcorrSed[self.elements.specie==self.num_srev] / salinitycorr[self.elements.specie==self.num_srev] # Updating sorption rates if transfer_setup=='organics': # Updating sorption rates according to local SPM concentration concSPM=self.environment.spm * 1e-6 # (Kg/L) from (g/m3) # Apply SPM concentration profile if SPM reader has not depth coordinate # SPM concentration is kept constant to surface value in the mixed layer # Exponentially decreasing with depth below the mixed layers if not self.SPM_vertical_levels_given: lowerMLD = self.elements.z < -self.environment.ocean_mixed_layer_thickness #concSPM[lowerMLD] = concSPM[lowerMLD]/2 concSPM[lowerMLD] = concSPM[lowerMLD] * np.exp( -(self.elements.z[lowerMLD]+self.environment.ocean_mixed_layer_thickness[lowerMLD]) *np.log(0.5)/self.get_config('chemical:particle_concentration_half_depth') ) self.elements.transfer_rates1D[self.elements.specie==self.num_lmm,self.num_prev] = \ self.k_ads * concSPM[self.elements.specie==self.num_lmm] # k13 if transfer_setup == 'metals': # Updating sorption rates according to local SPM concentration and salinity concSPM=self.environment.spm * 1e-3 # (Kg/m3) from (g/m3) salinity=self.environment.sea_water_salinity # Apply SPM concentration profile if SPM reader has not depth coordinate # SPM concentration is kept constant to surface value in the mixed layer # Exponentially decreasing with depth below the mixed layers if not self.SPM_vertical_levels_given: lowerMLD = self.elements.z < -self.environment.ocean_mixed_layer_thickness #concSPM[lowerMLD] = concSPM[lowerMLD]/2 concSPM[lowerMLD] = concSPM[lowerMLD] * np.exp( -(self.elements.z[lowerMLD]+self.environment.ocean_mixed_layer_thickness[lowerMLD]) *np.log(0.5)/self.get_config('chemical:particle_concentration_half_depth') ) Kd0 = self.get_config('chemical:transformations:Kd') # (m3/Kg) S0 = self.get_config('chemical:transformations:S0') # (PSU) Dc = self.get_config('chemical:transformations:Dc') # (1/s) sed_L = self.get_config('chemical:sediment:mixing_depth') # sediment mixing depth (m) sed_dens = self.get_config('chemical:sediment:density') # default particle density (kg/m3) sed_f = self.get_config('chemical:sediment:effective_fraction') # fraction of effective sorbents sed_phi = self.get_config('chemical:sediment:corr_factor') # sediment correction factor sed_poro = self.get_config('chemical:sediment:porosity') # sediment porosity sed_H = self.get_config('chemical:sediment:layer_thickness') # thickness of seabed interaction layer (m) # Adjust Kd for salinity according to Perianez 2018 https://doi.org/10.1016/j.jenvrad.2018.02.014 if S0>0: Kd=Kd0*(S0+salinity[self.elements.specie==self.num_lmm])/S0 self.elements.transfer_rates1D[self.elements.specie==self.num_lmm,self.num_prev] = \ Dc * Kd * concSPM[self.elements.specie==self.num_lmm] # k13 self.elements.transfer_rates1D[self.elements.specie==self.num_lmm,self.num_srev] = \ Dc * Kd * sed_L * sed_dens * (1.-sed_poro) * sed_f * sed_phi / sed_H if transfer_setup=='organics': # Updating sorption rates according to local DOC concentration concDOM = self.environment.doc * 12e-6 / 1.025 / 0.526 * 1e-3 # (Kg[OM]/L) from (umol[C]/Kg) # Apply DOC concentration profile if DOC reader has not depth coordinate # DOC concentration is kept constant to surface value in the mixed layer # Exponentially decreasing with depth below the mixed layers if not self.DOC_vertical_levels_given: lowerMLD = self.elements.z < -self.environment.ocean_mixed_layer_thickness #concDOM[lowerMLD] = concDOM[lowerMLD]/2 concDOM[lowerMLD] = concDOM[lowerMLD] * np.exp( -(self.elements.z[lowerMLD]+self.environment.ocean_mixed_layer_thickness[lowerMLD]) *np.log(0.5)/self.get_config('chemical:doc_concentration_half_depth') ) self.elements.transfer_rates1D[self.elements.specie==self.num_lmm,self.num_humcol] = \ self.k_ads * concDOM[self.elements.specie==self.num_lmm] # k12 if self.get_config('chemical:species:Sediment_reversible'): # Only LMM chemicals close to seabed are allowed to interact with sediments # minimum height/maximum depth for each particle Zmin = -1.*self.environment.sea_floor_depth_below_sea_level interaction_thick = self.get_config('chemical:sediment:layer_thickness') # thickness of seabed interaction layer (m) dist_to_seabed = self.elements.z - Zmin self.elements.transfer_rates1D[(self.elements.specie == self.num_lmm) & (dist_to_seabed > interaction_thick), self.num_srev] = 0. elif transfer_setup=='Sandnesfj_Al': sal = self.environment.sea_water_salinity sali = np.searchsorted(self.salinity_intervals, sal) - 1 self.elements.transfer_rates1D = self.transfer_rates[sali,self.elements.specie,:]
[docs] def update_partitioning(self): '''Check if transformation processes shall occur Do transformation (change value of self.elements.specie) Update element properties for the transformed elements ''' specie_in = self.elements.specie.copy() # for storage of the initial partitioning specie_out = self.elements.specie.copy() # for storage of the final partitioning deltat = self.time_step.total_seconds() # length of a time step phaseshift = np.array(self.num_elements_active()*[False]) # Denotes which trajectory that shall be transformed p = 1. - np.exp(-self.elements.transfer_rates1D*deltat) # Probability for transformation psum = np.sum(p,axis=1) ran1=np.random.random(self.num_elements_active()) # Transformation where ran1 < total probability for transformation phaseshift[ ran1 < psum ] = True logger.info('Number of transformations: %s' % sum(phaseshift)) if sum(phaseshift) == 0: return ran4 = np.random.random(sum(phaseshift)) # New random number to decide which specie to end up in ttmp=[] # list for storing the out specie # Loop through each trajectory for ii in range(sum(phaseshift)): # Compare random number to the relative probability for each transfer process ttmp.append(np.searchsorted(np.cumsum(p[phaseshift][ii]/psum[phaseshift][ii]),ran4[ii])) specie_out[phaseshift] = np.array(ttmp) # Set the new partitioning self.elements.specie = specie_out logger.debug('old species: %s' % specie_in[phaseshift]) logger.debug('new species: %s' % specie_out[phaseshift]) for iin in range(self.nspecies): for iout in range(self.nspecies): self.ntransformations[iin,iout]+=sum((specie_in[phaseshift]==iin) & (specie_out[phaseshift]==iout)) logger.debug('Number of transformations total:\n %s' % self.ntransformations ) # Update Chemical properties after transformations self.update_chemical_diameter(specie_in, specie_out) self.sorption_to_sediments(specie_in, specie_out) self.desorption_from_sediments(specie_in, specie_out)
[docs] def sorption_to_sediments(self,sp_in=None,sp_out=None): '''Update Chemical properties when sorption to sediments occurs''' # Set z to local sea depth if self.get_config('chemical:species:LMM'): self.elements.z[(sp_out==self.num_srev) & (sp_in==self.num_lmm)] = \ -1.*self.environment.sea_floor_depth_below_sea_level[(sp_out==self.num_srev) & (sp_in==self.num_lmm)] self.elements.moving[(sp_out==self.num_srev) & (sp_in==self.num_lmm)] = 0 if self.get_config('chemical:species:LMMcation'): self.elements.z[(sp_out==self.num_srev) & (sp_in==self.num_lmmcation)] = \ -1.*self.environment.sea_floor_depth_below_sea_level[(sp_out==self.num_srev) & (sp_in==self.num_lmmcation)] self.elements.moving[(sp_out==self.num_srev) & (sp_in==self.num_lmmcation)] = 0 # avoid setting positive z values if np.nansum(self.elements.z>0): logger.debug('Number of elements lowered down to sea surface: %s' % np.nansum(self.elements.z>0)) self.elements.z[self.elements.z > 0] = 0
[docs] def desorption_from_sediments(self,sp_in=None,sp_out=None): '''Update Chemical properties when desorption from sediments occurs''' desorption_depth = self.get_config('chemical:sediment:desorption_depth') std = self.get_config('chemical:sediment:desorption_depth_uncert') if self.get_config('chemical:species:LMM'): self.elements.z[(sp_out==self.num_lmm) & (sp_in==self.num_srev)] = \ -1.*self.environment.sea_floor_depth_below_sea_level[(sp_out==self.num_lmm) & (sp_in==self.num_srev)] + desorption_depth self.elements.moving[(sp_out==self.num_lmm) & (sp_in==self.num_srev)] = 1 if std > 0: logger.debug('Adding uncertainty for desorption from sediments: %s m' % std) self.elements.z[(sp_out==self.num_lmm) & (sp_in==self.num_srev)] += np.random.normal( 0, std, sum((sp_out==self.num_lmm) & (sp_in==self.num_srev))) if self.get_config('chemical:species:LMMcation'): self.elements.z[(sp_out==self.num_lmmcation) & (sp_in==self.num_srev)] = \ -1.*self.environment.sea_floor_depth_below_sea_level[(sp_out==self.num_lmmcation) & (sp_in==self.num_srev)] + desorption_depth self.elements.moving[(sp_out==self.num_lmmcation) & (sp_in==self.num_srev)] = 1 if std > 0: logger.debug('Adding uncertainty for desorption from sediments: %s m' % std) self.elements.z[(sp_out==self.num_lmmcation) & (sp_in==self.num_srev)] += np.random.normal( 0, std, sum((sp_out==self.num_lmmcation) & (sp_in==self.num_srev))) # avoid setting positive z values if np.nansum(self.elements.z>0): logger.debug('Number of elements lowered down to sea surface: %s' % np.nansum(self.elements.z>0)) self.elements.z[self.elements.z > 0] = 0
[docs] def update_chemical_diameter(self,sp_in=None,sp_out=None): '''Update the diameter of the chemicals when specie is changed''' dia_part=self.get_config('chemical:particle_diameter') dia_diss=self.get_config('chemical:dissolved_diameter') # Transfer to reversible particles self.elements.diameter[(sp_out==self.num_prev) & (sp_in!=self.num_prev)] = dia_part # TODO Choose a proper diameter for aggregated particles if self.get_config('chemical:species:Humic_colloid'): self.elements.diameter[(sp_out==self.num_prev) & (sp_in==self.num_humcol)] = dia_part/2 logger.debug('Updated particle diameter for %s elements' % len(self.elements.diameter[(sp_out==self.num_prev) & (sp_in!=self.num_prev)])) std = self.get_config('chemical:particle_diameter_uncertainty') if std > 0: logger.debug('Adding uncertainty for particle diameter: %s m' % std) self.elements.diameter[(sp_out==self.num_prev) & (sp_in!=self.num_prev)] += np.random.normal( 0, std, sum((sp_out==self.num_prev) & (sp_in!=self.num_prev))) # Transfer to slowly reversible particles if self.get_config('chemical:slowly_fraction'): self.elements.diameter[(sp_out==self.num_psrev) & (sp_in!=self.num_psrev)] = dia_part if std > 0: logger.debug('Adding uncertainty for slowly rev particle diameter: %s m' % std) self.elements.diameter[(sp_out==self.num_psrev) & (sp_in!=self.num_psrev)] += np.random.normal( 0, std, sum((sp_out==self.num_psrev) & (sp_in!=self.num_psrev))) # Transfer to irreversible particles if self.get_config('chemical:irreversible_fraction'): self.elements.diameter[(sp_out==self.num_pirrev) & (sp_in!=self.num_pirrev)] = dia_part if std > 0: logger.debug('Adding uncertainty for irrev particle diameter: %s m' % std) self.elements.diameter[(sp_out==self.num_pirrev) & (sp_in!=self.num_pirrev)] += np.random.normal( 0, std, sum((sp_out==self.num_pirrev) & (sp_in!=self.num_pirrev))) # Transfer to LMM if self.get_config('chemical:species:LMM'): self.elements.diameter[(sp_out==self.num_lmm) & (sp_in!=self.num_lmm)] = dia_diss if self.get_config('chemical:species:LMManion'): self.elements.diameter[(sp_out==self.num_lmmanion) & (sp_in!=self.num_lmmanion)] = dia_diss if self.get_config('chemical:species:LMMcation'): self.elements.diameter[(sp_out==self.num_lmmcation) & (sp_in!=self.num_lmmcation)] = dia_diss # Transfer to colloids if self.get_config('chemical:species:Colloid'): self.elements.diameter[(sp_out==self.num_col) & (sp_in!=self.num_col)] = dia_diss if self.get_config('chemical:species:Humic_colloid'): self.elements.diameter[(sp_out==self.num_humcol) & (sp_in!=self.num_humcol)] = dia_diss if self.get_config('chemical:species:Polymer'): self.elements.diameter[(sp_out==self.num_polymer) & (sp_in!=self.num_polymer)] = dia_diss
[docs] def bottom_interaction(self,Zmin=None): ''' Change partitioning of chemicals that reach bottom due to settling. particle specie -> sediment specie ''' if not ((self.get_config('chemical:species:Particle_reversible')) & (self.get_config('chemical:species:Sediment_reversible')) or (self.get_config('chemical:slowly_fraction')) or (self.get_config('chemical:irreversible_fraction'))): return bottom = np.array(np.where(self.elements.z <= Zmin)[0]) kktmp = np.array(np.where(self.elements.specie[bottom] == self.num_prev)[0]) self.elements.specie[bottom[kktmp]] = self.num_srev self.ntransformations[self.num_prev,self.num_srev]+=len(kktmp) self.elements.moving[bottom[kktmp]] = 0 if self.get_config('chemical:slowly_fraction'): kktmp = np.array(np.where(self.elements.specie[bottom] == self.num_psrev)[0]) self.elements.specie[bottom[kktmp]] = self.num_ssrev self.ntransformations[self.num_psrev,self.num_ssrev]+=len(kktmp) self.elements.moving[bottom[kktmp]] = 0 if self.get_config('chemical:irreversible_fraction'): kktmp = np.array(np.where(self.elements.specie[bottom] == self.num_pirrev)[0]) self.elements.specie[bottom[kktmp]] = self.num_sirrev self.ntransformations[self.num_pirrev,self.num_sirrev]+=len(kktmp) self.elements.moving[bottom[kktmp]] = 0
[docs] def resuspension(self): """ Simple method to estimate the resuspension of sedimented particles, checking whether the current speed near the bottom is above a critical velocity Sediment species -> Particle specie """ # Exit function if particles and sediments not are present if not ((self.get_config('chemical:species:Particle_reversible')) & (self.get_config('chemical:species:Sediment_reversible'))): return specie_in = self.elements.specie.copy() critvel = self.get_config('chemical:sediment:resuspension_critvel') resusp_depth = self.get_config('chemical:sediment:resuspension_depth') std = self.get_config('chemical:sediment:resuspension_depth_uncert') Zmin = -1.*self.environment.sea_floor_depth_below_sea_level x_vel = self.environment.x_sea_water_velocity y_vel = self.environment.y_sea_water_velocity speed = np.sqrt(x_vel*x_vel + y_vel*y_vel) bottom = (self.elements.z <= Zmin) resusp = ( (bottom) & (speed >= critvel) ) if self.get_config('chemical:slowly_fraction'): resusp = ( resusp & (self.elements.specie!=self.num_ssrev) ) # Prevent ssrev (buried) to be resuspended # TODO buried sediment should be a new specie logger.info('Number of resuspended particles: {}'.format(np.sum(resusp))) self.elements.moving[resusp] = 1 self.elements.z[resusp] = Zmin[resusp] + resusp_depth if std > 0: logger.debug('Adding uncertainty for resuspension from sediments: %s m' % std) self.elements.z[resusp] += np.random.normal( 0, std, sum(resusp)) # avoid setting positive z values if np.nansum(self.elements.z>0): logger.debug('Number of elements lowered down to sea surface: %s' % np.nansum(self.elements.z>0)) self.elements.z[self.elements.z > 0] = 0 self.ntransformations[self.num_srev,self.num_prev]+=sum((resusp) & (self.elements.specie==self.num_srev)) self.elements.specie[(resusp) & (self.elements.specie==self.num_srev)] = self.num_prev if self.get_config('chemical:slowly_fraction'): self.ntransformations[self.num_ssrev,self.num_psrev]+=sum((resusp) & (self.elements.specie==self.num_ssrev)) self.elements.specie[(resusp) & (self.elements.specie==self.num_ssrev)] = self.num_psrev if self.get_config('chemical:irreversible_fraction'): self.ntransformations[self.num_sirrev,self.num_pirrev]+=sum((resusp) & (self.elements.specie==self.num_sirrev)) self.elements.specie[(resusp) & (self.elements.specie==self.num_sirrev)] = self.num_pirrev specie_out = self.elements.specie.copy() self.update_chemical_diameter(specie_in, specie_out)
[docs] def degradation(self): '''degradation.''' if self.get_config('chemical:transformations:degradation') is True: if self.get_config('chemical:transformations:degradation_mode')=='OverallRateConstants': # TODO: Rearrange code. Calculations here are for overall degradation including # degradation, photodegradation, and hydrolysys logger.debug('Calculating overall degradation using overall rate constants') degraded_now = np.zeros(self.num_elements_active()) # Degradation in the water k_W_tot = -np.log(0.5)/(self.get_config('chemical:transformations:t12_W_tot')*(60*60)) # (1/s) Tref_kWt = self.get_config('chemical:transformations:Tref_kWt') DH_kWt = self.get_config('chemical:transformations:DeltaH_kWt') W = (self.elements.specie == self.num_lmm) \ + (self.elements.specie == self.num_humcol) TW=self.environment.sea_water_temperature[W] TW[TW==0]=np.median(TW) k_W_fin = k_W_tot * self.tempcorr("Arrhenius",DH_kWt,TW,Tref_kWt) degraded_now[W] = self.elements.mass[W] * (1-np.exp(-k_W_fin * self.time_step.total_seconds())) # Degradation in the sediments k_S_tot = -np.log(0.5)/(self.get_config('chemical:transformations:t12_S_tot')*(60*60)) # (1/s) Tref_kSt = self.get_config('chemical:transformations:Tref_kSt') DH_kSt = self.get_config('chemical:transformations:DeltaH_kSt') S = (self.elements.specie == self.num_srev) \ + (self.elements.specie == self.num_ssrev) TS=self.environment.sea_water_temperature[S] TS[TS==0]=np.median(TS) k_S_fin = k_S_tot * self.tempcorr("Arrhenius",DH_kSt,TS,Tref_kSt) degraded_now[S] = self.elements.mass[S] * (1-np.exp(-k_S_fin * self.time_step.total_seconds())) self.elements.mass_degraded_water[W] = self.elements.mass_degraded_water[W] + degraded_now[W] self.elements.mass_degraded_sediment[S] = self.elements.mass_degraded_sediment[S] + degraded_now[S] self.elements.mass_degraded = self.elements.mass_degraded + degraded_now self.elements.mass = self.elements.mass - degraded_now self.deactivate_elements(self.elements.mass < (self.elements.mass + self.elements.mass_degraded + self.elements.mass_volatilized)/100, reason='removed') #to_deactivate = self.elements.mass < (self.elements.mass + self.elements.mass_degraded + self.elements.mass_volatilized)/100 #vol_morethan_degr = self.elements.mass_degraded >= self.elements.mass_volatilized # #self.deactivate_elements(to_deactivate + vol_morethan_degr, reason='volatilized') #self.deactivate_elements(to_deactivate + ~vol_morethan_degr, reason='degraded') else: pass
[docs] def volatilization(self): if self.get_config('chemical:transformations:volatilization') is True: logger.debug('Calculating: volatilization') volatilized_now = np.zeros(self.num_elements_active()) MolWtCO2=44 MolWtH2O=18 MolWt=self.get_config('chemical:transformations:MolWt') wind=5 # (m/s) (to read from atmosferic forcing) mixedlayerdepth=50 # m (to read from ocean forcing) Undiss_n=1 # 1 for PAHs Henry=self.get_config('chemical:transformations:Henry') # (atm m3/mol) Vp=self.get_config('chemical:transformations:Vpress') Tref_Vp=self.get_config('chemical:transformations:Tref_Vpress') DH_Vp=self.get_config('chemical:transformations:DeltaH_Vpress') Slb=self.get_config('chemical:transformations:Solub') Tref_Slb=self.get_config('chemical:transformations:Tref_Solub') DH_Slb=self.get_config('chemical:transformations:DeltaH_Solub') R=8.206e-05 #(atm m3)/(mol K) mixedlayerdepth = self.environment.ocean_mixed_layer_thickness # mask of dissolved elements within mixed layer W = (self.elements.specie == self.num_lmm) \ * (-self.elements.z <= mixedlayerdepth) # does volatilization apply only to num_lmm? # check mixedlayerdepth = mixedlayerdepth[W] T=self.environment.sea_water_temperature[W] T[T==0]=np.median(T) # temporary fix for missing values S=self.environment.sea_water_salinity[W] wind=(self.environment.x_wind[W]**2 + self.environment.y_wind[W]**2)**.5 Henry=( (Vp * self.tempcorr("Arrhenius",DH_Vp,T,Tref_Vp))) \ / (Slb * self.tempcorr("Arrhenius",DH_Slb,T,Tref_Slb)) \ * MolWt / 101325. # atm m3 mol-1 #k_S_fin = k_S_tot * self.tempcorr("Arrhenius",DH_kSt,TS,Tref_kSt) # Calculate mass transfer coefficient water side # Schwarzenbach et al., 2016 Eq.(19-20) MTCw = ((9e-4)+(7.2e-6*wind**3)) * (MolWtCO2/MolWt)**0.25 / Undiss_n # Calculate mass transfer coefficient air side # Schwarzenbach et al., 2016 Eq.(19-17)(19-18)(19-19) # Simple #MTCaH2O = 0.1 + 0.11 * wind # More complex Sca_H2O = 0.62 # 0.6 in the book. check MTCaH2O = 0.1 + wind*(6.1+0.63*wind)**0.5 \ /(13.3*(Sca_H2O)**0.5 + (6.1e-4+(6.3e-5)*wind)**-0.5 -5 + 1.25*np.log(Sca_H2O) ) MTCa = MTCaH2O * (MolWtH2O/MolWt)**(1/3) # Calculate overall volatilization mass tansfer coefficient HenryLaw = Henry * (1 + 0.01143 * S) / ( R * (T+273.15) ) MTCvol = 1 / ( 1/MTCw + 1/(MTCa * HenryLaw)) # (cm/s) #mixedlayerdepth = self.environment.ocean_mixed_layer_thickness[W] #Thick = np.clip(self.environment.sea_floor_depth_below_sea_level[W],0,mixedlayerdepth) # (m) Thick = mixedlayerdepth # Degubbing information to screen #print('################### Volatilization-info ##################') #print('Mixed Layer ',len(mixedlayerdepth),min(mixedlayerdepth),max(mixedlayerdepth),'m') #print('Temperature ',len(T),min(T),max(T),'C') #print('Salinity ',len(S),min(S),max(S)) #print('Henry ',len(Henry),min(Henry),max(Henry),'atm m3 / mol') #print('HenryLaw ',len(HenryLaw),min(HenryLaw),max(HenryLaw)) #print('wind ',len(wind),min(wind),max(wind), 'm/s') #print('MTCa ',len(MTCa),min(MTCa),max(MTCa),'cm/s') #print('MTCw ',len(MTCw),min(MTCw),max(MTCw),'cm/s') #print('MTCa*HenryLaw ',len(MTCa*HenryLaw),min(MTCa*HenryLaw),max(MTCa*HenryLaw),'cm/s') #print('MTCvol ',len(MTCvol),min(MTCvol),max(MTCvol),'cm/s') K_volatilization = 0.01 * MTCvol / Thick # (1/s) #logger.debug('MTCa: %s cm/s' % MTCa) #logger.debug('MTCw: %s cm/s' % MTCw) #logger.debug('Henry: %s ' % HenryLaw) #logger.debug('MTCvol: %s cm/s' % MTCvol) #logger.debug('T: %s C' % T) #logger.debug('S: %s ' % S) #logger.debug('Thick: %s ' % Thick) volatilized_now[W] = self.elements.mass[W] * (1-np.exp(-K_volatilization * self.time_step.total_seconds())) self.elements.mass_volatilized = self.elements.mass_volatilized + volatilized_now self.elements.mass = self.elements.mass - volatilized_now self.deactivate_elements(self.elements.mass < (self.elements.mass + self.elements.mass_degraded + self.elements.mass_volatilized)/100, reason='removed') #to_deactivate = self.elements.mass < (self.elements.mass + self.elements.mass_degraded + self.elements.mass_volatilized)/100 #vol_morethan_degr = self.elements.mass_degraded >= self.elements.mass_volatilized # #self.deactivate_elements(to_deactivate + vol_morethan_degr, reason='volatilized') #self.deactivate_elements(to_deactivate + ~vol_morethan_degr, reason='degraded') else: pass
[docs] def update(self): """Update positions and properties of Chemical particles.""" # Workaround due to conversion of datatype self.elements.specie = self.elements.specie.astype(np.int32) # Degradation and Volatilization if self.get_config('chemical:transfer_setup')=='organics': self.degradation() self.volatilization() # Dynamic Partitioning if self.get_config('chemical:dynamic_partitioning') is True: self.update_transfer_rates() self.update_partitioning() # Turbulent Mixing if self.get_config('drift:vertical_mixing') is True: self.update_terminal_velocity() self.vertical_mixing() else: self.update_terminal_velocity() self.vertical_buoyancy() # Resuspension self.resuspension() logger.info('partitioning: {} {}'.format([sum(self.elements.specie==ii) for ii in range(self.nspecies)],self.name_species)) # Horizontal advection self.advect_ocean_current() # Vertical advection if self.get_config('drift:vertical_advection') is True: self.vertical_advection() # Update transfer rates after last time step if self.time == (self.expected_end_time - self.time_step) or \ self.time == (self.expected_end_time) or \ self.num_elements_active() == 0 : self.update_transfer_rates()
# ################ # POSTPROCESSING
[docs] def simulation_summary(self, chemical_compound): '''Print a summary of the simulation: number of elements, number of transformations and final speciation ''' print(chemical_compound) print('Final speciation:') for isp,sp in enumerate(self.name_species): print ('{:32}: {:>6}'.format(sp,sum(self.elements.specie==isp))) print('Number of transformations:') for isp in range(self.nspecies): print('{}'.format(['{:>9}'.format(np.int32(item)) for item in self.ntransformations[isp,:]]) ) m_pre = sum(self.elements.mass)+sum(self.elements_deactivated.mass) m_deg = sum(self.elements.mass_degraded)+sum(self.elements_deactivated.mass_degraded) m_deg_w = sum(self.elements.mass_degraded_water)+sum(self.elements_deactivated.mass_degraded_water) m_deg_s = sum(self.elements.mass_degraded_sediment)+sum(self.elements_deactivated.mass_degraded_sediment) m_vol = sum(self.elements.mass_volatilized)+sum(self.elements_deactivated.mass_volatilized) m_tot = m_pre + m_deg + m_vol print('Mass balance:') print('mass preserved :', m_pre * 1e-6,' g ',m_pre/m_tot*100,'%') print('mass degraded :', m_deg * 1e-6,' g ',m_deg/m_tot*100,'%') print(' in water column :', m_deg_w * 1e-6,' g ',m_deg_w/m_tot*100,'%') print(' in sediments :', m_deg_s * 1e-6,' g ',m_deg_s/m_tot*100,'%') print('mass volatilized :', m_vol * 1e-6,' g ',m_vol/m_tot*100,'%')
[docs] def write_netcdf_chemical_density_map(self, filename, pixelsize_m='auto', zlevels=None, deltat=None, density_proj=None, llcrnrlon=None, llcrnrlat=None, urcrnrlon=None, urcrnrlat=None, mass_unit=None, time_avg_conc=False, horizontal_smoothing=False, smoothing_cells=0, reader_sea_depth=None, landmask_shapefile=None, origin_marker=None): '''Write netCDF file with map of Chemical species densities and concentrations''' from netCDF4 import Dataset, date2num #, stringtochar if landmask_shapefile is not None: if 'shape' in self.env.readers.keys(): # removing previously stored landmask del self.env.readers['shape'] # Adding new landmask from opendrift.readers import reader_shape custom_landmask = reader_shape.Reader.from_shpfiles(landmask_shapefile) self.add_reader(custom_landmask) elif 'global_landmask' not in self.env.readers.keys(): from opendrift.readers import reader_global_landmask global_landmask = reader_global_landmask.Reader() self.add_reader(global_landmask) if reader_sea_depth is not None: from opendrift.readers import reader_netCDF_CF_generic reader_sea_depth = reader_netCDF_CF_generic.Reader(reader_sea_depth) else: print('A reader for ''sea_floor_depth_below_sea_level'' must be specified') import sys sys.exit() # Temporary workaround if self.nspecies and self.name_species are not defined # TODO Make sure that these are saved when the simulation data is saved to the ncdf file # Then this workaround can be removed if not hasattr(self,'nspecies'): self.nspecies=4 if not hasattr(self,'name_species'): self.name_species = ['dissolved', 'DOC', 'SPM', 'sediment'] logger.info('Postprocessing: Write density and concentration to netcdf file') # Default bathymetry resolution 500x500. Can be increased (carefully) if high-res data is available and needed grid=np.meshgrid(np.linspace(llcrnrlon,urcrnrlon,500), np.linspace(llcrnrlat,urcrnrlat,500)) self.conc_lon=grid[0] self.conc_lat=grid[1] self.conc_topo=reader_sea_depth.get_variables_interpolated_xy(['sea_floor_depth_below_sea_level'], x = self.conc_lon.flatten(), y = self.conc_lat.flatten(), time = reader_sea_depth.times[0] if reader_sea_depth.times is not None else None )[0]['sea_floor_depth_below_sea_level'].reshape(self.conc_lon.shape) if pixelsize_m == 'auto': lon, lat = self.get_lonlats() latspan = lat.max()-lat.min() pixelsize_m=30 if latspan > .05: pixelsize_m = 50 if latspan > .1: pixelsize_m = 300 if latspan > .3: pixelsize_m = 500 if latspan > .7: pixelsize_m = 1000 if latspan > 2: pixelsize_m = 2000 if latspan > 5: pixelsize_m = 4000 if density_proj is None: # add default projection with equal-area property density_proj = pyproj.Proj('+proj=moll +ellps=WGS84 +lon_0=0.0') if mass_unit==None: mass_unit='microgram' # default unit for chemicals z = self.get_property('z')[0] if not zlevels==None: zlevels = np.sort(zlevels) z_array = np.append(np.append(-10000, zlevels) , max(0,np.nanmax(z))) else: z_array = [min(-10000,np.nanmin(z)), max(0,np.nanmax(z))] logger.info('vertical grid boundaries: {}'.format( [str(item) for item in z_array] ) ) # # H is array containing number of elements within each box defined by lon_array, lat_array and z_array H, lon_array, lat_array = \ self.get_chemical_density_array(pixelsize_m, z_array, density_proj=density_proj, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat, weight='mass',origin_marker=origin_marker) # calculating center point for eacxh pixel lon_array = (lon_array[:-1,:-1] + lon_array[1:,1:])/2 lat_array = (lat_array[:-1,:-1] + lat_array[1:,1:])/2 landmask = np.zeros_like(H[0,0,0,:,:]) if landmask_shapefile is not None: landmask = self.env.readers['shape'].__on_land__(lon_array,lat_array) else: landmask = self.env.readers['global_landmask'].__on_land__(lon_array,lat_array) if horizontal_smoothing: # Compute horizontally smoother field logger.debug('H.shape: ' + str(H.shape)) Hsm = np.zeros_like(H) for zi in range(len(z_array)-1): for sp in range(self.nspecies): for ti in range(H.shape[0]): Hsm[ti,sp,zi,:,:] = self.horizontal_smooth(H[ti,sp,zi,:,:],n=smoothing_cells) # Compute mean depth and volume in each pixel grid cell pixel_mean_depth = self.get_pixel_mean_depth(lon_array, lat_array) pixel_volume = np.zeros_like(H[0,0,:,:,:]) for zi,zz in enumerate(z_array[:-1]): topotmp = -pixel_mean_depth.copy() topotmp[np.where(topotmp < zz)] = zz topotmp = z_array[zi+1] - topotmp topotmp[np.where(topotmp < .1)] = 0. pixel_volume[zi,:,:] = topotmp * pixelsize_m**2 pixel_volume[np.where(pixel_volume==0.)] = np.nan # Compute mass of dry sediment in each pixel grid cell sed_L = self.get_config('chemical:sediment:mixing_depth') sed_dens = self.get_config('chemical:sediment:density') sed_poro = self.get_config('chemical:sediment:porosity') pixel_sed_mass = (pixelsize_m**2 *sed_L)*(1-sed_poro)*sed_dens # mass in kg # TODO this should be multiplied for the fraction og grid cell are that is not on land #conc = np.zeros_like(H) #if horizontal_smoothing: # conc_sm = np.zeros_like(Hsm) for ti in range(H.shape[0]): for sp in range(self.nspecies): if not self.name_species[sp].lower().startswith('sed'): #print('divide by volume') H[ti,sp,:,:,:] = H[ti,sp,:,:,:] / pixel_volume if horizontal_smoothing: Hsm[ti,sp,:,:,:] = Hsm[ti,sp,:,:,:] / pixel_volume elif self.name_species[sp].lower().startswith('sed'): #print('divide by mass') #print(pixel_sed_mass) H[ti,sp,:,:,:] = H[ti,sp,:,:,:] / pixel_sed_mass if horizontal_smoothing: Hsm[ti,sp,:,:,:] = Hsm[ti,sp,:,:,:] / pixel_sed_mass times = np.array( self.get_time_array()[0] ) if time_avg_conc: conctmp = H[:-1,:,:,:,:] cshape = conctmp.shape mdt = np.mean(times[1:] - times[:-1]) # output frequency in opendrift output file if deltat==None: ndt = 1 else: ndt = int( deltat / (mdt.total_seconds()/3600.) ) times2 = times[::ndt] times2 = times2[1:] odt = int(cshape[0]/ndt) logger.debug ('ndt '+ str(ndt)) # number of time steps over which to average in conc file logger.debug ('odt '+ str(odt)) # number of average slices Landmask=np.zeros_like(H[0:odt,:,:,:,:]) for zi in range(len(z_array)-1): for sp in range(self.nspecies): for ti in range(odt): Landmask[ti,sp,zi,:,:] = landmask # This may probably be written more efficiently! mean_conc = np.zeros( [odt,cshape[1],cshape[2],cshape[3],cshape[4]] ) for ii in range(odt): meantmp = np.mean(conctmp[(ii*ndt):(ii+1)*ndt,:,:,:,:],axis=0) mean_conc[ii,:,:,:,:] = meantmp nc = Dataset(filename, 'w') nc.createDimension('x', lon_array.shape[0]) nc.createDimension('y', lon_array.shape[1]) nc.createDimension('depth', len(z_array)-1) nc.createDimension('specie', self.nspecies) nc.createDimension('time', H.shape[0]) # times = self.get_time_array()[0] timestr = 'seconds since 1970-01-01 00:00:00' nc.createVariable('time', 'f8', ('time',)) nc.variables['time'][:] = date2num(times, timestr) nc.variables['time'].units = timestr nc.variables['time'].standard_name = 'time' if time_avg_conc: nc.createDimension('avg_time', odt) nc.createVariable('avg_time', 'f8', ('avg_time',)) nc.variables['avg_time'][:] = date2num(times2, timestr) # np.arange(mean_conc.shape[0]) nc.variables['avg_time'].units = timestr # Projection nc.createVariable('projection', 'i8') nc.variables['projection'].proj4 = density_proj.definition_string() # Cell size nc.createVariable('cell_size','f8') nc.variables['cell_size'][:] = pixelsize_m nc.variables['cell_size'].long_name = 'Length of cell' nc.variables['cell_size'].unit = 'm' nc.createVariable('smoothing_cells','i8') nc.variables['smoothing_cells'][:] = smoothing_cells nc.variables['smoothing_cells'].long_name = 'Number of cells in each direction for horizontal smoothing' nc.variables['smoothing_cells'].units = '1' # Coordinates nc.createVariable('lon', 'f8', ('y','x')) nc.createVariable('lat', 'f8', ('y','x')) nc.createVariable('depth', 'f8', ('depth',)) nc.createVariable('specie', 'i4', ('specie',)) nc.variables['lon'][:] = lon_array.T nc.variables['lon'].long_name = 'longitude' nc.variables['lon'].short_name = 'longitude' nc.variables['lon'].units = 'degrees_east' nc.variables['lat'][:] = lat_array.T nc.variables['lat'].long_name = 'latitude' nc.variables['lat'].short_name = 'latitude' nc.variables['lat'].units = 'degrees_north' nc.variables['depth'][:] = z_array[1:] nc.variables['specie'][:] = np.arange(self.nspecies) nc.variables['specie'].long_name = ' '.join(['{}:{}'.format(isp,sp) for isp,sp in enumerate(self.name_species)]) # Density #nc.createVariable('density', 'i4', # ('time','specie','depth','y', 'x'),fill_value=99999) #H = np.swapaxes(H, 3, 4).astype('i4') ##H = np.ma.masked_where(H==0, H) #nc.variables['density'][:] = H #nc.variables['density'].long_name = 'Number of elements in grid cell' #nc.variables['density'].grid_mapping = 'projection' #nc.variables['density'].units = '1' #if horizontal_smoothing: # nc.createVariable('density_smooth', 'f8', # ('time','specie','depth','y', 'x'),fill_value=1.e36) # Hsm = np.swapaxes(Hsm, 3, 4).astype('f8') # #Hsm = np.ma.masked_where(Hsm==0, Hsm) # nc.variables['density_smooth'][:] = Hsm # nc.variables['density_smooth'].long_name = 'Horizontally smoothed number of elements in grid cell' # nc.variables['density_smooth'].comment = 'Smoothed over '+str(smoothing_cells)+' grid points in all horizontal directions' # Chemical concentration if 0: nc.createVariable('concentration', 'f8', ('time','specie','depth','y', 'x'),fill_value=1.e36) H = np.ma.masked_where(Landmask==1,H) H = np.swapaxes(H, 3, 4) #.astype('i4') nc.variables['concentration'][:] = H nc.variables['concentration'].long_name = self.get_config('chemical:compound') +' concentration ' + '\n' + 'specie '+ \ ' '.join(['{}:{}'.format(isp,sp) for isp,sp in enumerate(self.name_species)]) nc.variables['concentration'].grid_mapping = 'projection_lonlat' nc.variables['concentration'].units = mass_unit+'/m3'+' (sed '+mass_unit+'/Kg)' # Chemical concentration, horizontally smoothed if horizontal_smoothing: nc.createVariable('concentration_smooth', 'f8', ('time','specie','depth','y', 'x'),fill_value=1.e36) Hsm = np.ma.masked_where(Landmask==1, Hsm) Hsm = np.swapaxes(Hsm, 3, 4) #.astype('i4') nc.variables['concentration_smooth'][:] = Hsm nc.variables['concentration_smooth'].long_name = self.get_config('chemical:compound') +' horizontally smoothed concentration ' + '\n' + 'specie '+ \ ' '.join(['{}:{}'.format(isp,sp) for isp,sp in enumerate(self.name_species)]) nc.variables['concentration_smooth'].grid_mapping = 'projection_lonlat' nc.variables['concentration_smooth'].units = mass_unit+'/m3'+' (sed '+mass_unit+'/Kg)' nc.variables['concentration_smooth'].comment = 'Smoothed over '+str(smoothing_cells)+' grid points in all horizontal directions' # Chemical concentration, time averaged if time_avg_conc: nc.createVariable('concentration_avg', 'f8', ('avg_time','specie','depth','y', 'x'),fill_value=+1.e36) mean_conc = np.ma.masked_where(Landmask[0:odt,:,:,:,:]==1, mean_conc) conc2 = np.swapaxes(mean_conc, 3, 4) #.astype('i4') #conc2 = np.ma.masked_where(landmask==1, conc2) nc.variables['concentration_avg'][:] = conc2 nc.variables['concentration_avg'].long_name = self.get_config('chemical:compound') + ' time averaged concentration ' + '\n' + 'specie '+ \ ' '.join(['{}:{}'.format(isp,sp) for isp,sp in enumerate(self.name_species)]) nc.variables['concentration_avg'].grid_mapping = 'projection_lonlat' nc.variables['concentration_avg'].units = mass_unit+'/m3'+' (sed '+mass_unit+'/Kg)' # Volume of boxes nc.createVariable('volume', 'f8', ('depth','y', 'x'),fill_value=0) pixel_volume = np.swapaxes(pixel_volume, 1, 2) #.astype('i4') pixel_volume = np.ma.masked_where(pixel_volume==0, pixel_volume) nc.variables['volume'][:] = pixel_volume nc.variables['volume'].long_name = 'Volume of grid cell (' + str(pixelsize_m)+'x'+str(pixelsize_m)+'m)' nc.variables['volume'].grid_mapping = 'projection_lonlat' nc.variables['volume'].units = 'm3' # Topography nc.createVariable('topo', 'f8', ('y', 'x'),fill_value=0) pixel_mean_depth = np.ma.masked_where(landmask==1, pixel_mean_depth) nc.variables['topo'][:] = pixel_mean_depth.T nc.variables['topo'].long_name = 'Depth of grid point' nc.variables['topo'].grid_mapping = 'projection_lonlat' nc.variables['topo'].units = 'm' # Binary mask nc.createVariable('land', 'i4', ('y', 'x'),fill_value=-1) #landmask = np.ma.masked_where(landmask==0, landmask) nc.variables['land'][:] = np.swapaxes(landmask,0,1).astype('i4') nc.variables['land'].long_name = 'Binary land mask' nc.variables['land'].grid_mapping = 'projection_lonlat' nc.variables['land'].units = 'm' nc.close() logger.info('Wrote to '+filename)
[docs] def get_chemical_density_array(self, pixelsize_m, z_array, density_proj=None, llcrnrlon=None,llcrnrlat=None, urcrnrlon=None,urcrnrlat=None, weight=None, origin_marker=None): ''' compute a particle concentration map from particle positions Use user defined projection (density_proj=<proj4_string>) or create a lon/lat grid (density_proj=None) ''' lon = self.get_property('lon')[0] lat = self.get_property('lat')[0] times = self.get_time_array()[0] # Redundant ?? if density_proj is None: # add default projection with equal-area property density_proj = pyproj.Proj('+proj=moll +ellps=WGS84 +lon_0=0.0') # create a grid in the specified projection x,y = density_proj(lon, lat) if llcrnrlon is not None: llcrnrx,llcrnry = density_proj(llcrnrlon,llcrnrlat) urcrnrx,urcrnry = density_proj(urcrnrlon,urcrnrlat) else: llcrnrx,llcrnry = x.min()-pixelsize_m, y.min()-pixelsize_m urcrnrx,urcrnry = x.max()+pixelsize_m, y.max()+pixelsize_m x_array = np.arange(llcrnrx,urcrnrx, pixelsize_m) y_array = np.arange(llcrnry,urcrnry, pixelsize_m) bins=(x_array, y_array) outsidex, outsidey = max(x_array)*1.5, max(y_array)*1.5 z = self.get_property('z')[0] if weight is not None: weight_array = self.get_property(weight)[0] status = self.get_property('status')[0] specie = self.get_property('specie')[0] if origin_marker is not None: originmarker = self.get_property('origin_marker')[0] Nspecies = self.nspecies H = np.zeros((len(times), Nspecies, len(z_array) - 1, len(x_array) - 1, len(y_array) - 1 )) for sp in range(Nspecies): for i in range(len(times)): if weight is not None: weights = weight_array[i,:] if origin_marker is not None: weight_array[i,:] = weight_array[i,:] * (originmarker[i,:]==origin_marker) else: weights = None for zi in range(len(z_array)-1): kktmp = ( (specie[i,:]==sp) & (z[i,:]>z_array[zi]) & (z[i,:]<=z_array[zi+1]) ) H[i,sp,zi,:,:], dummy, dummy = \ np.histogram2d(x[i,kktmp], y[i,kktmp], weights=weight_array[i,kktmp], bins=bins) if density_proj is not None: Y,X = np.meshgrid(y_array, x_array) lon_array, lat_array = density_proj(X,Y,inverse=True) return H, lon_array, lat_array
[docs] def get_pixel_mean_depth(self,lons,lats): from scipy import interpolate # Ocean model depth and lat/lon h_grd = self.conc_topo h_grd[np.isnan(h_grd)] = 0. nx = h_grd.shape[0] ny = h_grd.shape[1] lat_grd = self.conc_lat[:nx,:ny] lon_grd = self.conc_lon[:nx,:ny] # Interpolate topography to new grid h = interpolate.griddata((lon_grd.flatten(),lat_grd.flatten()), h_grd.flatten(), (lons, lats), method='linear') return h
[docs] def horizontal_smooth(self, a, n=0): if n==0: num_coarse=a return num_coarse xdm=a.shape[1] ydm=a.shape[0] #msk = self.conc_mask b=np.zeros([ydm+2*n,xdm+2*n],dtype=int) b[n:-n,n:-n]=a num_coarse = np.zeros([ydm,xdm],dtype=float) smo_tmp1=np.zeros([ydm,xdm]) #smo_msk1=np.zeros([ydm-2*n,xdm-2*n],dtype=float) nlayers = 0 for ism in np.arange(-n,n+1): for jsm in np.arange(-n,n+1): smo_tmp = b[n+jsm:ydm+n+jsm, n+ism:xdm+n+ism] smo_tmp1+=smo_tmp # Must preferrably take care of land points # smo_msk = msk[n+jsm:ydm-n+jsm, n+ism:xdm-n+ism] # smo_msk1+=smo_msk nlayers+=1 if n>0: # num_coarse[n:-n,n:-n] = smo_tmp1 / smo_msk1 num_coarse[:,:] = smo_tmp1 / nlayers else: num_coarse = smo_tmp1 # num_coarse = num_coarse*msk return num_coarse
[docs] def emission_factors(self, scrubber_type, chemical_compound): """Emission factors for heavy metals and PAHs in open loop and closed loop scrubbers Hermansson et al 2021 https://doi.org/10.1016/j.trd.2021.102912 bilge water, gray water, anti fouling paint, sewage, food waster from EMERGE Deliverable 2.1 ash (atmospheric depositions) from EMERGE Deliverable 3.2 """ emission_factors_open_loop = { # mean +/-95% # ug/L ug/L "Arsenic": [6.8, 3.4], "Cadmium": [0.8, 0.3], "Chromium": [15., 6.5], "Copper": [36., 12.], "Iron": [260., 250.], "Lead": [8.8, 4.4], "Mercury": [0.09, 0.01], "Nickel": [48., 12.], "Vanadium": [170., 49.], "Zinc": [110., 59.], "Cobalt": [0.17, 0.14], "Selenium": [97., 38], # "Naphthalene": [2.81, 0.77], "Phenanthrene": [1.51, 0.29], "Fluoranthene": [0.16, 0.04], "Benzo-a-anthracene": [0.12, 0.05], "Benzo-a-pyrene": [0.05, 0.02], "Dibenzo-ah-anthracene": [0.03, 0.01], # "Acenaphthylene": [0.12, 0.07], "Acenaphthene": [0.19, 0.07], "Fluorene": [0.46, 0.10], "Anthracene": [0.08, 0.04], "Pyrene": [0.31, 0.11], "Chrysene": [0.19, 0.07], "Benzo-b-fluoranthene": [0.04, 0.02], "Benzo-k-fluoranthene": [0.01, 0.01], "Indeno-123cd-pyrene": [0.07, 0.06], "Benzo-ghi-perylene": [0.02, 0.01], # "Nitrate": [2830., 2060.], "Nitrite": [760., 680.], "Ammonium": [730., 30.], "Sulphur": [2200000., 446000.], } emission_factors_closed_loop = { # mean +/-95% # ug/L ug/L "Arsenic": [22., 9.4], "Cadmium": [0.55, 0.19], "Chromium": [1300., 1700.], "Copper": [480., 230.], "Iron": [490., 82.], "Lead": [7.7, 3.1], "Mercury": [0.07, 0.02], "Nickel": [2700., 860.], "Vanadium": [9100., 3200.], "Zinc": [370., 200.], "Cobalt": [0., 0.], "Selenium": [0., 0.], # "Naphthalene": [2.08, 1.05], "Phenanthrene": [5.00, 2.30], "Fluoranthene": [0.63, 0.41], "Benzo-a-anthracene": [0.30, 0.29], "Benzo-a-pyrene": [0.06, 0.05], "Dibenzo-ah-anthracene": [0.03, 0.02], # "Acenaphthylene": [0.09, 0.06], "Acenaphthene": [0.47, 0.31], "Fluorene": [1.32, 0.54], "Anthracene": [1.55, 2.00], "Pyrene": [0.76, 0.59], "Chrysene": [0.50, 0.45], "Benzo-b-fluoranthene": [0.14, 0.12], "Benzo-k-fluoranthene": [0.02, 0.02], "Indeno-123-cd-pyrene": [0.04, 0.03], "Benzo-ghi-perylene": [0.07, 0.07], # "Nitrate": [110980., 100000.], "Nitrite": [55760., 55000.], "Ammonium": [0., 0.], "Sulphur": [12280000., 10104000.], } emission_factors_grey_water = { # mean +/-95% # ug/L ug/L "Arsenic": [5.98, 3.17], "Cadmium": [0.16, 0.09], "Chromium": [7.28, 2.06], "Copper": [267., 97.], "Lead": [25.6, 21.01], "Mercury": [0.16, 0.09], "Nickel": [25.0, 19.36], "Selenium": [16.1, 10.64], "Zinc": [517., 112.], # "Nitrogen": [28900., 0.0], } emission_factors_bilge_water = { # mean +/-95% # ug/L ug/L "Arsenic": [35.9, 33.2], "Cadmium": [0.32, 0.07], "Chromium": [16.3, 15.4], "Copper": [49.7, 22.9], "Lead": [3.0, 1.24], "Nickel": [71.1, 11.8], "Selenium": [2.95, 1.01], "Vanadium": [76.5, 22.4], "Zinc": [949., 660.], # "Nitrate": [110980., 100000.], "Nitrite": [55760., 55000.], "Ammonium": [0., 0.], "Sulphur": [12280000., 10104000.], # "Naphthalene": [50.6, 34.3], "Phenanthrene": [3.67, 2.51], "Fluoranthene": [0.60, 0.96], "Benzo(a)anthracene": [0.10, 0.18], "Benzo(a)pyrene": [0.10, 0.15], "Dibenzo(a,h)anthracene": [0.02, 0.01], # "Acenaphthylene": [0.29, 0.17], "Acenaphthene": [1.42, 0.86], "Fluorene": [3.33, 2.43], "Anthracene": [0.22, 0.14], "Pyrene": [1.23, 1.33], "Chrysene": [0.17, 0.25], "Benzo(b)fluoranthene": [0.09, 0.13], "Benzo(k)fluoranthene": [0.03, 0.00], "Indeno(1,2,3-cd)pyrene": [0.05, 0.06], "Benzo(ghi)perylene": [0.13, 0.16], } emission_factors_sewage_water = { # mean +/-95% # ug/L ug/L "Arsenic": [22.9, 7.4], "Cadmium": [0.12, 0.10], "Chromium": [11.9, 8.2], "Copper": [319, 190], "Lead": [6.5, 3.1], "Mercury": [0.22, 0.12], "Nickel": [32.3, 21.3], "Selenium": [43.7, 18.3], "Zinc": [395., 174.], # "Nitrogen": [430., 0.], } emission_factors_AFP = { # Copper = 63.546 g/mol # Zinc = 65.38 g/mol # CuPyr = 315.86 g/mol = Copper(II) pyrithione = 0.2112 of Cu # CuO = 79.55 g/mol = Copper(II) oxide = 0.7989 of Cu # Zineb = 275.7 g/mol = Zinc ethylenebis(dithiocarbamate) = 0.2371 of Zn # ZnO = 81.38 g/mol = Zinc(II) oxide = 0.8033 of Zn # ZPyr = 317.70 g/mol = Zinc(II) pyrithione = 0.2058 of Zn # mean +/-95% # ug/L ug/L "CuO_AFP": [0.7989, 0.], "CuPyr_AFP": [0.2112, 0.], "Zineb_AFP": [0.2371, 0.], "ZnO_AFP": [0.8033, 0.], "ZnPyr_AFP": [0.2058, 0.], } emission_factors_SILAM_ash = { # g/g "Aresenic": [8.09E-5], "Cadmium": [6.30E-6], "Chromium": [2.10E-4], "Copper": [2.52E-4], "Iron": [2.52E-2], "Mercury": [6.30E-6], "Nickel": [4.10E-2], "Lead": [1.16E-4], "Vanadium": [8.30E-2], "Zinc": [2.42E-3], } if scrubber_type=="open_loop": Emission_factors = emission_factors_open_loop.get(chemical_compound)[0] elif scrubber_type=="closed_loop": Emission_factors = emission_factors_closed_loop.get(chemical_compound)[0] elif scrubber_type=="bilge_water": Emission_factors = emission_factors_bilge_water.get(chemical_compound)[0] elif scrubber_type=="grey_water": Emission_factors = emission_factors_grey_water.get(chemical_compound)[0] elif scrubber_type=="sewage_water": Emission_factors = emission_factors_sewage_water.get(chemical_compound)[0] elif scrubber_type=="AFP": # Copper and Zinc from antifouling paint Emission_factors = 1e6*emission_factors_AFP.get(chemical_compound)[0] # 1g = 1e6 ug: AFP is expressed as g elif scrubber_type=="AFP_metals_total": Emission_factors = 1e6 # g to ug elif scrubber_type=="N_sewage": # Nitrogen from sewage Emission_factors = 1e9 # 1kg = 1e9 ug: N_sewage is expressed as kg elif scrubber_type=="N_foodwaste": # Nitrogen from foodwaste Emission_factors = 1e9 # 1kg = 1e9 ug: N_sewage is expressed as kg elif scrubber_type=="SILAM_metals": Emission_factors = 1e9 #+ 1kg = 1e9 ug: Lead and Cadmium depositions given in kg elif scrubber_type=="SILAM_metals_from_ash": Emission_factors = 1e9*emission_factors_SILAM_ash.get(chemical_compound)[0] # 1kg=1e9ug: Ash depositions given in kg return Emission_factors
# TODO: Add emission uncertainty based on 95% confidence interval
[docs] def seed_from_DataArray(self, steam, lowerbound=0, higherbound=np.inf, radius=0, scrubber_type="open_loop", chemical_compound="Copper", mass_element_ug=100e3, number_of_elements=None, **kwargs): """Seed elements based on a dataarray with STEAM emission data Arguments: steam: dataarray with steam emission data, with coordinates * latitude (latitude) float32 * longitude (longitude) float32 * time (time) datetime64[ns] radius: scalar, unit: meters lowerbound: scalar, elements with lower values are discarded """ if chemical_compound is None: chemical_compound = self.get_config('chemical:compound') #mass_element_ug=1e3 # 1e3 - 1 element is 1mg chemical #mass_element_ug=20e3 # 100e3 - 1 element is 100mg chemical #mass_element_ug=100e3 # 100e3 - 1 element is 100mg chemical #mass_element_ug=1e6 # 1e6 - 1 element is 1g chemical sel=np.where((steam > lowerbound) & (steam < higherbound)) t=steam.time[sel[0]].data la=steam.latitude[sel[1]].data lo=steam.longitude[sel[2]].data data=np.array(steam.data) if number_of_elements is not None: total_volume = np.sum(data[sel]) total_mass = total_volume * self.emission_factors(scrubber_type, chemical_compound) mass_element_ug = total_mass / number_of_elements mass_element_ug_0 = total_mass / number_of_elements for i in range(0,t.size): scrubberwater_vol_l=data[sel][i] mass_ug=scrubberwater_vol_l * self.emission_factors(scrubber_type, chemical_compound) if number_of_elements is None: number=np.array(mass_ug / mass_element_ug).astype('int') else: number=np.ceil(np.array(mass_ug / mass_element_ug_0)).astype('int') mass_element_ug=mass_ug/number time = datetime.utcfromtimestamp((t[i] - np.datetime64('1970-01-01T00:00:00')) / np.timedelta64(1, 's')) if number>0: z = -1*np.random.uniform(0, 1, number) self.seed_elements(lon=lo[i]*np.ones(number), lat=la[i]*np.ones(number), radius=radius, number=number, time=time, mass=mass_element_ug,mass_degraded=0,mass_volatilized=0, z=z, origin_marker=1) mass_residual = mass_ug - number*mass_element_ug if mass_residual>0 and number_of_elements is None: z = -1*np.random.uniform(0, 1, 1) self.seed_elements(lon=lo[i], lat=la[i], radius=radius, number=1, time=time, mass=mass_residual,mass_degraded=0,mass_volatilized=0, z=z, origin_marker=1)
seed_from_STEAM = seed_from_DataArray ''' Alias of seed_from_DataArray method for backward compatibility '''
[docs] def init_chemical_compound(self, chemical_compound = None): ''' Chemical parameters for a selection of PAHs: Naphthalene, Phenanthrene, Fluorene, Benzo-a-anthracene, Benzo-a-pyrene, Dibenzo-ah-anthracene Data collected from literature by Isabel Hanstein (University of Heidelberg / Norwegian Meteorological Insitute) Mattia Boscherini, Loris Calgaro (University Ca' Foscari, Venezia) Manuel Aghito (Norwegian Meteorological Institute / University of Bergen) ''' if chemical_compound is not None: self.set_config('chemical:compound',chemical_compound) if self.get_config('chemical:compound') is None: raise ValueError("Chemical compound not defined") if self.get_config('chemical:compound') in ["Naphthalene",'C1-Naphthalene','Acenaphthene','Acenaphthylene','Fluorene']: #partitioning self.set_config('chemical:transfer_setup','organics') self.set_config('chemical:transformations:dissociation','nondiss') self.set_config('chemical:transformations:LogKOW',3.361) self.set_config('chemical:transformations:TrefKOW',25) self.set_config('chemical:transformations:DeltaH_KOC_Sed',-21036) self.set_config('chemical:transformations:DeltaH_KOC_DOM',-25900) ### Phenanthrene value self.set_config('chemical:transformations:Setchenow', 0.2503) #degradation self.set_config('chemical:transformations:t12_W_tot', 224.08) self.set_config('chemical:transformations:Tref_kWt', 25) self.set_config('chemical:transformations:DeltaH_kWt', 50000) ### generic self.set_config('chemical:transformations:t12_S_tot', 5012.4) self.set_config('chemical:transformations:Tref_kSt', 25) self.set_config('chemical:transformations:DeltaH_kSt', 50000) ### generic #volatilization self.set_config('chemical:transformations:MolWt', 128.1705) self.set_config('chemical:transformations:Henry', 4.551e-4) self.set_config('chemical:transformations:Vpress', 11.2) self.set_config('chemical:transformations:Tref_Vpress', 25) self.set_config('chemical:transformations:DeltaH_Vpress', 55925) self.set_config('chemical:transformations:Solub', 31.4) self.set_config('chemical:transformations:Tref_Solub', 25) self.set_config('chemical:transformations:DeltaH_Solub', 25300) elif self.get_config('chemical:compound') in ["Phenanthrene",'Dibenzothiophene','C2-Naphthalene','Anthracene','C3-Naphthalene','C1-Dibenzothiophene']: #partitioning self.set_config('chemical:transfer_setup','organics') self.set_config('chemical:transformations:dissociation','nondiss') self.set_config('chemical:transformations:LogKOW',4.505) self.set_config('chemical:transformations:TrefKOW',25) self.set_config('chemical:transformations:DeltaH_KOC_Sed',-24900) self.set_config('chemical:transformations:DeltaH_KOC_DOM',-25900) self.set_config('chemical:transformations:Setchenow', 0.3026) #degradation self.set_config('chemical:transformations:t12_W_tot', 1125.79) self.set_config('chemical:transformations:Tref_kWt', 25) self.set_config('chemical:transformations:DeltaH_kWt', 50000) ### generic self.set_config('chemical:transformations:t12_S_tot', 29124.96) self.set_config('chemical:transformations:Tref_kSt', 25) self.set_config('chemical:transformations:DeltaH_kSt', 50000) ### generic #volatilization self.set_config('chemical:transformations:MolWt', 178.226) self.set_config('chemical:transformations:Henry', 4.294e-5) self.set_config('chemical:transformations:Vpress', 0.0222) self.set_config('chemical:transformations:Tref_Vpress', 25) self.set_config('chemical:transformations:DeltaH_Vpress', 71733) self.set_config('chemical:transformations:Solub', 1.09) self.set_config('chemical:transformations:Tref_Solub', 25) self.set_config('chemical:transformations:DeltaH_Solub', 34800) elif self.get_config('chemical:compound') in ["Fluoranthene",'Pyrene','C1-Phenanthrene','C2-Dibenzothiophene']: #partitioning self.set_config('chemical:transfer_setup','organics') self.set_config('chemical:transformations:dissociation','nondiss') self.set_config('chemical:transformations:LogKOW',5.089) self.set_config('chemical:transformations:TrefKOW',25) self.set_config('chemical:transformations:DeltaH_KOC_Sed',-47413) self.set_config('chemical:transformations:DeltaH_KOC_DOM',-27900) self.set_config('chemical:transformations:Setchenow', 0.2885) #degradation self.set_config('chemical:transformations:t12_W_tot', 1705.02) self.set_config('chemical:transformations:Tref_kWt', 25) self.set_config('chemical:transformations:DeltaH_kWt', 50000) ### generic self.set_config('chemical:transformations:t12_S_tot', 55000) self.set_config('chemical:transformations:Tref_kSt', 25) self.set_config('chemical:transformations:DeltaH_kSt', 50000) ### generic #volatilization self.set_config('chemical:transformations:MolWt', 202.26) self.set_config('chemical:transformations:Henry', 1.439e-5) self.set_config('chemical:transformations:Vpress', 0.00167) self.set_config('chemical:transformations:Tref_Vpress', 25) self.set_config('chemical:transformations:DeltaH_Vpress', 79581) self.set_config('chemical:transformations:Solub', 0.231) self.set_config('chemical:transformations:Tref_Solub', 25) self.set_config('chemical:transformations:DeltaH_Solub', 30315) elif self.get_config('chemical:compound') in ["Benzo-a-anthracene",'C2-Phenanthrene','Benzo-b-fluoranthene','Chrysene']: #partitioning self.set_config('chemical:transfer_setup','organics') self.set_config('chemical:transformations:dissociation','nondiss') self.set_config('chemical:transformations:LogKOW',5.724) self.set_config('chemical:transformations:TrefKOW',25) self.set_config('chemical:transformations:DeltaH_KOC_Sed', -38000) ### Pyrene value self.set_config('chemical:transformations:DeltaH_KOC_DOM', -25400) ### Pyrene value self.set_config('chemical:transformations:Setchenow', 0.3605) #degradation self.set_config('chemical:transformations:t12_W_tot', 1467.62) self.set_config('chemical:transformations:Tref_kWt', 25) self.set_config('chemical:transformations:DeltaH_kWt', 50000) ### generic self.set_config('chemical:transformations:t12_S_tot', 46600) self.set_config('chemical:transformations:Tref_kSt', 25) self.set_config('chemical:transformations:DeltaH_kSt', 50000) ### generic #volatilization self.set_config('chemical:transformations:MolWt', 228.29) self.set_config('chemical:transformations:Henry', 6.149e-6) self.set_config('chemical:transformations:Vpress', 0.0000204) self.set_config('chemical:transformations:Tref_Vpress', 25) self.set_config('chemical:transformations:DeltaH_Vpress', 100680) self.set_config('chemical:transformations:Solub', 0.011) self.set_config('chemical:transformations:Tref_Solub', 25) self.set_config('chemical:transformations:DeltaH_Solub', 46200) elif self.get_config('chemical:compound') in ["Benzo-a-pyrene",'C3-Dibenzothiophene','C3-Phenanthrene']: #partitioning self.set_config('chemical:transfer_setup','organics') self.set_config('chemical:transformations:dissociation','nondiss') self.set_config('chemical:transformations:LogKOW', 6.124) self.set_config('chemical:transformations:TrefKOW',25) self.set_config('chemical:transformations:DeltaH_KOC_Sed', -43700) ### mean value 16 PAHs self.set_config('chemical:transformations:DeltaH_KOC_DOM', -31280) self.set_config('chemical:transformations:Setchenow', 0.171) #degradation self.set_config('chemical:transformations:t12_W_tot', 1491.42) self.set_config('chemical:transformations:Tref_kWt', 25) self.set_config('chemical:transformations:DeltaH_kWt', 50000) ### generic self.set_config('chemical:transformations:t12_S_tot', 44934.76) self.set_config('chemical:transformations:Tref_kSt', 25) self.set_config('chemical:transformations:DeltaH_kSt', 50000) ### generic #volatilization self.set_config('chemical:transformations:MolWt', 252.32) self.set_config('chemical:transformations:Henry', 6.634e-7) self.set_config('chemical:transformations:Vpress', 0.00000136) self.set_config('chemical:transformations:Tref_Vpress', 25) self.set_config('chemical:transformations:DeltaH_Vpress', 107887) self.set_config('chemical:transformations:Solub', 0.00229) self.set_config('chemical:transformations:Tref_Solub', 25) self.set_config('chemical:transformations:DeltaH_Solub', 38000) elif self.get_config('chemical:compound') in ["Dibenzo-ah-anthracene",'Benzo-k-fluoranthene','Benzo-ghi-perylene','Indeno-123cd-pyrene']: #partitioning self.set_config('chemical:transfer_setup','organics') self.set_config('chemical:transformations:dissociation','nondiss') self.set_config('chemical:transformations:LogKOW', 6.618) self.set_config('chemical:transformations:TrefKOW',25) self.set_config('chemical:transformations:DeltaH_KOC_Sed', -43700) ### mean value 16 PAHs self.set_config('chemical:transformations:DeltaH_KOC_DOM', -30900) self.set_config('chemical:transformations:Setchenow', 0.338) #degradation self.set_config('chemical:transformations:t12_W_tot', 1464.67) self.set_config('chemical:transformations:Tref_kWt', 25) self.set_config('chemical:transformations:DeltaH_kWt', 50000) ### generic self.set_config('chemical:transformations:t12_S_tot', 40890.08) self.set_config('chemical:transformations:Tref_kSt', 25) self.set_config('chemical:transformations:DeltaH_kSt', 50000) ### generic #volatilization self.set_config('chemical:transformations:MolWt', 278.35) self.set_config('chemical:transformations:Henry', 4.894e-8) self.set_config('chemical:transformations:Vpress', 0.0000000427) self.set_config('chemical:transformations:Tref_Vpress', 25) self.set_config('chemical:transformations:DeltaH_Vpress', 112220) self.set_config('chemical:transformations:Solub', 0.00142) self.set_config('chemical:transformations:Tref_Solub', 25) self.set_config('chemical:transformations:DeltaH_Solub', 38000) ### Benzo-a-pyrene value elif self.get_config('chemical:compound') == "Copper": self.set_config('chemical:transfer_setup','metals') self.set_config('chemical:transformations:Kd', 60.1) # Tomczak et Al 2019 #self.set_config('chemical:transformations:Kd', 50) # Merlin Expo, high confidence self.set_config('chemical:transformations:S0', 17.0) # note below elif self.get_config('chemical:compound') == "Zinc": self.set_config('chemical:transfer_setup','metals') self.set_config('chemical:transformations:Kd', 173) # Tomczak et Al 2019 #self.set_config('chemical:transformations:Kd', 100) # Merlin Expo, high confidence self.set_config('chemical:transformations:S0', 17.0) # note below elif self.get_config('chemical:compound') == "Lead": self.set_config('chemical:transfer_setup','metals') self.set_config('chemical:transformations:Kd', 369) # Tomczak et Al 2019 #self.set_config('chemical:transformations:Kd', 500) # Merlin Expo, strong confidence self.set_config('chemical:transformations:S0', 17.0) # note below elif self.get_config('chemical:compound') == "Vanadium": self.set_config('chemical:transfer_setup','metals') self.set_config('chemical:transformations:Kd', 42.9) # Tomczak et Al 2019 #self.set_config('chemical:transformations:Kd', 5) # Merlin Expo, weak confidence self.set_config('chemical:transformations:S0', 17.0) # note below elif self.get_config('chemical:compound') == "Cadmium": self.set_config('chemical:transfer_setup','metals') self.set_config('chemical:transformations:Kd', 134) # Tomczak et Al 2019 #self.set_config('chemical:transformations:Kd', 79) # Merlin Expo, strong confidence #self.set_config('chemical:transformations:Kd', 6.6) # Turner Millward 2002 self.set_config('chemical:transformations:S0', 17.0) # note below elif self.get_config('chemical:compound') == "Chromium": self.set_config('chemical:transfer_setup','metals') self.set_config('chemical:transformations:Kd', 124) # Tomczak et Al 2019 #self.set_config('chemical:transformations:Kd', 130) # Cr(III) Merlin Expo, moderate confidence #self.set_config('chemical:transformations:Kd', 180) # Turner Millward 2002 self.set_config('chemical:transformations:S0', 17.0) # note below elif self.get_config('chemical:compound') == "Nickel": self.set_config('chemical:transfer_setup','metals') self.set_config('chemical:transformations:Kd', 31.1) # Tomczak et Al 2019 #self.set_config('chemical:transformations:Kd', 25) # Merlin Expo, strong confidence #self.set_config('chemical:transformations:Kd', 5.3) # Turner Millward 2002 self.set_config('chemical:transformations:S0', 17.0) # note below
# Default value for S0 is set to 17.0. This correspond to a Kd at salinity 35 being 32.7% # of the fresh water value, which was the average reduction obtained comparing the values # in Tomczak et Al 2019 to the "ocean margins" recommended values in IAEA TRS no.422, for a # selection of metals (Cd, Cr, Hg, Ni, Pb, Zn). This gives very similar results the value # 15.8, suggested in Perianez 2018. # https://doi.org/10.1016/j.apgeochem.2019.04.003 # https://www-pub.iaea.org/MTCD/Publications/PDF/TRS422_web.pdf # https://doi.org/10.1016/j.jenvrad.2018.02.014 # # Merlin Expo Kd values are mean values from Allison and Allison 2005 # https://cfpub.epa.gov/si/si_public_record_report.cfm?dirEntryId=135783
[docs] def plot_mass(self, legend=['dissolved','SPM','sediment'], mass_unit='g', time_unit='hours', title=None, filename=None, start_date=None): """Plot chemical mass distribution between the different species legend list of specie labels, for example ['dissolved','SPM','sediment'] mass_unit 'g','mg','ug' time_unit 'seconds', 'minutes', 'hours' , 'days' title figure title string """ import matplotlib.pyplot as plt fig, ax = plt.subplots() if not title == []: fig.suptitle(title) mass=self.get_property('mass') sp=self.get_property('specie') steps=self.steps_output bars=np.zeros((steps,5)) mass_conversion_factor=1e-6 if mass_unit=='g' and self.elements.variables['mass']['units']=='ug': mass_conversion_factor=1e-6 if mass_unit=='mg' and self.elements.variables['mass']['units']=='ug': mass_conversion_factor=1e-3 if mass_unit=='ug' and self.elements.variables['mass']['units']=='ug': mass_conversion_factor=1 if mass_unit=='kg' and self.elements.variables['mass']['units']=='ug': mass_conversion_factor=1e-9 time_conversion_factor = self.time_step_output.total_seconds() / (60*60) if time_unit=='seconds': time_conversion_factor = self.time_step_output.total_seconds() if time_unit=='minutes': time_conversion_factor = self.time_step_output.total_seconds() / 60 if time_unit=='hours': time_conversion_factor = self.time_step_output.total_seconds() / (60*60) if time_unit=='days': time_conversion_factor = self.time_step_output.total_seconds() / (24*60*60) for i in range(steps): bars[i]=[np.sum(mass[0][i,:]*(sp[0][i,:]==0))*mass_conversion_factor, np.sum(mass[0][i,:]*(sp[0][i,:]==1))*mass_conversion_factor, np.sum(mass[0][i,:]*(sp[0][i,:]==2))*mass_conversion_factor, np.sum(mass[0][i,:]*(sp[0][i,:]==3))*mass_conversion_factor, np.sum(mass[0][i,:]*(sp[0][i,:]==4))*mass_conversion_factor] bottom=np.zeros_like(bars[:,0]) if 'dissolved' in legend: ax.bar(np.arange(steps),bars[:,self.num_lmm],width=1.25,color='midnightblue') bottom=bars[:,self.num_lmm] print('dissolved' + ' : ' + str(bars[-1,self.num_lmm]) + mass_unit +' ('+ str(100*bars[-1,self.num_lmm]/np.sum(bars[-1,:]))+'%)') if 'DOC' in legend: ax.bar(np.arange(steps),bars[:,self.num_humcol],bottom=bottom,width=1.25,color='royalblue') bottom=bottom+bars[:,self.num_humcol] print('DOC' + ' : ' + str(bars[-1,self.num_humcol]) + mass_unit +' ('+ str(100*bars[-1,self.num_humcol]/np.sum(bars[-1,:]))+'%)') if 'SPM' in legend: ax.bar(np.arange(steps),bars[:,self.num_prev],bottom=bottom,width=1.25,color='palegreen') bottom=bottom+bars[:,self.num_prev] print('SPM' + ' : ' + str(bars[-1,self.num_prev]) + mass_unit +' ('+ str(100*bars[-1,self.num_prev]/np.sum(bars[-1,:]))+'%)') if 'sediment' in legend: ax.bar(np.arange(steps),bars[:,self.num_srev],bottom=bottom,width=1.25,color='orange') bottom=bottom+bars[:,self.num_srev] print('sediment' + ' : ' + str(bars[-1,self.num_srev]) + mass_unit +' ('+ str(100*bars[-1,self.num_srev]/np.sum(bars[-1,:]))+'%)') ax.legend(list(filter(None, legend))) ax.set_ylabel('mass (' + mass_unit + ')') if start_date is None: ax.axes.get_xaxis().set_ticklabels(np.round(ax.axes.get_xticks() * time_conversion_factor)) ax.set_xlabel('time (' + time_unit + ')') else: date_values = [datetime.strptime(start_date,"%Y-%m-%d") + i*self.time_step for i in range(steps)] # Set a fraction of datetime values as labels for the x-axis fraction = 24*7 # Show every second datetime value ax.set_xticks(np.arange(0, steps, fraction)) ax.set_xticklabels([date.strftime('%m-%d') for date in date_values[::fraction]]) ax.set_xlabel('time (month-day)') fig.show() if filename is not None: plt.savefig(filename, format=filename[-3:], transparent=True, bbox_inches="tight", dpi=300)