From 070b1ae67f0e653de6b6220f2e76c8d747e9c9e7 Mon Sep 17 00:00:00 2001 From: Magne Simonsen Date: Wed, 5 Jun 2024 11:39:44 +0200 Subject: [PATCH 1/6] Major upgrade - change isotope and transfer setup config - remove the species as config variables - remove radionuclide:slowly_fraction and radionuclide:irreversible_fraction config keywords - refine the initialization and update of particle size (distribution) - include min/max particle size as config variables - enable release of slowly reversible particles (new config variable) - set depth intervals used in conc file as config settings - initialize species and specie names according to pre-defined specie setup. - enable more isotopes, which impact the Kd value - pick particles/dissolved elements randomly during initialization - --- opendrift/models/radionuclides.py | 602 ++++++++++++++++-------------- 1 file changed, 322 insertions(+), 280 deletions(-) diff --git a/opendrift/models/radionuclides.py b/opendrift/models/radionuclides.py index 009e860eb..e167b61f7 100644 --- a/opendrift/models/radionuclides.py +++ b/opendrift/models/radionuclides.py @@ -42,14 +42,14 @@ class Radionuclide(Lagrangian3DArray): # ('transfer_rates1D', {'dtype':np.array(3, dtype=np.float32), # 'units': '1/s', # 'default': 0.}) - ('LMM_fraction', {'dtype':np.float32, - 'units':'', - 'default':0, - 'seed':False}), - ('particle_fraction', {'dtype':np.float32, - 'units':'', - 'default':0, - 'seed':False}) + # ('LMM_fraction', {'dtype':np.float32, + # 'units':'', + # 'default':0, + # 'seed':False}), + # ('particle_fraction', {'dtype':np.float32, + # 'units':'', + # 'default':0, + # 'seed':False}) ]) @@ -107,13 +107,6 @@ def __init__(self, *args, **kwargs): # TODO: descriptions and units must be added in config setting below self._add_config({ - 'radionuclide:transfer_setup': {'type': 'enum', - 'enum': ['Sandnesfj_Al','Bokna_137Cs', '137Cs_rev', 'custom'], 'default': '137Cs_rev', - 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, - 'radionuclide:slowly_fraction': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - 'radionuclide:irreversible_fraction': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'radionuclide:dissolved_diameter': {'type': 'float', 'default': 0, 'min': 0, 'max': 100e-6, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, @@ -123,6 +116,12 @@ def __init__(self, *args, **kwargs): 'radionuclide:particle_diameter_uncertainty': {'type': 'float', 'default': 1e-7, 'min': 0, 'max': 100e-6, 'units': 'm', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Standard deviation of particle size distribution.'}, + 'radionuclide:particle_diameter_minimum': {'type': 'float', 'default': 0.45e-6, + 'min': 0, 'max': 100e-6, 'units': 'm', + 'level': CONFIG_LEVEL_ADVANCED, 'description': 'Mimimum particle size.'}, + 'radionuclide:particle_diameter_maximum': {'type': 'float', 'default': 63.e-6, + 'min': 0, 'max': 100e-6, 'units': 'm', + 'level': CONFIG_LEVEL_ADVANCED, 'description': 'Maximum particle size.'}, 'radionuclide:particlesize_distribution': {'type':'enum', 'enum': ['normal','lognormal'], 'default':'normal', 'level':CONFIG_LEVEL_ADVANCED, @@ -133,35 +132,24 @@ def __init__(self, *args, **kwargs): 'seed:particle_fraction': {'type': 'float','default': 0.9, 'min': 0, 'max': 1, 'units': '1', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Fraction of initial discharge released as particle species'}, + 'seed:slowly_fraction': {'type': 'float','default': 0., + 'min': 0, 'max': 1, 'units': '1', + 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Fraction of PARTICLE discharge released as slowly reversible particle species'}, 'seed:total_release': {'type': 'float','default': 100.e9, 'min': 0, 'max': 1e36, 'units': 'Bq', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Total release (Bq)'}, - # Species - 'radionuclide:species:LMM': {'type': 'bool', 'default': True, - 'level': CONFIG_LEVEL_BASIC, 'description': 'Toggle LMM species'}, - 'radionuclide:species:LMMcation': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:LMManion': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Colloid': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Humic_colloid': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - 'radionuclide:species:Polymer': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - 'radionuclide:species:Particle_reversible': {'type': 'bool', 'default': True, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Particle_slowly_reversible': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Particle_irreversible': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - 'radionuclide:species:Sediment_reversible': {'type': 'bool', 'default': True, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Sediment_slowly_reversible': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Sediment_irreversible': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - # Transformations + + # ISOTOPES + 'radionuclide:isotope': {'type': 'enum', 'default': '137Cs', + 'enum': ['Al', '137Cs', '129I', '241Am', 'manual'], + 'level':CONFIG_LEVEL_ESSENTIAL, 'description':'Isotope'}, + 'radionuclide:specie_setup': {'type': 'enum', 'default': 'LMM + Rev', + 'enum': ['LMM + Rev', 'LMM + Rev + Slow rev', 'LMM + Rev + Irrev', + 'LMM + Rev + Slow rev + Irrev', 'LMM + Colloid + Rev'], + 'level':CONFIG_LEVEL_ESSENTIAL, 'description':'Species enabled'}, + + + # Transformation parameters 'radionuclide:transformations:Kd': {'type': 'float', 'default': 2.0, 'min': 0, 'max': 1e9, 'units': 'm3/kg', 'level': CONFIG_LEVEL_BASIC, 'description': ''}, @@ -205,7 +193,14 @@ def __init__(self, *args, **kwargs): 'radionuclide:sediment:resuspension_critvel': {'type': 'float', 'default': .01, 'min': 0, 'max': 1, 'units': 'm/s', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - }) + # OUTPUT config + 'radionuclide:output:depthintervals': {'type':'enum', + 'enum': ['-50.', '-25.', '-10.', '-5.', '-2.', '-1.', + '-25., -5.', '-50., -10., -1.', + '-10., -1.'], 'default':'-10.', + 'level':CONFIG_LEVEL_ESSENTIAL, + 'description':'Depth intervals for computation of concentration'}, + }) @@ -218,80 +213,100 @@ def prepare_run(self): logger.info( '{:>3} {}'.format( i, sp ) ) - logger.info( 'transfer setup: %s' % self.get_config('radionuclide:transfer_setup')) + logger.info( 'Isotope: %s' % self.get_config('radionuclide:isotope')) + logger.info( 'specie setup: %s' % self.get_config('radionuclide:specie_setup')) logger.info('nspecies: %s' % self.nspecies) logger.info('Transfer rates:\n %s' % self.transfer_rates) + tmp = self.get_config('radionuclide:output:depthintervals') + self.depthintervals = [float(item) for item in tmp.split(',')] + logger.info('DEPTH_INTERVALS: %s',self.depthintervals) + + + super(RadionuclideDrift, self).prepare_run() def init_species(self): - # Initialize specie types - if self.get_config('radionuclide:transfer_setup')=='Bokna_137Cs': - self.set_config('radionuclide:species:LMM',True) - self.set_config('radionuclide:species:Particle_reversible', True) - self.set_config('radionuclide:species:Particle_slowly_reversible', True) - self.set_config('radionuclide:species:Sediment_reversible', True) - self.set_config('radionuclide:species:Sediment_slowly_reversible', True) - elif self.get_config('radionuclide:transfer_setup')=='137Cs_rev': - self.set_config('radionuclide:species:LMM',True) - self.set_config('radionuclide:species:Particle_reversible', True) - self.set_config('radionuclide:species:Sediment_reversible', True) - elif self.get_config('radionuclide:transfer_setup')=='Sandnesfj_Al': - self.set_config('radionuclide:species:LMM', False) - self.set_config('radionuclide:species:LMMcation', True) - self.set_config('radionuclide:species:LMManion', True) - self.set_config('radionuclide:species:Humic_colloid', True) - self.set_config('radionuclide:species:Polymer', True) - self.set_config('radionuclide:species:Particle_reversible', True) - self.set_config('radionuclide:species:Sediment_reversible', True) - elif self.get_config('radionuclide:transfer_setup')=='custom': - # Do nothing, species must be set manually - pass - else: - logger.error('No valid transfer_setup {}'.format(self.get_config('radionuclide:transfer_setup'))) + + # Initialize slowly and irrev fractions to False, set True later if enabled + self.species_slowly_fraction = False + self.species_irreversible_fraction = False + # Initialize species, add enabled species to name_species list self.name_species=[] - if self.get_config('radionuclide:species:LMM'): + self.specie_setup = self.get_config('radionuclide:specie_setup') + + if 'lmm + rev' in self.specie_setup.casefold(): self.name_species.append('LMM') - if self.get_config('radionuclide:species:LMMcation'): + self.name_species.append('Particle reversible') + self.name_species.append('Sediment reversible') + else: + logger.error('No valid specie setup {}'.format(self.get_config('radionuclide:specie_setup'))) + + + if '+ slow rev' in self.specie_setup.casefold(): + self.name_species.append('Particle slowly reversible') + self.name_species.append('Sediment slowly reversible') + self.species_slowly_fraction = True + if '+ irrev' in self.specie_setup.casefold(): + self.name_species.append('Particle irreversible') + self.name_species.append('Sediment irreversible') + self.species_irreversible_fraction = True + if self.specie_setup.casefold() == 'lmm + colloid + rev': self.name_species.append('LMMcation') - if self.get_config('radionuclide:species:LMManion'): self.name_species.append('LMManion') - if self.get_config('radionuclide:species:Colloid'): - self.name_species.append('Colloid') - if self.get_config('radionuclide:species:Humic_colloid'): self.name_species.append('Humic colloid') - if self.get_config('radionuclide:species:Polymer'): self.name_species.append('Polymer') - if self.get_config('radionuclide:species:Particle_reversible'): self.name_species.append('Particle reversible') - if self.get_config('radionuclide:species:Particle_slowly_reversible'): - self.name_species.append('Particle slowly reversible') - if self.get_config('radionuclide:species:Particle_irreversible'): - self.name_species.append('Particle irreversible') - if self.get_config('radionuclide:species:Sediment_reversible'): self.name_species.append('Sediment reversible') - if self.get_config('radionuclide:species:Sediment_slowly_reversible'): - self.name_species.append('Sediment slowly reversible') - if self.get_config('radionuclide:species:Sediment_irreversible'): - self.name_species.append('Sediment irreversible') - if self.get_config('radionuclide:species:Sediment_slowly_reversible') and \ - self.get_config('radionuclide:species:Particle_slowly_reversible'): - self.set_config('radionuclide:slowly_fraction', True) - if self.get_config('radionuclide:species:Sediment_irreversible') and \ - self.get_config('radionuclide:species:Particle_irreversible'): - self.set_config('radionuclide: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 ) ) - 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 ) ) + + + + def set_init_diameter(self, num, idxs,diam): + """ Initialize diameter for particles, according to size distribution """ + + distr = self.get_config('radionuclide:particlesize_distribution') + + init_diam = np.zeros(num) + npart = len(idxs) + + mu = 0. + uncert = self.get_config('radionuclide:particle_diameter_uncertainty') + if diam>0: + uncert_ln = uncert/diam * 3. + else: + uncert_ln = 0. + + rng = np.random.default_rng() + + if distr=='normal': + truncs = rng.normal(mu, uncert, npart) + init_diam[idxs] = diam + truncs + + elif distr=='lognormal': + truncs = rng.lognormal(mu, uncert_ln, size=npart ) # * diameter + init_diam[idxs] = diam * truncs + else: + logger.error('Not available distribution: %s',distr) + + pmin = self.get_config('radionuclide:particle_diameter_minimum') + pmax = self.get_config('radionuclide:particle_diameter_maximum') + init_diam[idxs][init_diam[idxs]pmax] = pmax + + + return init_diam @@ -302,11 +317,9 @@ def init_species(self): def seed_elements(self, *args, **kwargs): self.init_species() - self.init_transfer_rates() - if 'number' in kwargs: num_elements = kwargs['number'] else: @@ -314,6 +327,8 @@ def seed_elements(self, *args, **kwargs): + + # Initialize species if 'specie' in kwargs: print('num_elements', num_elements) try: @@ -325,9 +340,6 @@ def seed_elements(self, *args, **kwargs): init_specie[:] = kwargs['specie'] kwargs['specie'] = init_specie[:] - - - else: # Set initial speciation @@ -341,7 +353,6 @@ def seed_elements(self, *args, **kwargs): else: lmm_frac = self.get_config('seed:LMM_fraction') - shift = int(num_elements * (1-particle_frac)) 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)) @@ -349,11 +360,23 @@ def seed_elements(self, *args, **kwargs): raise ValueError('Illegal specie fraction combination : ' + str(lmm_frac) + ' '+ str(particle_frac) ) init_specie = np.ones(num_elements, int) - if self.get_config('radionuclide:transfer_setup')=='Sandnesfj_Al': - init_specie[:shift] = self.num_lmmcation + dissolved = np.random.rand(num_elements)0 and 'Particle slowly reversible' in self.name_species: + ndiss = np.sum(dissolved) + prtidx = np.where(~dissolved)[0] + slowprt = np.random.rand(num_elements-ndiss)9} {:>3} {:24} '.format( np.sum(init_specie==i), i, sp ) ) + + if all(init_specie == self.num_srev): + print( 'ALL ELEMENTS ARE SEDIMENTS') + kwargs['z'] = 'seafloor' + # TODO SET MOVING = 0 + + elif any(init_specie == self.num_srev): + print( 'SOME ELEMENTS ARE SEDIMENTS') + exit() + + + # Set initial particle size according to speciation if 'diameter' in kwargs: diameter = kwargs['diameter'] else: diameter = self.get_config('radionuclide:particle_diameter') - init_diam =np.zeros(num_elements,float) - if self.get_config('radionuclide:particlesize_distribution')=='normal': - uncert = self.get_config('radionuclide:particle_diameter_uncertainty') - init_diam[init_specie==self.num_prev] = diameter - init_diam[init_specie==self.num_prev] += np.random.normal( - 0, uncert, sum(init_specie==self.num_prev)) - elif self.get_config('radionuclide:particlesize_distribution')=='lognormal': - std = 0.5 - mu = 0. - if self.get_config('radionuclide:species:Particle_reversible'): - init_diam[init_specie==self.num_prev] = \ - np.random.lognormal(mu, std, size=sum(init_specie==self.num_prev) ) * diameter - if self.get_config('radionuclide:species:Particle_slowly_reversible'): - init_diam[init_specie==self.num_psrev] = \ - np.random.lognormal(mu, std, size=sum(init_specie==self.num_psrev) ) * diameter - if self.get_config('radionuclide:species:Particle_irreversible'): - init_diam[init_specie==self.num_pirrev] = \ - np.random.lognormal(mu, std, size=sum(init_specie==self.num_pirrev) ) * diameter - - - init_diam[(init_specie==self.num_prev) & (init_diam<0.45e-6)] = 0.45e-6 - init_diam[(init_specie==self.num_prev) & (init_diam>63.5e-6)] = 63.5e-6 - if self.get_config('radionuclide:species:Particle_slowly_reversible'): - init_diam[(init_specie==self.num_psrev) & (init_diam<0.45e-6)] = 0.45e-6 - init_diam[(init_specie==self.num_psrev) & (init_diam>63.5e-6)] = 63.5e-6 - if self.get_config('radionuclide:species:Particle_irreversible'): - init_diam[(init_specie==self.num_pirrev) & (init_diam<0.45e-6)] = 0.45e-6 - init_diam[(init_specie==self.num_pirrev) & (init_diam>63.5e-6)] = 63.5e-6 - - logger.info('Min: {:.2e}, Max: {:.2e}, Numer of particles seeded: {}, at {} distribution'.format( - np.min(init_diam[init_diam>0.]), np.max(init_diam[init_diam>0.]), sum(init_diam>0.), - self.get_config('radionuclide:particlesize_distribution'))) + logger.info('PARTICLE DIAMETER %s',diameter) + + particles = [] + for i,sp in enumerate(self.name_species): + if 'particle' in sp.casefold(): + prt = np.where(init_specie==i)[0] + particles.extend(prt) + + init_diam = self.set_init_diameter(num_elements, particles, diameter) + + + try: + logger.info('Min: {:.2e}, Max: {:.2e}, Numer of particles seeded: {}, at {} distribution'.format( + np.min(init_diam[init_diam>0.]), np.max(init_diam[init_diam>0.]), sum(init_diam>0.), + self.get_config('radionuclide:particlesize_distribution'))) + except Exception: + logger.info('All diameters are 0 (min:{}, max: {})'.format( np.min(init_diam), np.max(init_diam)) ) kwargs['diameter'] = init_diam + + + # Set z to seabed for sediments + # sediments = [] + # for i,sp in enumerate(self.name_species): + # if 'sediment' in sp.casefold(): + # sed = np.where(init_specie==i)[0] + # sediments.extend(sed) + + # kwargs['z'] = 'seafloor' + # logger.debug('Seafloor is selected, neglecting z') +# init_z = np.ones(num_elements,dtype=float) +# init_z[:] = kwargs['z'] +# init_z[sediments] = -10000. #\ +# #-1.*self.environment.sea_floor_depth_below_sea_level[sediments] +# # self.elements.moving[sediments] = 0 +# kwargs['z'] = init_z super(RadionuclideDrift, self).seed_elements(*args, **kwargs) + def init_kd(self): + ''' Initialization of Kd value, dependent on simulated radionuclide + ''' + + # Values from IAEA (2004) + kd_values = { '137Cs' : 4.0e0, + '129I' : 7.0e-2, + '241Am' : 2.0e3, + 'Al' : None, + } + + + print('ISOTOPE',self.isotope) + + if self.isotope == 'manual': + self.kd = self.get_config('radionuclide:transformations:Kd') + return + + if not self.isotope in kd_values.keys(): + logger.error('No Kd value implemented for %s ', self.isotope) + + else: + self.kd = kd_values[self.isotope] + + + + + def init_transfer_rates(self): ''' Initialization of background values in the transfer rates 2D array. ''' - transfer_setup=self.get_config('radionuclide:transfer_setup') # logger.info( 'transfer setup: %s' % transfer_setup) + self.isotope = self.get_config('radionuclide:isotope') + + self.init_kd() + + logger.info('ISOTOPE %s',self.isotope) + logger.info('SPECIE SETUP %s',self.specie_setup) + logger.info('Kd: %s',self.kd) self.transfer_rates = np.zeros([self.nspecies,self.nspecies]) self.ntransformations = np.zeros([self.nspecies,self.nspecies]) - if transfer_setup == 'Bokna_137Cs': + + if 'lmm + rev' in self.specie_setup.casefold(): 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('radionuclide:transformations:Kd') + # Simpler version of Values from Simonsen et al (2019a) + # Only consider the reversible fraction + Kd = self.kd # depends on isotope Dc = self.get_config('radionuclide:transformations:Dc') - slow_coeff = self.get_config('radionuclide:transformations:slow_coeff') susp_mat = 1.e-3 # concentration of available suspended particulate matter (kg/m3) sedmixdepth = self.get_config('radionuclide:sediment:sedmixdepth') # sediment mixing depth (m) default_density = self.get_config('radionuclide:sediment:sediment_density') # default particle density (kg/m3) @@ -448,72 +518,38 @@ def init_transfer_rates(self): self.transfer_rates[self.num_lmm,self.num_srev] = \ Dc * Kd * sedmixdepth * default_density * (1.-poro) * f * phi / layer_thick self.transfer_rates[self.num_srev,self.num_lmm] = Dc * phi + + else: + logger.error('No transfer setup available') + + + if '+ slow rev' in self.specie_setup.casefold(): + + self.num_psrev = self.specie_name2num('Particle slowly reversible') + self.num_ssrev = self.specie_name2num('Sediment slowly reversible') + + slow_coeff = self.get_config('radionuclide:transformations:slow_coeff') + 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 + if '+ irrev' in self.specie_setup.casefold(): - elif transfer_setup == '137Cs_rev': + self.num_pirrev = self.specie_name2num('Particle irreversible') + self.num_sirrev = self.specie_name2num('Sediment irreversible') - self.num_lmm = self.specie_name2num('LMM') - self.num_prev = self.specie_name2num('Particle reversible') - self.num_srev = self.specie_name2num('Sediment reversible') + slow_coeff = self.get_config('radionuclide:transformations:slow_coeff') + + self.transfer_rates[self.num_ssrev,self.num_sirrev] = slow_coeff + self.transfer_rates[self.num_psrev,self.num_pirrev] = slow_coeff - # Simpler version of Values from Simonsen et al (2019a) - # Only consider the reversible fraction - Kd = self.get_config('radionuclide:transformations:Kd') - Dc = self.get_config('radionuclide:transformations:Dc') - susp_mat = 1.e-3 # concentration of available suspended particulate matter (kg/m3) - sedmixdepth = self.get_config('radionuclide:sediment:sedmixdepth') # sediment mixing depth (m) - default_density = self.get_config('radionuclide:sediment:sediment_density') # default particle density (kg/m3) - f = self.get_config('radionuclide:sediment:effective_fraction') # fraction of effective sorbents - phi = self.get_config('radionuclide:sediment:corr_factor') # sediment correction factor - poro = self.get_config('radionuclide:sediment:porosity') # sediment porosity - layer_thick = self.get_config('radionuclide:sediment:layer_thick') # thickness of seabed interaction layer (m) - self.transfer_rates[self.num_lmm,self.num_prev] = Dc * Kd * susp_mat - self.transfer_rates[self.num_prev,self.num_lmm] = Dc - self.transfer_rates[self.num_lmm,self.num_srev] = \ - Dc * Kd * sedmixdepth * default_density * (1.-poro) * f * phi / layer_thick - self.transfer_rates[self.num_srev,self.num_lmm] = Dc * phi - elif transfer_setup=='custom': - # Set of custom values for testing/development - - self.num_lmm = self.specie_name2num('LMM') - if self.get_config('radionuclide:species:Colloid'): - self.num_col = self.specie_name2num('Colloid') - if self.get_config('radionuclide:species:Particle_reversible'): - self.num_prev = self.specie_name2num('Particle reversible') - if self.get_config('radionuclide:species:Sediment_reversible'): - self.num_srev = self.specie_name2num('Sediment reversible') - if self.get_config('radionuclide:slowly_fraction'): - self.num_psrev = self.specie_name2num('Particle slowly reversible') - self.num_ssrev = self.specie_name2num('Sediment slowly reversible') - if self.get_config('radionuclide:irreversible_fraction'): - self.num_pirrev = self.specie_name2num('Particle irreversible') - self.num_sirrev = self.specie_name2num('Sediment irreversible') - - - if self.get_config('radionuclide: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('radionuclide:transformations:Dc') - if self.get_config('radionuclide: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('radionuclide:transformations:Dc') * self.get_config('radionuclide:sediment:corr_factor') -# self.transfer_rates[self.num_srev,self.num_lmm] = 5.e-6 - - if self.get_config('radionuclide: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': + if self.specie_setup.casefold() == 'lmm + colloid + rev' and \ + self.isotope.casefold() == 'al': # Use values from Simonsen et al (2019b) self.num_lmmanion = self.specie_name2num('LMManion') self.num_lmmcation = self.specie_name2num('LMMcation') @@ -576,10 +612,6 @@ def init_transfer_rates(self): - 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]): @@ -587,8 +619,6 @@ def init_transfer_rates(self): else: np.fill_diagonal(self.transfer_rates,0.) -# self.transfer_rates[:] = 0. -# print ('\n ###### \n IMPORTANT:: \n transfer rates are all set to zero for testing purposes! \n#### \n ') @@ -658,17 +688,18 @@ def update_terminal_velocity(self, Tprofiles=None, self.elements.terminal_velocity = W * self.elements.moving + + + + 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('radionuclide:transfer_setup') - if transfer_setup == 'Bokna_137Cs' or \ - transfer_setup=='custom' or \ - transfer_setup=='137Cs_rev': + if not self.specie_setup.casefold() == 'lmm + colloid + rev': self.elements.transfer_rates1D = self.transfer_rates[self.elements.specie,:] - if self.get_config('radionuclide:species:Sediment_reversible'): + if 'Sediment reversible' in self.name_species: # Only LMM radionuclides 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 @@ -677,8 +708,7 @@ def update_transfer_rates(self): self.elements.transfer_rates1D[(self.elements.specie == self.num_lmm) & (dist_to_seabed > interaction_thick), self.num_srev] = 0. - - if self.get_config('radionuclide:species:Particle_reversible'): + if 'Particle reversible' in self.name_species: # Modify particle adsorption according to local particle concentration # (LMM -> reversible particles) kktmp = self.elements.specie == self.num_lmm @@ -687,7 +717,7 @@ def update_transfer_rates(self): self.environment.conc3[kktmp] / 1.e-3 # self.environment.particle_conc[kktmp] / 1.e-3 - elif transfer_setup=='Sandnesfj_Al': + elif self.specie_setup.casefold() == 'lmm + colloid + rev' and self.isotope=='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,:] @@ -748,18 +778,16 @@ def update_speciation(self): - - def sorption_to_sediments(self,sp_in=None,sp_out=None): '''Update radionuclide properties when sorption to sediments occurs''' # Set z to local sea depth - if self.get_config('radionuclide:species:LMM'): + if 'LMM' in self.name_species: 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('radionuclide:species:LMMcation'): + if 'LMMcation' in self.name_species: 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 @@ -777,7 +805,7 @@ def desorption_from_sediments(self,sp_in=None,sp_out=None): std = self.get_config('radionuclide:sediment:desorption_depth_uncert') - if self.get_config('radionuclide:species:LMM'): + if 'LMM' in self.name_species: 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 @@ -785,7 +813,7 @@ def desorption_from_sediments(self,sp_in=None,sp_out=None): 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('radionuclide:species:LMMcation'): + if 'LMMcation' in self.name_species: 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 @@ -810,45 +838,40 @@ def update_radionuclide_diameter(self,sp_in=None,sp_out=None): dia_diss=self.get_config('radionuclide:dissolved_diameter') - # Transfer to reversible particles - self.elements.diameter[(sp_out==self.num_prev) & (sp_in!=self.num_prev)] = dia_part - std = self.get_config('radionuclide: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('radionuclide: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))) + # Transform to particle species + to_particles = np.where( (sp_out==self.num_prev) & (sp_in!=self.num_prev) ) [0] + if self.species_slowly_fraction: + to_particles = np.append(to_particles, np.where((sp_out==self.num_psrev) & (sp_in==self.num_ssrev))[0] ) + if self.species_irreversible_fraction: + to_particles = np.append(to_particles, np.where((sp_out==self.num_pirrev) & (sp_in==self.num_sirrev))[0] ) + + new_diam = self.set_init_diameter(self.num_elements_total(), to_particles, dia_part) + self.elements.diameter[to_particles] = new_diam[to_particles] + + + + + # Transform to LMM and colloid species + if 'LMM' in self.name_species: + to_diss = np.where( (sp_out==self.num_lmm) & (sp_in!=self.num_lmm) )[0] + elif 'LMMcation' and 'LMManion' in self.name_species: + to_diss = np.where( (sp_out==self.num_lmmcation) & (sp_in!=self.num_lmmcation) )[0] + to_diss = np.append( to_diss, np.where( (sp_out==self.num_lmmanion) & (sp_in!=self.num_lmmanion) )[0] ) + if 'Colloid' in self.name_species: + to_diss = np.append(to_diss, np.where( (sp_out==self.num_col) & (sp_in!=self.num_col) )[0] ) + if 'Humic_colloid' in self.name_species: + to_diss = np.append( to_diss, np.where( (sp_out==self.num_humcol) & (sp_in!=self.num_humcol) )[0] ) + if 'Polymer' in self.name_species: + to_diss = np.append( to_diss, np.where( (sp_out==self.num_polymer) & (sp_in!=self.num_polymer) )[0] ) + + new_diam = self.set_init_diameter(self.num_elements_total(), to_diss, dia_diss) + self.elements.diameter[to_diss] = new_diam[to_diss] + + - # Transfer to irreversible particles - if self.get_config('radionuclide: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('radionuclide:species:LMM'): - self.elements.diameter[(sp_out==self.num_lmm) & (sp_in!=self.num_lmm)] = dia_diss - if self.get_config('radionuclide:species:LMManion'): - self.elements.diameter[(sp_out==self.num_lmmanion) & (sp_in!=self.num_lmmanion)] = dia_diss - if self.get_config('radionuclide:species:LMMcation'): - self.elements.diameter[(sp_out==self.num_lmmcation) & (sp_in!=self.num_lmmcation)] = dia_diss - # Transfer to colloids - if self.get_config('radionuclide:species:Colloid'): - self.elements.diameter[(sp_out==self.num_col) & (sp_in!=self.num_col)] = dia_diss - if self.get_config('radionuclide:species:Humic_colloid'): - self.elements.diameter[(sp_out==self.num_humcol) & (sp_in!=self.num_humcol)] = dia_diss - if self.get_config('radionuclide:species:Polymer'): - self.elements.diameter[(sp_out==self.num_polymer) & (sp_in!=self.num_polymer)] = dia_diss @@ -856,10 +879,14 @@ def update_radionuclide_diameter(self,sp_in=None,sp_out=None): def bottom_interaction(self,Zmin=None): ''' Change speciation of radionuclides that reach bottom due to settling. particle specie -> sediment specie ''' - if not ((self.get_config('radionuclide:species:Particle_reversible')) & - (self.get_config('radionuclide:species:Sediment_reversible')) or - (self.get_config('radionuclide:slowly_fraction')) or - (self.get_config('radionuclide:irreversible_fraction'))): + + if not ( ('Particle reversible' in self.name_species) & + ('Sediment reversible' in self.name_species) or + self.species_slowly_fraction or self.species_irreversible_fraction): + print('RETURN: ', ('Particle reversible' in self.name_species), + ('Sediment reversible' in self.name_species), + self.species_slowly_fraction, + self.species_irreversible_fraction) return bottom = np.array(np.where(self.elements.z <= Zmin)[0]) @@ -867,12 +894,15 @@ def bottom_interaction(self,Zmin=None): 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('radionuclide:slowly_fraction'): +# self.elements.moving[bottom] = 0 + + if self.species_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('radionuclide:irreversible_fraction'): + + if self.species_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) @@ -886,8 +916,11 @@ def resuspension(self): Sediment species -> Particle specie """ # Exit function if particles and sediments not are present - if not ((self.get_config('radionuclide:species:Particle_reversible')) & - (self.get_config('radionuclide:species:Sediment_reversible'))): + # if not ((self.get_config('radionuclide:species:Particle_reversible')) & + # (self.get_config('radionuclide:species:Sediment_reversible'))): + # return + if not ('Particle reversible' in self.name_species and 'Sediment reversible' in self.name_species): + logger.info('No sediments and particles present...') return specie_in = self.elements.specie.copy() @@ -918,10 +951,12 @@ def resuspension(self): 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('radionuclide:slowly_fraction'): +# if self.get_config('radionuclide:slowly_fraction'): + if self.species_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('radionuclide:irreversible_fraction'): + if self.species_irreversible_fraction: +# if self.get_config('radionuclide: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 @@ -958,7 +993,9 @@ def update(self): logger.info('Speciation: {} {}'.format([sum(self.elements.specie==ii) for ii in range(self.nspecies)],self.name_species)) + #self.elements.moving[self.elements.z <= -1*self.environment.sea_floor_depth_below_sea_level] = 0. + # Horizontal advection self.advect_ocean_current() @@ -973,8 +1010,6 @@ def update(self): - - # ################ # POSTPROCESSING @@ -1407,29 +1442,36 @@ def gui_postproc(self): logger.info ('{:32}: {:>6}'.format(sp,sum(self.elements.specie==isp))) homefolder = expanduser("~") - filename = homefolder+'/conc_radio.nc' + filename = homefolder+'/conc_radio_gui.nc' self.guipp_saveconcfile(filename) + """ zlayer = [-1,-2] self.guipp_plotandsaveconc(filename=filename, outfilename=homefolder+'/radio_plots/RadioConc', zlayers=zlayer, specie= ['Total', 'LMM', 'Particle reversible'], ) + """ - def guipp_saveconcfile(self,filename='none'): + def guipp_saveconcfile(self,filename='none', pixelsize_m=200., reader_sea_depth=None): - for readerName in self.readers: - reader = self.readers[readerName] - if 'sea_floor_depth_below_sea_level' in reader.variables: - break + + 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: + for readerName in self.env.readers: + reader = self.env.readers[readerName] + if 'sea_floor_depth_below_sea_level' in reader.variables: + break self.conc_lon = reader.get_variables('longitude', @@ -1445,8 +1487,8 @@ def guipp_saveconcfile(self,filename='none'): y=[reader.ymin,reader.ymax])['sea_floor_depth_below_sea_level'][:] - self.write_netcdf_radionuclide_density_map(filename, pixelsize_m=200., - zlevels=[-25, -2.], + self.write_netcdf_radionuclide_density_map(filename, pixelsize_m=pixelsize_m, + zlevels=self.depthintervals, activity_unit='Bq', horizontal_smoothing=True, smoothing_cells=1, @@ -1495,14 +1537,14 @@ def guipp_plotandsaveconc(self, filename=None, outfilename=None, zlayers=None, t else: specie_arr = specie - for readerName in self.readers: - reader = self.readers[readerName] - if 'sea_floor_depth_below_sea_level' in reader.variables: - break + # for readerName in self.readers: + # reader = self.readers[readerName] + # if 'sea_floor_depth_below_sea_level' in reader.variables: + # break - gtopo = reader.get_variables('sea_floor_depth_below_sea_level', - x=[reader.xmin,reader.xmax], - y=[reader.ymin,reader.ymax])['sea_floor_depth_below_sea_level'][:] + # gtopo = reader.get_variables('sea_floor_depth_below_sea_level', + # x=[reader.xmin,reader.xmax], + # y=[reader.ymin,reader.ymax])['sea_floor_depth_below_sea_level'][:] From 64622afdc6a76c107039228fa9767e9960b60916 Mon Sep 17 00:00:00 2001 From: Magne Simonsen Date: Wed, 5 Jun 2024 11:39:44 +0200 Subject: [PATCH 2/6] Major upgrade - change isotope and transfer setup config - remove the species as config variables - remove radionuclide:slowly_fraction and radionuclide:irreversible_fraction config keywords - refine the initialization and update of particle size (distribution) - include min/max particle size as config variables - enable release of slowly reversible particles (new config variable) - set depth intervals used in conc file as config settings - initialize species and specie names according to pre-defined specie setup. - enable more isotopes, which impact the Kd value - pick particles/dissolved elements randomly during initialization - --- opendrift/models/radionuclides.py | 602 ++++++++++++++++-------------- 1 file changed, 322 insertions(+), 280 deletions(-) diff --git a/opendrift/models/radionuclides.py b/opendrift/models/radionuclides.py index 009e860eb..e167b61f7 100644 --- a/opendrift/models/radionuclides.py +++ b/opendrift/models/radionuclides.py @@ -42,14 +42,14 @@ class Radionuclide(Lagrangian3DArray): # ('transfer_rates1D', {'dtype':np.array(3, dtype=np.float32), # 'units': '1/s', # 'default': 0.}) - ('LMM_fraction', {'dtype':np.float32, - 'units':'', - 'default':0, - 'seed':False}), - ('particle_fraction', {'dtype':np.float32, - 'units':'', - 'default':0, - 'seed':False}) + # ('LMM_fraction', {'dtype':np.float32, + # 'units':'', + # 'default':0, + # 'seed':False}), + # ('particle_fraction', {'dtype':np.float32, + # 'units':'', + # 'default':0, + # 'seed':False}) ]) @@ -107,13 +107,6 @@ def __init__(self, *args, **kwargs): # TODO: descriptions and units must be added in config setting below self._add_config({ - 'radionuclide:transfer_setup': {'type': 'enum', - 'enum': ['Sandnesfj_Al','Bokna_137Cs', '137Cs_rev', 'custom'], 'default': '137Cs_rev', - 'level': CONFIG_LEVEL_ESSENTIAL, 'description': ''}, - 'radionuclide:slowly_fraction': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - 'radionuclide:irreversible_fraction': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'radionuclide:dissolved_diameter': {'type': 'float', 'default': 0, 'min': 0, 'max': 100e-6, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, @@ -123,6 +116,12 @@ def __init__(self, *args, **kwargs): 'radionuclide:particle_diameter_uncertainty': {'type': 'float', 'default': 1e-7, 'min': 0, 'max': 100e-6, 'units': 'm', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Standard deviation of particle size distribution.'}, + 'radionuclide:particle_diameter_minimum': {'type': 'float', 'default': 0.45e-6, + 'min': 0, 'max': 100e-6, 'units': 'm', + 'level': CONFIG_LEVEL_ADVANCED, 'description': 'Mimimum particle size.'}, + 'radionuclide:particle_diameter_maximum': {'type': 'float', 'default': 63.e-6, + 'min': 0, 'max': 100e-6, 'units': 'm', + 'level': CONFIG_LEVEL_ADVANCED, 'description': 'Maximum particle size.'}, 'radionuclide:particlesize_distribution': {'type':'enum', 'enum': ['normal','lognormal'], 'default':'normal', 'level':CONFIG_LEVEL_ADVANCED, @@ -133,35 +132,24 @@ def __init__(self, *args, **kwargs): 'seed:particle_fraction': {'type': 'float','default': 0.9, 'min': 0, 'max': 1, 'units': '1', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Fraction of initial discharge released as particle species'}, + 'seed:slowly_fraction': {'type': 'float','default': 0., + 'min': 0, 'max': 1, 'units': '1', + 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Fraction of PARTICLE discharge released as slowly reversible particle species'}, 'seed:total_release': {'type': 'float','default': 100.e9, 'min': 0, 'max': 1e36, 'units': 'Bq', 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Total release (Bq)'}, - # Species - 'radionuclide:species:LMM': {'type': 'bool', 'default': True, - 'level': CONFIG_LEVEL_BASIC, 'description': 'Toggle LMM species'}, - 'radionuclide:species:LMMcation': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:LMManion': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Colloid': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Humic_colloid': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - 'radionuclide:species:Polymer': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - 'radionuclide:species:Particle_reversible': {'type': 'bool', 'default': True, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Particle_slowly_reversible': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Particle_irreversible': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - 'radionuclide:species:Sediment_reversible': {'type': 'bool', 'default': True, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Sediment_slowly_reversible': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_BASIC, 'description': ''}, - 'radionuclide:species:Sediment_irreversible': {'type': 'bool', 'default': False, - 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - # Transformations + + # ISOTOPES + 'radionuclide:isotope': {'type': 'enum', 'default': '137Cs', + 'enum': ['Al', '137Cs', '129I', '241Am', 'manual'], + 'level':CONFIG_LEVEL_ESSENTIAL, 'description':'Isotope'}, + 'radionuclide:specie_setup': {'type': 'enum', 'default': 'LMM + Rev', + 'enum': ['LMM + Rev', 'LMM + Rev + Slow rev', 'LMM + Rev + Irrev', + 'LMM + Rev + Slow rev + Irrev', 'LMM + Colloid + Rev'], + 'level':CONFIG_LEVEL_ESSENTIAL, 'description':'Species enabled'}, + + + # Transformation parameters 'radionuclide:transformations:Kd': {'type': 'float', 'default': 2.0, 'min': 0, 'max': 1e9, 'units': 'm3/kg', 'level': CONFIG_LEVEL_BASIC, 'description': ''}, @@ -205,7 +193,14 @@ def __init__(self, *args, **kwargs): 'radionuclide:sediment:resuspension_critvel': {'type': 'float', 'default': .01, 'min': 0, 'max': 1, 'units': 'm/s', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, - }) + # OUTPUT config + 'radionuclide:output:depthintervals': {'type':'enum', + 'enum': ['-50.', '-25.', '-10.', '-5.', '-2.', '-1.', + '-25., -5.', '-50., -10., -1.', + '-10., -1.'], 'default':'-10.', + 'level':CONFIG_LEVEL_ESSENTIAL, + 'description':'Depth intervals for computation of concentration'}, + }) @@ -218,80 +213,100 @@ def prepare_run(self): logger.info( '{:>3} {}'.format( i, sp ) ) - logger.info( 'transfer setup: %s' % self.get_config('radionuclide:transfer_setup')) + logger.info( 'Isotope: %s' % self.get_config('radionuclide:isotope')) + logger.info( 'specie setup: %s' % self.get_config('radionuclide:specie_setup')) logger.info('nspecies: %s' % self.nspecies) logger.info('Transfer rates:\n %s' % self.transfer_rates) + tmp = self.get_config('radionuclide:output:depthintervals') + self.depthintervals = [float(item) for item in tmp.split(',')] + logger.info('DEPTH_INTERVALS: %s',self.depthintervals) + + + super(RadionuclideDrift, self).prepare_run() def init_species(self): - # Initialize specie types - if self.get_config('radionuclide:transfer_setup')=='Bokna_137Cs': - self.set_config('radionuclide:species:LMM',True) - self.set_config('radionuclide:species:Particle_reversible', True) - self.set_config('radionuclide:species:Particle_slowly_reversible', True) - self.set_config('radionuclide:species:Sediment_reversible', True) - self.set_config('radionuclide:species:Sediment_slowly_reversible', True) - elif self.get_config('radionuclide:transfer_setup')=='137Cs_rev': - self.set_config('radionuclide:species:LMM',True) - self.set_config('radionuclide:species:Particle_reversible', True) - self.set_config('radionuclide:species:Sediment_reversible', True) - elif self.get_config('radionuclide:transfer_setup')=='Sandnesfj_Al': - self.set_config('radionuclide:species:LMM', False) - self.set_config('radionuclide:species:LMMcation', True) - self.set_config('radionuclide:species:LMManion', True) - self.set_config('radionuclide:species:Humic_colloid', True) - self.set_config('radionuclide:species:Polymer', True) - self.set_config('radionuclide:species:Particle_reversible', True) - self.set_config('radionuclide:species:Sediment_reversible', True) - elif self.get_config('radionuclide:transfer_setup')=='custom': - # Do nothing, species must be set manually - pass - else: - logger.error('No valid transfer_setup {}'.format(self.get_config('radionuclide:transfer_setup'))) + + # Initialize slowly and irrev fractions to False, set True later if enabled + self.species_slowly_fraction = False + self.species_irreversible_fraction = False + # Initialize species, add enabled species to name_species list self.name_species=[] - if self.get_config('radionuclide:species:LMM'): + self.specie_setup = self.get_config('radionuclide:specie_setup') + + if 'lmm + rev' in self.specie_setup.casefold(): self.name_species.append('LMM') - if self.get_config('radionuclide:species:LMMcation'): + self.name_species.append('Particle reversible') + self.name_species.append('Sediment reversible') + else: + logger.error('No valid specie setup {}'.format(self.get_config('radionuclide:specie_setup'))) + + + if '+ slow rev' in self.specie_setup.casefold(): + self.name_species.append('Particle slowly reversible') + self.name_species.append('Sediment slowly reversible') + self.species_slowly_fraction = True + if '+ irrev' in self.specie_setup.casefold(): + self.name_species.append('Particle irreversible') + self.name_species.append('Sediment irreversible') + self.species_irreversible_fraction = True + if self.specie_setup.casefold() == 'lmm + colloid + rev': self.name_species.append('LMMcation') - if self.get_config('radionuclide:species:LMManion'): self.name_species.append('LMManion') - if self.get_config('radionuclide:species:Colloid'): - self.name_species.append('Colloid') - if self.get_config('radionuclide:species:Humic_colloid'): self.name_species.append('Humic colloid') - if self.get_config('radionuclide:species:Polymer'): self.name_species.append('Polymer') - if self.get_config('radionuclide:species:Particle_reversible'): self.name_species.append('Particle reversible') - if self.get_config('radionuclide:species:Particle_slowly_reversible'): - self.name_species.append('Particle slowly reversible') - if self.get_config('radionuclide:species:Particle_irreversible'): - self.name_species.append('Particle irreversible') - if self.get_config('radionuclide:species:Sediment_reversible'): self.name_species.append('Sediment reversible') - if self.get_config('radionuclide:species:Sediment_slowly_reversible'): - self.name_species.append('Sediment slowly reversible') - if self.get_config('radionuclide:species:Sediment_irreversible'): - self.name_species.append('Sediment irreversible') - if self.get_config('radionuclide:species:Sediment_slowly_reversible') and \ - self.get_config('radionuclide:species:Particle_slowly_reversible'): - self.set_config('radionuclide:slowly_fraction', True) - if self.get_config('radionuclide:species:Sediment_irreversible') and \ - self.get_config('radionuclide:species:Particle_irreversible'): - self.set_config('radionuclide: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 ) ) - 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 ) ) + + + + def set_init_diameter(self, num, idxs,diam): + """ Initialize diameter for particles, according to size distribution """ + + distr = self.get_config('radionuclide:particlesize_distribution') + + init_diam = np.zeros(num) + npart = len(idxs) + + mu = 0. + uncert = self.get_config('radionuclide:particle_diameter_uncertainty') + if diam>0: + uncert_ln = uncert/diam * 3. + else: + uncert_ln = 0. + + rng = np.random.default_rng() + + if distr=='normal': + truncs = rng.normal(mu, uncert, npart) + init_diam[idxs] = diam + truncs + + elif distr=='lognormal': + truncs = rng.lognormal(mu, uncert_ln, size=npart ) # * diameter + init_diam[idxs] = diam * truncs + else: + logger.error('Not available distribution: %s',distr) + + pmin = self.get_config('radionuclide:particle_diameter_minimum') + pmax = self.get_config('radionuclide:particle_diameter_maximum') + init_diam[idxs][init_diam[idxs]pmax] = pmax + + + return init_diam @@ -302,11 +317,9 @@ def init_species(self): def seed_elements(self, *args, **kwargs): self.init_species() - self.init_transfer_rates() - if 'number' in kwargs: num_elements = kwargs['number'] else: @@ -314,6 +327,8 @@ def seed_elements(self, *args, **kwargs): + + # Initialize species if 'specie' in kwargs: print('num_elements', num_elements) try: @@ -325,9 +340,6 @@ def seed_elements(self, *args, **kwargs): init_specie[:] = kwargs['specie'] kwargs['specie'] = init_specie[:] - - - else: # Set initial speciation @@ -341,7 +353,6 @@ def seed_elements(self, *args, **kwargs): else: lmm_frac = self.get_config('seed:LMM_fraction') - shift = int(num_elements * (1-particle_frac)) 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)) @@ -349,11 +360,23 @@ def seed_elements(self, *args, **kwargs): raise ValueError('Illegal specie fraction combination : ' + str(lmm_frac) + ' '+ str(particle_frac) ) init_specie = np.ones(num_elements, int) - if self.get_config('radionuclide:transfer_setup')=='Sandnesfj_Al': - init_specie[:shift] = self.num_lmmcation + dissolved = np.random.rand(num_elements)0 and 'Particle slowly reversible' in self.name_species: + ndiss = np.sum(dissolved) + prtidx = np.where(~dissolved)[0] + slowprt = np.random.rand(num_elements-ndiss)9} {:>3} {:24} '.format( np.sum(init_specie==i), i, sp ) ) + + if all(init_specie == self.num_srev): + print( 'ALL ELEMENTS ARE SEDIMENTS') + kwargs['z'] = 'seafloor' + # TODO SET MOVING = 0 + + elif any(init_specie == self.num_srev): + print( 'SOME ELEMENTS ARE SEDIMENTS') + exit() + + + # Set initial particle size according to speciation if 'diameter' in kwargs: diameter = kwargs['diameter'] else: diameter = self.get_config('radionuclide:particle_diameter') - init_diam =np.zeros(num_elements,float) - if self.get_config('radionuclide:particlesize_distribution')=='normal': - uncert = self.get_config('radionuclide:particle_diameter_uncertainty') - init_diam[init_specie==self.num_prev] = diameter - init_diam[init_specie==self.num_prev] += np.random.normal( - 0, uncert, sum(init_specie==self.num_prev)) - elif self.get_config('radionuclide:particlesize_distribution')=='lognormal': - std = 0.5 - mu = 0. - if self.get_config('radionuclide:species:Particle_reversible'): - init_diam[init_specie==self.num_prev] = \ - np.random.lognormal(mu, std, size=sum(init_specie==self.num_prev) ) * diameter - if self.get_config('radionuclide:species:Particle_slowly_reversible'): - init_diam[init_specie==self.num_psrev] = \ - np.random.lognormal(mu, std, size=sum(init_specie==self.num_psrev) ) * diameter - if self.get_config('radionuclide:species:Particle_irreversible'): - init_diam[init_specie==self.num_pirrev] = \ - np.random.lognormal(mu, std, size=sum(init_specie==self.num_pirrev) ) * diameter - - - init_diam[(init_specie==self.num_prev) & (init_diam<0.45e-6)] = 0.45e-6 - init_diam[(init_specie==self.num_prev) & (init_diam>63.5e-6)] = 63.5e-6 - if self.get_config('radionuclide:species:Particle_slowly_reversible'): - init_diam[(init_specie==self.num_psrev) & (init_diam<0.45e-6)] = 0.45e-6 - init_diam[(init_specie==self.num_psrev) & (init_diam>63.5e-6)] = 63.5e-6 - if self.get_config('radionuclide:species:Particle_irreversible'): - init_diam[(init_specie==self.num_pirrev) & (init_diam<0.45e-6)] = 0.45e-6 - init_diam[(init_specie==self.num_pirrev) & (init_diam>63.5e-6)] = 63.5e-6 - - logger.info('Min: {:.2e}, Max: {:.2e}, Numer of particles seeded: {}, at {} distribution'.format( - np.min(init_diam[init_diam>0.]), np.max(init_diam[init_diam>0.]), sum(init_diam>0.), - self.get_config('radionuclide:particlesize_distribution'))) + logger.info('PARTICLE DIAMETER %s',diameter) + + particles = [] + for i,sp in enumerate(self.name_species): + if 'particle' in sp.casefold(): + prt = np.where(init_specie==i)[0] + particles.extend(prt) + + init_diam = self.set_init_diameter(num_elements, particles, diameter) + + + try: + logger.info('Min: {:.2e}, Max: {:.2e}, Numer of particles seeded: {}, at {} distribution'.format( + np.min(init_diam[init_diam>0.]), np.max(init_diam[init_diam>0.]), sum(init_diam>0.), + self.get_config('radionuclide:particlesize_distribution'))) + except Exception: + logger.info('All diameters are 0 (min:{}, max: {})'.format( np.min(init_diam), np.max(init_diam)) ) kwargs['diameter'] = init_diam + + + # Set z to seabed for sediments + # sediments = [] + # for i,sp in enumerate(self.name_species): + # if 'sediment' in sp.casefold(): + # sed = np.where(init_specie==i)[0] + # sediments.extend(sed) + + # kwargs['z'] = 'seafloor' + # logger.debug('Seafloor is selected, neglecting z') +# init_z = np.ones(num_elements,dtype=float) +# init_z[:] = kwargs['z'] +# init_z[sediments] = -10000. #\ +# #-1.*self.environment.sea_floor_depth_below_sea_level[sediments] +# # self.elements.moving[sediments] = 0 +# kwargs['z'] = init_z super(RadionuclideDrift, self).seed_elements(*args, **kwargs) + def init_kd(self): + ''' Initialization of Kd value, dependent on simulated radionuclide + ''' + + # Values from IAEA (2004) + kd_values = { '137Cs' : 4.0e0, + '129I' : 7.0e-2, + '241Am' : 2.0e3, + 'Al' : None, + } + + + print('ISOTOPE',self.isotope) + + if self.isotope == 'manual': + self.kd = self.get_config('radionuclide:transformations:Kd') + return + + if not self.isotope in kd_values.keys(): + logger.error('No Kd value implemented for %s ', self.isotope) + + else: + self.kd = kd_values[self.isotope] + + + + + def init_transfer_rates(self): ''' Initialization of background values in the transfer rates 2D array. ''' - transfer_setup=self.get_config('radionuclide:transfer_setup') # logger.info( 'transfer setup: %s' % transfer_setup) + self.isotope = self.get_config('radionuclide:isotope') + + self.init_kd() + + logger.info('ISOTOPE %s',self.isotope) + logger.info('SPECIE SETUP %s',self.specie_setup) + logger.info('Kd: %s',self.kd) self.transfer_rates = np.zeros([self.nspecies,self.nspecies]) self.ntransformations = np.zeros([self.nspecies,self.nspecies]) - if transfer_setup == 'Bokna_137Cs': + + if 'lmm + rev' in self.specie_setup.casefold(): 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('radionuclide:transformations:Kd') + # Simpler version of Values from Simonsen et al (2019a) + # Only consider the reversible fraction + Kd = self.kd # depends on isotope Dc = self.get_config('radionuclide:transformations:Dc') - slow_coeff = self.get_config('radionuclide:transformations:slow_coeff') susp_mat = 1.e-3 # concentration of available suspended particulate matter (kg/m3) sedmixdepth = self.get_config('radionuclide:sediment:sedmixdepth') # sediment mixing depth (m) default_density = self.get_config('radionuclide:sediment:sediment_density') # default particle density (kg/m3) @@ -448,72 +518,38 @@ def init_transfer_rates(self): self.transfer_rates[self.num_lmm,self.num_srev] = \ Dc * Kd * sedmixdepth * default_density * (1.-poro) * f * phi / layer_thick self.transfer_rates[self.num_srev,self.num_lmm] = Dc * phi + + else: + logger.error('No transfer setup available') + + + if '+ slow rev' in self.specie_setup.casefold(): + + self.num_psrev = self.specie_name2num('Particle slowly reversible') + self.num_ssrev = self.specie_name2num('Sediment slowly reversible') + + slow_coeff = self.get_config('radionuclide:transformations:slow_coeff') + 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 + if '+ irrev' in self.specie_setup.casefold(): - elif transfer_setup == '137Cs_rev': + self.num_pirrev = self.specie_name2num('Particle irreversible') + self.num_sirrev = self.specie_name2num('Sediment irreversible') - self.num_lmm = self.specie_name2num('LMM') - self.num_prev = self.specie_name2num('Particle reversible') - self.num_srev = self.specie_name2num('Sediment reversible') + slow_coeff = self.get_config('radionuclide:transformations:slow_coeff') + + self.transfer_rates[self.num_ssrev,self.num_sirrev] = slow_coeff + self.transfer_rates[self.num_psrev,self.num_pirrev] = slow_coeff - # Simpler version of Values from Simonsen et al (2019a) - # Only consider the reversible fraction - Kd = self.get_config('radionuclide:transformations:Kd') - Dc = self.get_config('radionuclide:transformations:Dc') - susp_mat = 1.e-3 # concentration of available suspended particulate matter (kg/m3) - sedmixdepth = self.get_config('radionuclide:sediment:sedmixdepth') # sediment mixing depth (m) - default_density = self.get_config('radionuclide:sediment:sediment_density') # default particle density (kg/m3) - f = self.get_config('radionuclide:sediment:effective_fraction') # fraction of effective sorbents - phi = self.get_config('radionuclide:sediment:corr_factor') # sediment correction factor - poro = self.get_config('radionuclide:sediment:porosity') # sediment porosity - layer_thick = self.get_config('radionuclide:sediment:layer_thick') # thickness of seabed interaction layer (m) - self.transfer_rates[self.num_lmm,self.num_prev] = Dc * Kd * susp_mat - self.transfer_rates[self.num_prev,self.num_lmm] = Dc - self.transfer_rates[self.num_lmm,self.num_srev] = \ - Dc * Kd * sedmixdepth * default_density * (1.-poro) * f * phi / layer_thick - self.transfer_rates[self.num_srev,self.num_lmm] = Dc * phi - elif transfer_setup=='custom': - # Set of custom values for testing/development - - self.num_lmm = self.specie_name2num('LMM') - if self.get_config('radionuclide:species:Colloid'): - self.num_col = self.specie_name2num('Colloid') - if self.get_config('radionuclide:species:Particle_reversible'): - self.num_prev = self.specie_name2num('Particle reversible') - if self.get_config('radionuclide:species:Sediment_reversible'): - self.num_srev = self.specie_name2num('Sediment reversible') - if self.get_config('radionuclide:slowly_fraction'): - self.num_psrev = self.specie_name2num('Particle slowly reversible') - self.num_ssrev = self.specie_name2num('Sediment slowly reversible') - if self.get_config('radionuclide:irreversible_fraction'): - self.num_pirrev = self.specie_name2num('Particle irreversible') - self.num_sirrev = self.specie_name2num('Sediment irreversible') - - - if self.get_config('radionuclide: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('radionuclide:transformations:Dc') - if self.get_config('radionuclide: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('radionuclide:transformations:Dc') * self.get_config('radionuclide:sediment:corr_factor') -# self.transfer_rates[self.num_srev,self.num_lmm] = 5.e-6 - - if self.get_config('radionuclide: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': + if self.specie_setup.casefold() == 'lmm + colloid + rev' and \ + self.isotope.casefold() == 'al': # Use values from Simonsen et al (2019b) self.num_lmmanion = self.specie_name2num('LMManion') self.num_lmmcation = self.specie_name2num('LMMcation') @@ -576,10 +612,6 @@ def init_transfer_rates(self): - 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]): @@ -587,8 +619,6 @@ def init_transfer_rates(self): else: np.fill_diagonal(self.transfer_rates,0.) -# self.transfer_rates[:] = 0. -# print ('\n ###### \n IMPORTANT:: \n transfer rates are all set to zero for testing purposes! \n#### \n ') @@ -658,17 +688,18 @@ def update_terminal_velocity(self, Tprofiles=None, self.elements.terminal_velocity = W * self.elements.moving + + + + 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('radionuclide:transfer_setup') - if transfer_setup == 'Bokna_137Cs' or \ - transfer_setup=='custom' or \ - transfer_setup=='137Cs_rev': + if not self.specie_setup.casefold() == 'lmm + colloid + rev': self.elements.transfer_rates1D = self.transfer_rates[self.elements.specie,:] - if self.get_config('radionuclide:species:Sediment_reversible'): + if 'Sediment reversible' in self.name_species: # Only LMM radionuclides 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 @@ -677,8 +708,7 @@ def update_transfer_rates(self): self.elements.transfer_rates1D[(self.elements.specie == self.num_lmm) & (dist_to_seabed > interaction_thick), self.num_srev] = 0. - - if self.get_config('radionuclide:species:Particle_reversible'): + if 'Particle reversible' in self.name_species: # Modify particle adsorption according to local particle concentration # (LMM -> reversible particles) kktmp = self.elements.specie == self.num_lmm @@ -687,7 +717,7 @@ def update_transfer_rates(self): self.environment.conc3[kktmp] / 1.e-3 # self.environment.particle_conc[kktmp] / 1.e-3 - elif transfer_setup=='Sandnesfj_Al': + elif self.specie_setup.casefold() == 'lmm + colloid + rev' and self.isotope=='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,:] @@ -748,18 +778,16 @@ def update_speciation(self): - - def sorption_to_sediments(self,sp_in=None,sp_out=None): '''Update radionuclide properties when sorption to sediments occurs''' # Set z to local sea depth - if self.get_config('radionuclide:species:LMM'): + if 'LMM' in self.name_species: 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('radionuclide:species:LMMcation'): + if 'LMMcation' in self.name_species: 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 @@ -777,7 +805,7 @@ def desorption_from_sediments(self,sp_in=None,sp_out=None): std = self.get_config('radionuclide:sediment:desorption_depth_uncert') - if self.get_config('radionuclide:species:LMM'): + if 'LMM' in self.name_species: 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 @@ -785,7 +813,7 @@ def desorption_from_sediments(self,sp_in=None,sp_out=None): 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('radionuclide:species:LMMcation'): + if 'LMMcation' in self.name_species: 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 @@ -810,45 +838,40 @@ def update_radionuclide_diameter(self,sp_in=None,sp_out=None): dia_diss=self.get_config('radionuclide:dissolved_diameter') - # Transfer to reversible particles - self.elements.diameter[(sp_out==self.num_prev) & (sp_in!=self.num_prev)] = dia_part - std = self.get_config('radionuclide: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('radionuclide: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))) + # Transform to particle species + to_particles = np.where( (sp_out==self.num_prev) & (sp_in!=self.num_prev) ) [0] + if self.species_slowly_fraction: + to_particles = np.append(to_particles, np.where((sp_out==self.num_psrev) & (sp_in==self.num_ssrev))[0] ) + if self.species_irreversible_fraction: + to_particles = np.append(to_particles, np.where((sp_out==self.num_pirrev) & (sp_in==self.num_sirrev))[0] ) + + new_diam = self.set_init_diameter(self.num_elements_total(), to_particles, dia_part) + self.elements.diameter[to_particles] = new_diam[to_particles] + + + + + # Transform to LMM and colloid species + if 'LMM' in self.name_species: + to_diss = np.where( (sp_out==self.num_lmm) & (sp_in!=self.num_lmm) )[0] + elif 'LMMcation' and 'LMManion' in self.name_species: + to_diss = np.where( (sp_out==self.num_lmmcation) & (sp_in!=self.num_lmmcation) )[0] + to_diss = np.append( to_diss, np.where( (sp_out==self.num_lmmanion) & (sp_in!=self.num_lmmanion) )[0] ) + if 'Colloid' in self.name_species: + to_diss = np.append(to_diss, np.where( (sp_out==self.num_col) & (sp_in!=self.num_col) )[0] ) + if 'Humic_colloid' in self.name_species: + to_diss = np.append( to_diss, np.where( (sp_out==self.num_humcol) & (sp_in!=self.num_humcol) )[0] ) + if 'Polymer' in self.name_species: + to_diss = np.append( to_diss, np.where( (sp_out==self.num_polymer) & (sp_in!=self.num_polymer) )[0] ) + + new_diam = self.set_init_diameter(self.num_elements_total(), to_diss, dia_diss) + self.elements.diameter[to_diss] = new_diam[to_diss] + + - # Transfer to irreversible particles - if self.get_config('radionuclide: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('radionuclide:species:LMM'): - self.elements.diameter[(sp_out==self.num_lmm) & (sp_in!=self.num_lmm)] = dia_diss - if self.get_config('radionuclide:species:LMManion'): - self.elements.diameter[(sp_out==self.num_lmmanion) & (sp_in!=self.num_lmmanion)] = dia_diss - if self.get_config('radionuclide:species:LMMcation'): - self.elements.diameter[(sp_out==self.num_lmmcation) & (sp_in!=self.num_lmmcation)] = dia_diss - # Transfer to colloids - if self.get_config('radionuclide:species:Colloid'): - self.elements.diameter[(sp_out==self.num_col) & (sp_in!=self.num_col)] = dia_diss - if self.get_config('radionuclide:species:Humic_colloid'): - self.elements.diameter[(sp_out==self.num_humcol) & (sp_in!=self.num_humcol)] = dia_diss - if self.get_config('radionuclide:species:Polymer'): - self.elements.diameter[(sp_out==self.num_polymer) & (sp_in!=self.num_polymer)] = dia_diss @@ -856,10 +879,14 @@ def update_radionuclide_diameter(self,sp_in=None,sp_out=None): def bottom_interaction(self,Zmin=None): ''' Change speciation of radionuclides that reach bottom due to settling. particle specie -> sediment specie ''' - if not ((self.get_config('radionuclide:species:Particle_reversible')) & - (self.get_config('radionuclide:species:Sediment_reversible')) or - (self.get_config('radionuclide:slowly_fraction')) or - (self.get_config('radionuclide:irreversible_fraction'))): + + if not ( ('Particle reversible' in self.name_species) & + ('Sediment reversible' in self.name_species) or + self.species_slowly_fraction or self.species_irreversible_fraction): + print('RETURN: ', ('Particle reversible' in self.name_species), + ('Sediment reversible' in self.name_species), + self.species_slowly_fraction, + self.species_irreversible_fraction) return bottom = np.array(np.where(self.elements.z <= Zmin)[0]) @@ -867,12 +894,15 @@ def bottom_interaction(self,Zmin=None): 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('radionuclide:slowly_fraction'): +# self.elements.moving[bottom] = 0 + + if self.species_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('radionuclide:irreversible_fraction'): + + if self.species_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) @@ -886,8 +916,11 @@ def resuspension(self): Sediment species -> Particle specie """ # Exit function if particles and sediments not are present - if not ((self.get_config('radionuclide:species:Particle_reversible')) & - (self.get_config('radionuclide:species:Sediment_reversible'))): + # if not ((self.get_config('radionuclide:species:Particle_reversible')) & + # (self.get_config('radionuclide:species:Sediment_reversible'))): + # return + if not ('Particle reversible' in self.name_species and 'Sediment reversible' in self.name_species): + logger.info('No sediments and particles present...') return specie_in = self.elements.specie.copy() @@ -918,10 +951,12 @@ def resuspension(self): 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('radionuclide:slowly_fraction'): +# if self.get_config('radionuclide:slowly_fraction'): + if self.species_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('radionuclide:irreversible_fraction'): + if self.species_irreversible_fraction: +# if self.get_config('radionuclide: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 @@ -958,7 +993,9 @@ def update(self): logger.info('Speciation: {} {}'.format([sum(self.elements.specie==ii) for ii in range(self.nspecies)],self.name_species)) + #self.elements.moving[self.elements.z <= -1*self.environment.sea_floor_depth_below_sea_level] = 0. + # Horizontal advection self.advect_ocean_current() @@ -973,8 +1010,6 @@ def update(self): - - # ################ # POSTPROCESSING @@ -1407,29 +1442,36 @@ def gui_postproc(self): logger.info ('{:32}: {:>6}'.format(sp,sum(self.elements.specie==isp))) homefolder = expanduser("~") - filename = homefolder+'/conc_radio.nc' + filename = homefolder+'/conc_radio_gui.nc' self.guipp_saveconcfile(filename) + """ zlayer = [-1,-2] self.guipp_plotandsaveconc(filename=filename, outfilename=homefolder+'/radio_plots/RadioConc', zlayers=zlayer, specie= ['Total', 'LMM', 'Particle reversible'], ) + """ - def guipp_saveconcfile(self,filename='none'): + def guipp_saveconcfile(self,filename='none', pixelsize_m=200., reader_sea_depth=None): - for readerName in self.readers: - reader = self.readers[readerName] - if 'sea_floor_depth_below_sea_level' in reader.variables: - break + + 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: + for readerName in self.env.readers: + reader = self.env.readers[readerName] + if 'sea_floor_depth_below_sea_level' in reader.variables: + break self.conc_lon = reader.get_variables('longitude', @@ -1445,8 +1487,8 @@ def guipp_saveconcfile(self,filename='none'): y=[reader.ymin,reader.ymax])['sea_floor_depth_below_sea_level'][:] - self.write_netcdf_radionuclide_density_map(filename, pixelsize_m=200., - zlevels=[-25, -2.], + self.write_netcdf_radionuclide_density_map(filename, pixelsize_m=pixelsize_m, + zlevels=self.depthintervals, activity_unit='Bq', horizontal_smoothing=True, smoothing_cells=1, @@ -1495,14 +1537,14 @@ def guipp_plotandsaveconc(self, filename=None, outfilename=None, zlayers=None, t else: specie_arr = specie - for readerName in self.readers: - reader = self.readers[readerName] - if 'sea_floor_depth_below_sea_level' in reader.variables: - break + # for readerName in self.readers: + # reader = self.readers[readerName] + # if 'sea_floor_depth_below_sea_level' in reader.variables: + # break - gtopo = reader.get_variables('sea_floor_depth_below_sea_level', - x=[reader.xmin,reader.xmax], - y=[reader.ymin,reader.ymax])['sea_floor_depth_below_sea_level'][:] + # gtopo = reader.get_variables('sea_floor_depth_below_sea_level', + # x=[reader.xmin,reader.xmax], + # y=[reader.ymin,reader.ymax])['sea_floor_depth_below_sea_level'][:] From caaf63e8074262ae0019d3e8f9bb4fbbfe4f5752 Mon Sep 17 00:00:00 2001 From: Magne Simonsen Date: Fri, 14 Jun 2024 11:23:05 +0200 Subject: [PATCH 3/6] New check for isotope-specific speciation, improve conc computation --- opendrift/models/radionuclides.py | 78 +++++++++++++++++++++++++++++-- 1 file changed, 73 insertions(+), 5 deletions(-) diff --git a/opendrift/models/radionuclides.py b/opendrift/models/radionuclides.py index e167b61f7..620e1e80b 100644 --- a/opendrift/models/radionuclides.py +++ b/opendrift/models/radionuclides.py @@ -112,10 +112,10 @@ def __init__(self, *args, **kwargs): 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, 'radionuclide:particle_diameter': {'type': 'float', 'default': 5e-6, 'min': .45e-6, 'max': 63.e-6, 'units': 'm', - 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Mean particle diameter. Determines the settling velocity.'}, + 'level': CONFIG_LEVEL_ADVANCED, 'description': 'Mean particle diameter. Determines the settling velocity.'}, 'radionuclide:particle_diameter_uncertainty': {'type': 'float', 'default': 1e-7, 'min': 0, 'max': 100e-6, 'units': 'm', - 'level': CONFIG_LEVEL_ESSENTIAL, 'description': 'Standard deviation of particle size distribution.'}, + 'level': CONFIG_LEVEL_ADVANCED, 'description': 'Standard deviation of particle size distribution.'}, 'radionuclide:particle_diameter_minimum': {'type': 'float', 'default': 0.45e-6, 'min': 0, 'max': 100e-6, 'units': 'm', 'level': CONFIG_LEVEL_ADVANCED, 'description': 'Mimimum particle size.'}, @@ -243,8 +243,11 @@ def init_species(self): self.name_species.append('LMM') self.name_species.append('Particle reversible') self.name_species.append('Sediment reversible') - else: - logger.error('No valid specie setup {}'.format(self.get_config('radionuclide:specie_setup'))) + # else: + # if not self.isotope == 'Al': + # logger.error('No valid specie setup {}'.format(self.get_config('radionuclide:specie_setup'))) + # else: + # pass if '+ slow rev' in self.specie_setup.casefold(): @@ -311,11 +314,39 @@ def set_init_diameter(self, num, idxs,diam): + def check_speciation(self): + isotop = self.get_config('radionuclide:isotope') + specie_setup = self.get_config('radionuclide:specie_setup') + + legal_species = { '137Cs' : ['LMM + Rev', + 'LMM + Rev + Slow rev', + 'LMM + Rev + Irrev', + 'LMM + Rev + Slow rev + Irrev'], + '129I' : ['LMM + Rev', + 'LMM + Rev + Slow rev + Irrev'], + '241Am' : ['LMM + Rev', + 'LMM + Rev + Slow rev', + 'LMM + Rev + Slow rev + Irrev'], + 'Al' : ['LMM + Colloid + Rev'], + } + + + if specie_setup in legal_species[isotop]: + logger.info(f'All good, isotop: {isotop}, species: {specie_setup}') + else: + logger.info(f'Not a legal combination of isotop: {isotop} and speciation: {specie_setup}') + logger.error('Illegal speciation for %s: %s',isotop, specie_setup) + exit() + + + def seed_elements(self, *args, **kwargs): + + self.check_speciation() self.init_species() self.init_transfer_rates() @@ -1023,6 +1054,7 @@ def write_netcdf_radionuclide_density_map(self, filename, pixelsize_m='auto', zl time_avg_conc=False, horizontal_smoothing=False, smoothing_cells=0, + reader_sea_depth=None, ): '''Write netCDF file with map of radionuclide species densities and concentrations''' @@ -1048,6 +1080,42 @@ def write_netcdf_radionuclide_density_map(self, filename, pixelsize_m='auto', zl pixelsize_m = 4000 + + + if reader_sea_depth is not None: + from opendrift.readers import reader_netCDF_CF_generic,reader_ROMS_native + reader = reader_sea_depth + else: + for readerName in self.env.readers: + reader = self.env.readers[readerName] + if 'sea_floor_depth_below_sea_level' in reader.variables: + break + + + try: + self.conc_lat = reader.get_variables('latitude', + x=[reader.xmin,reader.xmax], + y=[reader.ymin,reader.ymax])['latitude'][:] + except: + self.conc_lat = reader.get_variables('grid_latitude_at_cell_center', + x=[reader.xmin,reader.xmax], + y=[reader.ymin,reader.ymax])['grid_latitude_at_cell_center'][:] + try: + self.conc_lon = reader.get_variables('longitude', + x=[reader.xmin,reader.xmax], + y=[reader.ymin,reader.ymax])['longitude'][:] + except: + self.conc_lon = reader.get_variables('grid_longitude_at_cell_center', + x=[reader.xmin,reader.xmax], + y=[reader.ymin,reader.ymax])['grid_longitude_at_cell_center'][:] + + self.conc_topo = reader.get_variables('sea_floor_depth_below_sea_level', + x=[reader.xmin,reader.xmax], + y=[reader.ymin,reader.ymax])['sea_floor_depth_below_sea_level'][:] + + + + if density_proj is None: # add default projection with equal-area property density_proj = pyproj.Proj('+proj=moll +ellps=WGS84 +lon_0=0.0') @@ -1385,7 +1453,7 @@ def get_pixel_mean_depth(self,lons,lats): 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') + h = interpolate.griddata((lon_grd[h_grd.mask==False].flatten(),lat_grd[h_grd.mask==False].flatten()), h_grd[h_grd.mask==False].flatten(), (lons, lats), method='linear') return h From fd5d20493ebfeca0e185e869889fe6b2addae11c Mon Sep 17 00:00:00 2001 From: Magne Simonsen Date: Mon, 17 Jun 2024 10:43:49 +0200 Subject: [PATCH 4/6] switch the depthintervals config to string datatype --- opendrift/models/radionuclides.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/opendrift/models/radionuclides.py b/opendrift/models/radionuclides.py index 620e1e80b..541115aba 100644 --- a/opendrift/models/radionuclides.py +++ b/opendrift/models/radionuclides.py @@ -194,10 +194,9 @@ def __init__(self, *args, **kwargs): 'min': 0, 'max': 1, 'units': 'm/s', 'level': CONFIG_LEVEL_ADVANCED, 'description': ''}, # OUTPUT config - 'radionuclide:output:depthintervals': {'type':'enum', - 'enum': ['-50.', '-25.', '-10.', '-5.', '-2.', '-1.', - '-25., -5.', '-50., -10., -1.', - '-10., -1.'], 'default':'-10.', + 'radionuclide:output:depthintervals': {'type':'str', + 'default': '-25, -10., -5., -1.', + 'min_length':0, 'max_length':60, 'level':CONFIG_LEVEL_ESSENTIAL, 'description':'Depth intervals for computation of concentration'}, }) From 8061083c9c4b6cfa3fd859d7f505501e72b68271 Mon Sep 17 00:00:00 2001 From: Magne Simonsen Date: Mon, 17 Jun 2024 13:29:25 +0200 Subject: [PATCH 5/6] Apply new isotope config to the test case --- tests/radionuclide/test_model.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/tests/radionuclide/test_model.py b/tests/radionuclide/test_model.py index 9d363eb5a..ac34f2e69 100644 --- a/tests/radionuclide/test_model.py +++ b/tests/radionuclide/test_model.py @@ -18,11 +18,9 @@ def test_radio_nuclide(test_data): o.set_config('vertical_mixing:timestep', 600.) # seconds o.set_config('drift:horizontal_diffusivity', 10) o.set_config('radionuclide:particle_diameter',5.e-6) # m - o.set_config('radionuclide:transformations:Kd',2.e0) # (m3/kg) - o.set_config('radionuclide:sediment:resuspension_depth',2.) - o.set_config('radionuclide:sediment:resuspension_depth_uncert',0.1) - o.set_config('radionuclide:sediment:resuspension_critvel',0.15) - o.set_config('radionuclide:transfer_setup','Bokna_137Cs') + o.set_config('radionuclide:isotope','137Cs') + o.set_config('radionuclide:specie_setup','LMM + Rev') + o.set_config('general:coastline_action', 'previous') o.set_config('seed:LMM_fraction',.45) o.set_config('seed:particle_fraction',.55) From 9ddab2ebdcc2def5336861f41b87e599fa4ff73c Mon Sep 17 00:00:00 2001 From: Magne Simonsen Date: Mon, 17 Jun 2024 13:39:12 +0200 Subject: [PATCH 6/6] update example --- examples/example_radionuclides.py | 54 ++++--------------------------- 1 file changed, 7 insertions(+), 47 deletions(-) diff --git a/examples/example_radionuclides.py b/examples/example_radionuclides.py index 2432f679a..652bc86e7 100755 --- a/examples/example_radionuclides.py +++ b/examples/example_radionuclides.py @@ -22,33 +22,14 @@ # Adjusting some configuration o.set_config('drift:vertical_mixing', True) -#o.set_config('environment:constant:ocean_vertical_diffusivity', 0) #o.set_config('vertical_mixing:diffusivitymodel','constant') # include settling without vertical turbulent mixing -o.set_config('vertical_mixing:diffusivitymodel','environment') # include settling without vertical turbulent mixing +o.set_config('vertical_mixing:diffusivitymodel','environment') # apply vertical diffusivity from ocean model # Vertical mixing requires fast time step o.set_config('vertical_mixing:timestep', 600.) # seconds o.set_config('drift:horizontal_diffusivity', 10) #%% -# Activate the desired species -#o.set_config('radionuclide:species:LMM', True) -#o.set_config('radionuclide:species:LMMcation', True) -#o.set_config('radionuclide:species:LMManion', True) -#o.set_config('radionuclide:species:Colloid', True) -#o.set_config('radionuclide:species:Humic_colloid', True) -#o.set_config('radionuclide:species:Polymer', True) -#o.set_config('radionuclide:species:Particle_reversible', True) -#o.set_config('radionuclide:species:Particle_irreversible', True) -#o.set_config('radionuclide:species:Particle_slowly_reversible', True) - -#o.set_config('radionuclide:species:Sediment_reversible', True) -#o.set_config('radionuclide:species:Sediment_irreversible', True) -#o.set_config('radionuclide:species:Sediment_slowly_reversible', True) - o.set_config('radionuclide:particle_diameter',5.e-6) # m -#o.set_config('radionuclide:particle_diameter_uncertainty',1.e-7) # m -o.set_config('radionuclide:transformations:Kd',2.e0) # (m3/kg) -#o.set_config('radionuclide:transformations:slow_coeff',1.e-6) o.set_config('radionuclide:sediment:resuspension_depth',2.) o.set_config('radionuclide:sediment:resuspension_depth_uncert',0.1) @@ -56,16 +37,12 @@ # -#o.set_config('radionuclide:transfer_setup','custom') -o.set_config('radionuclide:transfer_setup','Bokna_137Cs') -#o.set_config('radionuclide:transfer_setup','137Cs_rev') -#o.set_config('radionuclide:transfer_setup','Sandnesfj_Al') +o.set_config('radionuclide:isotope', '137Cs') +o.set_config('radionuclide:specie_setup','LMM + Rev') # By default, radionuclides do not strand towards coastline o.set_config('general:coastline_action', 'previous') o.set_config('general:seafloor_action','lift_to_seafloor') -#o.set_config('general:seafloor_action','previous') -#o.set_config('general:use_auto_landmask',False) o.set_config('seed:LMM_fraction',.45) @@ -80,20 +57,13 @@ td=datetime.today() time = datetime(td.year, td.month, td.day, 0) -#latseed= 61.2; lonseed= 4.3 # Sognesjen -#latseed= 59.0; lonseed= 10.75 # Hvaler/Koster -#latseed= 57.5; lonseed= 9.3 # Kattegat -latseed= 60.0; lonseed= 4.5 # Bergen (?) +latseed= 60.0; lonseed= 4.5 ntraj=5000 iniz=np.random.rand(ntraj) * -10. # seeding the radionuclides in the upper 10m o.seed_elements(lonseed, latseed, z=iniz, radius=1000,number=ntraj, time=time, -# LMM_fraction=0.25, -# particle_fraction=0.75 -# diameter=99.e-6#diam, - #specie=init_speciation ) @@ -137,17 +107,7 @@ -# Postprocessing: write to concentration netcdf file -o.conc_lat = reader_norkyst.get_variables('latitude', - x=[reader_norkyst.xmin,reader_norkyst.xmax], - y=[reader_norkyst.ymin,reader_norkyst.ymax])['latitude'][:] -o.conc_lon = reader_norkyst.get_variables('longitude', - x=[reader_norkyst.xmin,reader_norkyst.xmax], - y=[reader_norkyst.ymin,reader_norkyst.ymax])['longitude'][:] -o.conc_topo = reader_norkyst.get_variables('sea_floor_depth_below_sea_level', - x=[reader_norkyst.xmin,reader_norkyst.xmax], - y=[reader_norkyst.ymin,reader_norkyst.ymax])['sea_floor_depth_below_sea_level'][:] -#o.conc_mask = reader_norkyst.land_binary_mask +# # Postprocessing: write to concentration netcdf file #%% # @@ -160,8 +120,8 @@ smoothing_cells=1, time_avg_conc=True, deltat=2., # hours - llcrnrlon=4.4, llcrnrlat=59.9, - urcrnrlon=4.8, urcrnrlat=60.2, +# llcrnrlon=4.4, llcrnrlat=59.9, +# urcrnrlon=4.8, urcrnrlat=60.2, )