Source code for opendrift.models.sealice

# 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
# 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 <>.
# Copyright 2021, Julien Moreau, Plastic@Bay CIC

from datetime import timedelta, datetime
import numpy as np
import logging; logger = logging.getLogger(__name__)
from opendrift.models.oceandrift import Lagrangian3DArray, OceanDrift
from opendrift.models.physics_methods import hour_angle

[docs] class SeaLiceElement(Lagrangian3DArray): """ Extending Lagrangian3DArray with specific properties for larval and juvenile stages of sea lice into super individuals """ variables = Lagrangian3DArray.add_variables([ ('LicePerFish',{'dtype':np.float32, 'units':'', 'default':0.5}), ('AvFishW8',{'dtype':np.float32, 'units':'kg', 'default':4.5}), ('particle_biomass',{'dtype': np.float32, 'units': 'kg', 'default': 1000.}), ('hatched',{'dtype': np.float32, 'units': '', 'default': 0.}), ('nauplii', {'dtype': np.float32, 'units': '', 'default': 0.}), ('copepodid', {'dtype': np.float32, 'units': '', 'default': 0.}), ('dead', {'dtype': np.float32, 'units': '', 'default': 0.}), ('eliminated',{'dtype': np.int32, 'units': '', 'default': 0}), ('degree_days', {'dtype': np.float32, 'units': '', 'default': 0}), #range 40-170 ('safe_salinity_above', {'dtype': np.int32, 'units': '', 'default': 0}), ('temperature_above', {'dtype': np.float32, 'units': '', 'default': 10}), ('temperature_below', {'dtype': np.float32, 'units': '', 'default': 10}), ('light', {'dtype': np.float32, 'units': 'µmol photon s−1 m−2', 'default': 0.}), ])
[docs] class SeaLice(OceanDrift): """ Particle trajectory model based on the OpenDrift framework. Developed by Julien Moreau (Plastics@Bay CIC) Generic module for particles that are subject to vertical turbulent mixing with the possibility for positive or negative buoyancy Particles are sea-lice (Lepeophtheirus salmonis). """ ElementType = SeaLiceElement required_variables = { 'x_sea_water_velocity': {'fallback': 0}, 'y_sea_water_velocity': {'fallback': 0}, # 'sea_surface_wave_significant_height': {'fallback': 0}, # 'x_wind': {'fallback': 0}, # 'y_wind': {'fallback': 0}, 'land_binary_mask': {'fallback': None}, 'sea_floor_depth_below_sea_level': {'fallback': 50}, 'surface_net_downward_radiative_flux':{'fallback': 0}, 'ocean_vertical_diffusivity': {'fallback': 0.01}, 'sea_water_temperature': {'fallback': 10}, 'sea_water_salinity': {'fallback': 34} } # required_profiles_z_range = [0, -50] # The depth range (in m) which profiles should cover def __init__(self,*args, **kwargs): # Calling general constructor of parent class # super(SeaLice, self).__init__(*args, **kwargs) super().__init__(*args, **kwargs) # Configuration options self._add_config({ 'general:duration':{'type':'float', 'default':0., 'min': None, 'max': None, 'units': 'seconds', 'description': 'Experiment time in seconds', 'level': CONFIG_LEVEL_ESSENTIAL}, 'lice:seeding_time_step':{'type':'float', 'default':None, 'min': None, 'max': None, 'units': 'seconds', 'description': 'Time between particle release', 'level': CONFIG_LEVEL_ESSENTIAL}, 'lice:death_rate':{'type':'float', 'default':0.01/3600, 'min': 0., 'max': None, 'units': 's-1', 'description': 'Rate of Larvae death per seconds', 'level': CONFIG_LEVEL_BASIC}, 'lice:maturation_rate':{'type':'float', 'default':0.1/3600, 'min': 0., 'max': None, 'units': 's-1', 'description': 'Rate of Nauplii maturation in Copepodids', 'level': CONFIG_LEVEL_BASIC}, 'lice:maturity_date':{'type':'float', 'default':3.63, 'min': 0., 'max': None, 'units': 'days', 'description': 'Days to start maturing into Copepodids', 'level': CONFIG_LEVEL_BASIC}, 'lice:sinking_velocity':{'type':'float', 'default':0.00025, 'min': 0., 'max': 0.01, 'units': 'm.s-1', 'description': 'Larvae sinking velocity', 'level': CONFIG_LEVEL_BASIC}, 'lice:vertical_migration_speed':{'type':'float', 'default':0.00075, 'min': 0., 'max': 0.01, 'units': 'm.s-1', 'description': 'Larvae vertical speed', 'level': CONFIG_LEVEL_BASIC}, 'lice:freezing_salinity':{'type':'float', 'default':27., 'min': 0., 'max': 35., 'units': 'PSU', 'description': 'Salinity immobilising larvae', 'level': CONFIG_LEVEL_BASIC}, 'lice:avoided_salinity':{'type':'float', 'default':32., 'min': 0., 'max': 50., 'units': 'PSU', 'description': 'Salinity actively avoided', 'level': CONFIG_LEVEL_BASIC}, 'lice:nu':{'type':'float', 'default':500., 'min': 100, 'max': 1000, 'units': 'nm', 'description': 'Wavelength used to calculate irradiance', 'level': CONFIG_LEVEL_ADVANCED}, 'lice:k_water':{'type':'float', 'default':0.2, 'min': -10., 'max': 0., 'units': '', 'description': 'coefficient of exponential decay of light in water', 'level': CONFIG_LEVEL_ADVANCED}, 'lice:Nauplii_light_trigger':{'type':'float', 'default':2.E-5, 'min': 0., 'max': 1., 'units': 'µmol photon s−1 m−2', 'description': 'light detection threshold of Nauplii', 'level': CONFIG_LEVEL_BASIC}, 'lice:Copepodid_light_trigger':{'type':'float', 'default':0.392, 'min': 0.0, 'max': 1.0, 'units': 'µmol photon s−1 m−2', 'description': 'light detection threshold of copepodids', 'level': CONFIG_LEVEL_BASIC}, 'lice:twilight':{'type':'float', 'default':15, 'min': 0., 'max': 90., 'units': 'degrees', 'description': 'angle below the horizon for twilight', 'level': CONFIG_LEVEL_ADVANCED}, })
[docs] def prepare_run(self):"preparing the run...") ### need to find a better way of initialising the variables... self.prefix="lice:" self.vertical_migration_speed=self.get_config(self.prefix+'vertical_migration_speed') \ *self.time_step.total_seconds() self.sensing_distance=2*self.vertical_migration_speed self.tau= np.pi+2*np.deg2rad(self.get_config(self.prefix+'twilight')) #width of solar irridation distribution'nu') self.sinking_velocity=self.get_config(self.prefix+'sinking_velocity') \ *self.time_step.total_seconds() self.freezing_salinity=self.get_config(self.prefix+'freezing_salinity') self.avoided_salinity=self.get_config(self.prefix+'avoided_salinity') self.k_water=self.get_config(self.prefix+'k_water') self.Nauplii_light_trigger=self.get_config(self.prefix+'Nauplii_light_trigger') self.Copepodid_light_trigger=self.get_config(self.prefix+'Copepodid_light_trigger') self.new_born() self.population() self.ref_date = datetime(1,1,1,0,0)
[docs] def new_born(self): """ Approach by Rittenhouse et al., (2016) [Rittenhouse, M.A., C.W. Revie and A. Hurford, 2016. A model for sea lice (Lepeophtheirus salmonis) dynamics in a seasonally changing environment. Energetics, 16, 8–16., whereby the numbers of lice nauplii released every day (NP) are estimated according to the following equation: NP(t)=ηενA(t) where η = 592 is the number of eggs per clutch (Heuch et al., 2000; Rittenhouse et al., 2016), ε = 0.0476 is the egg string production rate per day (Heuch et al., 2000; Rittenhouse et al., 2016), ν = 0.6 is the hatching success i.e. the proportion of eggs which produce viable nauplii (Johnson and Albright, 1991), and A(t) is the total number of adult female lice on each farm, derived here from lice target levels and number of fish on site. """ logger.debug("Calculating standard spawn of nauplii") Ne, nu = 592 , 0.6 Eps=0.0476/timedelta(days=1).total_seconds()*self.get_config( \ self.prefix+'seeding_time_step') self.spawn=Ne* Eps* nu
[docs] def population(self): """ Simulate the population evolution according to the biological parameters Parameters ---------- t: numpy array (float) Time abscissa in seconds through the experiment self.Mat: int Age in timestep when the nauplii start to become adult copepodids maturation_rate: float Probability of maturation by time-step (t>= Mat) death_rate: float Probability of death by time-step Number of lice per particle when generated is standardised to 1 Returns ------- juv, adults, dead: np.arrays(int) Juveniles(Nauplii), adults(Copepodid), dead lices arrays through the duration of the experiment """ logger.debug("Building global population model") death_rate=self.get_config(self.prefix+'death_rate')* self.time_step.total_seconds() maturation_rate=self.get_config(self.prefix+'maturation_rate')* self.time_step.total_seconds() duration = self.get_config('general:duration')/ self.time_step.total_seconds() Mat = int(np.ceil(self.get_config(self.prefix+'maturity_date')*24*3600/ \ self.time_step.total_seconds())) # maturity age in timestep t=np.arange(0,duration+1,dtype=np.int32) self.juv=np.exp(-1*death_rate*t) self.juv[0]=1 self.dead=1-self.juv, dtype=np.float32) if Mat<duration: decayjuv=self.juv[Mat]*np.exp(-1*(maturation_rate+death_rate)*(t[t>=Mat]-Mat)) self.juv[t>=Mat] =decayjuv[t>=Mat]=1-self.dead[t>=Mat]-self.juv[t>=Mat]
[docs] def SI_pop(self): """ distribute the age fractions in the particles """ logger.debug("Aging the super_individuals") scaling = self.elements.particle_biomass*self.elements.LicePerFish/ \ self.elements.AvFishW8 #identify new spawns New_release = self.elements.age_seconds<=self.time_step.total_seconds() self.elements.hatched[New_release]=self.spawn*scaling[New_release] self.elements.nauplii[New_release]=self.elements.hatched[New_release] Free = ~New_release time_in_step=(self.elements.age_seconds[Free]/ \ self.time_step.total_seconds()).astype(np.int32) self.elements.nauplii[Free]=self.elements.hatched[Free]* \ self.juv[time_in_step] self.elements.copepodid[Free]=self.elements.hatched[Free]* \[time_in_step] self.elements.dead[Free]=self.elements.hatched[Free]* \ self.dead[time_in_step] # desactivate dead particles ## this test is not good for seeds with few lice Dying=self.elements.nauplii+ self.elements.copepodid<1 self.elements.eliminated[Dying]=1 self.deactivate_elements(self.elements.eliminated.astype(np.bool), reason="All dead")
[docs] def sensing(self): """ Lice sensing if above or bellow the conditions are better for them within a specified distance. moving the particles up and then down then back in their initial position """ logger.debug("sensing temperature and salinity") Backup= np.copy(self.elements.z[:]) self.elements.z +=self.sensing_distance self.elements.safe_salinity_above[self.environment.sea_water_salinity> \ self.avoided_salinity]=1 self.elements.temperature_above= self.environment.sea_water_temperature self.elements.z -= 2*self.sensing_distance self.elements.temperature_below= self.environment.sea_water_temperature self.elements.z = Backup
[docs] def degree_days(self): """ Calculate the degree days of a particles XXX: under development """ logger.debug("accumulate degree_days **Experimental**") ### define active elements self.elements.degree_days+=self.environment.sea_water_temperature* \ self.time_step.total_seconds()/timedelta(days=1).total_seconds()
# def solar_noon(self): # """ # Search solar noon at the longitude # """ #,minute=0,second=0) # Av_lon= np.mean(self.elements.lon) # Av_lat= np.mean( # minutes=np.arange(0,60*24, dtype=np.float32) # angles = np.empty_like(minutes) # for i in range(len(minutes)): # angles[i]=hour_angle(self.ref_date+timedelta(minutes=float(minutes[i])),Av_lon) # # solar_noon=o.ref_date+timedelta(minutes=float(minutes[np.argmin(np.abs(angles))])) # self.noon_elevation=solar_elevation(solar_noon, Av_lon, Av_lat)
[docs] def irradiance(self): """ Distribute the daily energy from irradiance with a gaussian distribution. We use the twilight times for high sensitivity organisms Convert irradiance from W.m-2 to micromol photon.s-1.m-2 calculate the photon flux in the water according to exponential decay from the sea surface """ self.elements.solar_angle=np.deg2rad(hour_angle(self.time,np.mean(self.elements.lon))) solar_coeff= np.sqrt(self.tau/(2*np.pi))*np.exp(-self.tau/2*self.elements.solar_angle**2) self.elements.light = solar_coeff* \ self.environment.surface_net_downward_radiative_flux* \*0.00836* \ np.exp(self.k_water*self.elements.z)
[docs] def depth_test(self): """ The natural range of the larvae is 0-50m """ self.elements.z[self.elements.z>0]=0 self.elements.z[self.elements.z<-50]=-50
[docs] def Lice_vertical_migration(self): """ Make larval sea lice to migrate vertically according to salinity, light and temperature triggers. """ ### all lice sink at the same speed self.elements.z -= self.sinking_velocity ### filter actively avoiding salt # Frozen=self.environment.sea_water_salinity<self.freezing_salinity Avoiding = np.logical_and((self.freezing_salinity< \ self.environment.sea_water_salinity), (self.environment.sea_water_salinity<self.avoided_salinity)) Normal_salt = self.environment.sea_water_salinity>=self.avoided_salinity self.sensing() ### Irradiance computation self.irradiance() ### identify the elements involved in the different scenarios Filter_N= self.elements.copepodid < self.elements.nauplii Filter_C=~Filter_N safe_up_salt=Normal_salt&self.elements.safe_salinity_above.astype(np.bool) #### need to be able to deal with false arrays of Filter_N and C... ### They generate empty arrays light_mig_N= safe_up_salt&Filter_N&(self.elements.light \ > self.Nauplii_light_trigger) light_mig_C= safe_up_salt&Filter_C&(self.elements.light \ > self.Copepodid_light_trigger) up_temp_mig= safe_up_salt&~light_mig_N&~light_mig_C& \ (self.elements.temperature_above \ > self.environment.sea_water_temperature) down_temp_mig=Normal_salt&~up_temp_mig&~light_mig_N& ~light_mig_C& \ (self.elements.temperature_below> \ self.environment.sea_water_temperature) ### Active migration going_down=np.logical_or(Avoiding,down_temp_mig) going_up=np.logical_or(np.logical_or(light_mig_N,light_mig_C),up_temp_mig) logger.debug("{} going down, {} going up".format(going_down.sum(), going_up.sum())) self.elements.z[going_down] -= self.vertical_migration_speed self.elements.z[going_up] += self.vertical_migration_speed
[docs] def update(self): self.SI_pop() self.degree_days() self.advect_ocean_current() self.vertical_mixing() self.Lice_vertical_migration() self.depth_test()