Modèle de la variation saisonnière de la concentration de Diclofenac dans un étang de stabilisation d’eaux usées en Arctique
Application à une communauté arctique
Jérémie Boucher-Fontaine
06 December 2020
Contexte
La gestion des eaux usées est particulièrement compliquée en Arctique du fait du sol rocheux ou du pergélisol qui empêche la mise en place de systèmes communs dans les régions du sud. Le gel, les inondations printanière, les quantité limité d’eau, des coût important de l’électricité, de l’essence et du transport, et l’accessibilité limité à ces communautés contribue à en faire un défi de taille qui n’est toujours pas réglé à l’heure actuelle (Gunnarsdóttir et al. 2013; Johnson et Inuit Tapiriit Kanatami 2008). Par conséquent, nombre de ces systèmes ne répondent pas aux normes appliquées dans les autres régions du Canada (Daley et al. 2015).
Malgré l’amélioration constante de la situation, et malgré la mise en place de plan de gestions considérant les réalités des populations Inuits (Johnson et Inuit Tapiriit Kanatami 2008), cet enjeu demeure préoccupant pour ces communautés (Daley et al. 2015). Parmi les préoccupations figure la contamination possible des animaux sauvages, qui constituent une part importante de l’alimentation de ces communautés. Il est ainsi pertinent de se pencher sur les conséquences des eaux usées sur les écosystèmes côtiers dans lesquelles se déversent éventuellement les eaux plus ou moins traités (Daley et al. 2015).
Selon la pente du terrain et la nature du sol, des systèmes de canalisation ont été mis en place dans différentes communautés. À Iqaluit, la capitale du Nunavut, est la seule ville ayant un système de canalisation acheminant directement une partie des eaux usées au site de traitement, combiné avec un système de traitement par camion. Dans la majorité des communautés, les systèmes de canalisation des habitations aboutissent dans un réservoir d’eaux usée, régulièrement vidée par des camions. Ces camions acheminent les eaux usées à une décharges ou une zone de traitement d’eaux usées désignées, soit généralement des étangs de stabilisation des déchets. (Daley et al. 2015; Gunnarsdóttir et al. 2013; Johnson et Inuit Tapiriit Kanatami 2008).
Dans les étangs de stabilisation, les eaux usées gèlent en hiver, et au printemps se déversent dans les zones humides avoisinantes et finissent dans la mer (Daley et al. 2015). Ainsi, toutes les eaux usées accumulées durant l’hiver se mobilisent et se déversent en peu de temps dans l’écosystème environnant. Ce système entre dans la catégorie des systèmes de traitements passifs.Malheureusement, ce type de système comporte des failles importantes.
Les eaux usées issues des zones humides atteignant la mer contiennent de nombreuses substance potentiellement nocives pour l’environnement et la santé humaine (Gunnarsdóttir et al. 2013). Une grande partie de ces composés sont des produits pharmaceutiques et de soins personnels (PPCS), dont les eaux usée des villes constituent la plus grande source (Gunnarsdóttir et al. 2013; Kallenborn et al. 2018). De faits, les composés pharmaceutiques administrés aux humains et aux animaux sont en grande partie excrété sans avoir été transformés. Ils se retrouvent ainsi dans les eaux usées des maisons et des hôpitaux. Biens qu’ils soient généralement éventuellement dégradé, le peu de traitements des eaux usées en Arctique contribue à rejeter des quantité notable de ces produits et des produits de dégradation dans l’environnement (Gunnarsdóttir et al. 2013; Kallenborn et al. 2018). Entre autres, des antibiotiques, des agents antimicrobiens, des stimulants et des antidépresseurs sont retrouvés jusque dans les écosystèmes marins (Kallenborn et al. 2018).
Ce phénomène est exacerbé en Arctique du fait des conditions environnementale. La température de l’eau étant basse toute l’année, le taux de dégradation de ces composés pharmaceutique est faible. Dans un même ordre d’idée, leur dégradation par des phénomène chimique et microbiologique est atténuée par la faible température et l’absence de lumière et de processus photosynthétique en hiver (Gunnarsdóttir et al. 2013; Kallenborn et al. 2018). Cela mène à des concentrations de composés supérieures à ce que l’on retrouve à proximité de ville du Sud, ce qui indique que l’émission incontrôlée de ces composés, associée à des conditions environnementales favorisant la pérennité de ces composés dans l’environnement prédomine sur la densité des populations humaines (Kallenborn et al. 2018). Ainsi, malgré la faible population en Arctique, la libération de PPCP et l’accumulation d’en l’environnement est importante.
Or, ces composés peuvent avoir des conséquences variées sur l’environnement. Les effets des drogues, stimulants et médicaments aboutissant dans l’océan demeure en grande partie méconnue et probablement très variés (Gunnarsdóttir et al. 2013). À titre d’exemple, certains composés pharmaceutique ont été associés à des perturbations endocrines chez des poissons, nuisant par exemple à la production d’œufs par des poissons (Rodgers-Gray et al. 2000). Il a aussi été démontré que l’exposition à des produits chimiques variés retrouvés dans l’environnement pouvait moduler la relation entre hôte et microbiome (Chiu et al. 2020; Jin et al. 2017). De fait, une exposition à des polluants chimiques dans l’eau modifierait la composition du microbiome et pourrait ainsi contribuer à la toxicité de ces composés (Claus, Guillou, et Ellero-Simatos 2016) constituant une menace pour ces organismes (Jin et al. 2017).
Ainsi, l’objectif de ce projet est de modéliser la concentration d’un de ces composés pharmaceutique courant en fonction du temps au sein de étang de stabilisation et par conséquent dans les effluents atteingnant les écosystèmes côtier. Le composés choisi est le Diclofenac, aussi connu sous le nom de Voltaren au Canada. Cette modélisation de la variation saisonnière de la concentration de ce composé permettera notamment d’évaluer l’importance relative de la population par rapport aux conditions environnementales quant à l’accumulation de ce produit dans les effluents atteingnants écosystèmes côtiers à proximité des villes, permettant de vérifier l’importance de considérer ce phénomène en Arctique.
Modèle conceptuel
Pour ce modèle, on considère donc 3 variables d’état, et 6 flux
Variables d’états
DFC : Concentration de Diclofenac dans le bassin de rétention
TPphoto : Concentration des produits de phototransformation, principalement du 8-chlorocarbazole-1-acetic acid
TPbio : Concentration des produits de biotransformation, principalement du 5-hydroxydiclofenac
Flux
f0 = apport de DFC par les eaux usées ()
f1 = Photodégradation du DFC
f2 = Biodégradation du DFC
f3 = photodégradation secondaire : dégradation des produits de phototransformation (8-chlorocarbazole-1-acetic acid)
f4 = Adsorption par sédiments : dégradation des produits de biotransformation (5-hydroxydiclofenac)
f5 = Dilution des produits de photodégradation par les eaux usées
f6 = Dilution des produits de biodégradation par les eaux usées
Taux de changement
\(\frac{dDFC}{dt}\) = f0 - f1 - f2 \(\frac{dTPphoto}{dt}\) = f1 - f3 + f5 \(\frac{dTPbio}{dt}\) = f2 - f4 + f6
Formulation mathématique
f0, f5 et f6 : Arrivée d’eaux usées
L’apport total de composés par les eaux usées dépends de trois choses, soit de la cocentration moyenne de ces composés dans les eaux usées non traités (Cinfluent), du volume d’eaux usées généré par personne (Vinfluent) et du nombre de personne (POP).
Afin de le traduire en augmentation de la concentration de DFC par les bassins de rétention, il faut considérer le volume du bassin de rétention (Vbassin)
On peut décrire le changement de concentration de ces composé dans le bassin de rétention par l’équation C1V1=C2V2,
où :
C1 = variation de concentration = Cinfluent - DFC
V1 = Vinfluent * POP
C2 = f0 (Augmentation de Concentration de DFC dans le bassin de rétention)
V2 = Vbassin + Vinfluent * pop * t
Ainsi,
f0 = \(\frac{(Cinfluent - DFC) * Vinfluent * POP}{Vbassin + Vinfluent*POP*t}\)
f0 = (Cinfluent* Vinfluent * POP) / (Vbassin +VinfluentPOPt)
Pour F5 et F6, on considère qu’il n’y a pas de produits de dégradation dans les eaux usées et donc que Cinfluent = 0. Ainsi,:
f5 = \(\frac{-TPphoto * Vinfluent * POP}{Vbassin + Vinfluent*POP*t}\)
f6 = \(\frac{-TPbio * Vinfluent * POP}{Vbassin + Vinfluent*POP*t}\) où :
Paramètres :
Cinfluent : la cocentration moyenne de DFC dans les eaux usées non traités = 340 ng/L
Vinfluent : volume d’eaux usées généré par personne par heure = 140 L/p*d
POP : nombre de personne dans la communautés = 1000 p
Vbassin : Volume du bassin = 75 000 000 L = superficie en ha * 10 000 * profondeur du bassin
Tous ces éléments sont pour l’instant considérés comme des paramètres. Il aurait été possible de faire de Vinfluent un forçage si l’on considère le dégèle au printemps et les variation de débit. POP pourrait aussi devenir un forçage si on considère une variation de population pendant l’année
Afin de pourvoir comparé appliquer ce modéle à des communautés de tailles différente, il serait nécessaire de faire de Vbassin un paramètre dépendant la taille de la population, puisque les villes ont souvent un plus grand bassin de stabilisation. Toutefois, afin d’étudier l’importance des différents paramètres sur ce modèle basé sur un bassin en particulier, il est plus utile de considérer la taille du bassin comme fixe, puisque dans la réalité il ne s’agrandit pas avec le nombre de personne. Voir page 441 de Smith 1996 pour équation linéaire de taille de bassin en fonction de la population et des émissions. Cependant, selon Ragush et al (2015), la taille des bassins en place en Arctique sont plusieurs fois plus petite que la taille surggérée. J’utiliserai ainsi une taille de bassin théorique idéalisée.
Pour le volume du bassin, étant donné que la plupart des bassins de rétentions sont de type facultatif, on prend une profondeur moyenne de 1.5 m. Selon Smith 1996, pour une population de 1000 personnes et un temps de rétention de 12 mois, les bassins couvrent environ 5 ha. Cela fait 5 * 10 000 m^2 x 1,5 m de profondeur, soit 75 000 m^3, soit 75 000 000 Litre
Le volume totale est pour sa part calculé en additionnant au volume du bassin le volume d’eau arrivée depuis t(0), ce qui donne Vtotal <- Vbassin + VinfluentPOPt
Cela mène à un volume toujours croissant. Ceci est voulu étant donnnée que les eaux usées des bassins de l’arctique sont généralement vidée une fois par année durant l’été, en générale début septembre. Cela permet donc déterminer la concentration dans des effluents qui atteingnent les écosystèmes, soit la variable qui nous intéresse véritablement.
f1 : Photodégradation
La photodégradation dépend de la concentration de DFC, d’un coefficient de photodégradation lui-même déterminé par l’éclairement incident, et limité par la pénétration de la lumière et par la température.
f1 = kphoto * [DFC] * limitation de la température * limitation par pénétration de la lumière
Détermination du coefficient de photodégradation
Le taux de photodégradation selon la lumière est déterminé par le rendement quantique (Quantum Yield), qui détermine le taux de dégradation par mol de photon (en Einstein). Selon Andreozzi, = 0.0375 mol Einstein-1 par la lumière du soleil à un PH de 5.5.
Le coefficient de photodégradation de basse est donnée par l’équation suivante (1- exp(-QYI04.57 * 3600)). Le I0 qui représente l’éclairement incident en W m-2. Pour convertir l’éclairement incident en photons afin de l’utiliser par le rendement quantique, il faut, selon http://www.egc.com/useful_info_lighting.php, multiplier l’élcairement par 4,57 pour obtenir des uEinstein, puis par 3600 pour trouver l’éclairement par heure.
Paramètres :
QY : Quantum Yield, ou rendement quantique, en µMol/µEinstein.
Forçage :
I0: est l’incidence lumineuse, en Watt/ m -2, convertie en µEinstein
Il faut également considérer que lorsqu’il y a une couvert de glace lorsque le bassin est gelé, et qu’alors il n’y a plus de lumière. On peut alors considérer que I0 est = 0. Pour ce faire, j’utiliserai la température de l’eau de la mer déterminée dans la baie de Frobisher, à proximité d’Iqaluit.
Donc : I0 = 0 lorsque temp < 273 kelvin
limitation de la température
L’effet de la température se fait au niveau du coefficient. Lorsque la température augemente, le coefficient de vitesse de la réaction augmente. Le coefficient se trouve donc en fonction de la température, en suivant la loi d’Arrhenius, selon des paramètres de E et de A trouvés expérimentalement. Afin de combiner cette limitation à celle de la lumière, le A est remplacé par le coefficient de vitesse trouver en fonction de la lumière, en h-1.
Selon la loi d’Arrhenius:
- e^\(\frac{-E}{RTemp}\))
Où
Paramètres :
A :facteur de fréquence, remplacé par la constante de vitesse de photodégradation déterminé en fonction de la lumière, en heure -1
E : énergie d’activation d’Arrhenius donnée en J·mol-1 = 5302 J mol−1
R : la constante des gaz parfait = 8.314462 J mol−1 K−1
Forcings :
Temp est la température de l’eau en Kelvin
ce qui donne une équation de flux de
Limitation par la pénétration de la lumière :
L’effet de la profondeur pourrait être calculé selon un coeffficent d’atténuation et une intrégrale selon la profondeur.
Toutefois, une étude de Nurmi (2019) suggère de considérer un taux de 100 % pour une profondeur de 0 - 0.1 m, un taux de 10 % de 0.1 à 1m, et un taux de 1 % entre 1 et 5 m. La profondeur standards des bassins étant d’environ de 1.5 m, et considérant que la couche d’eau est entièrement mélangé, il est possible de faire une moyenne pondérée de cette atténuation par la pénétration de la lumière, soit (1 * 0.10 + 0.1 * 0.90 + 0.01 * 0.5 ) / 1.5 = 0.13.
Pour pouvoir faire de la profondeur un paramètre testable, dans on peut utiliser dans l’équation le terme (1 * 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas)), en considérant que la profondeur du bassin est d’au moins 1 mètres, comme le suggère la littérature. Je suis conscient que la profondeur s’accentue probablement avec l’augmentation du niveau total de l’eau, mais je toruvais qu’il n’était pas nécessaire d’ajouter cette modulation au modèle étant donnée une complexité inutile.
On peu donc considérer que la photodégradation n’a lieux qu’à 13 % de son taux maximale en moyenne L’équation devient alors :
f1 = [DFC] * (1- e^\(\-QY*I0*4.57\)) - e^\(\frac{-E}{RTemp}\)) * (1 * 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas))
f1 = DFC * (1- exp(-QYI04.57)) * exp((- E/(R * Temp))) (1 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas))
f2 : Biodégradation
La biodégradation ne représente qu’une faible partie de la dégradation du DFC. Elle a lieux dans les sédiements, et elle dépend surtout de la température.
Je n’ai pas trouvé d’étude décrivant l’effet de la température en terme assez claire pour avoir des valeurs précises de constante d’Arrhenius. J’utiliserai donc un terme de limitation de q10 et une valeur de q10 standard de 2, qui pourrait idéalement être préciser par la suite.
f2 = Kbio * [DFC] * e^(\(\frac{Temp-Tref}{10}\)*ln(q10))
f2 = Kbio * DFC * exp(((Temp-Tref)/10)*log(q10)
où :
Paramètres :
q10 = Facteur d’augmentation pour chaque augmentation de 10°C = 2 (valeur standard à modifier éventuellement)
Kbio= constante de vitesse de biodégradation = Entre 0.064 et 0.18 d-1 selon groning ou Nurmi
Tref = température de référence à laquelle a été calculé les taux de biodégradation : 298 selon groning
Forcings :
Temp = Température de l’eau en kelvin
La biodégradation est inhibé lorsque l’étang de stabilisation est gelée. Or, Smith (1996) montre que dans une bonne partie du Nunavut, incluant coral harbour, l’épaisseur de la glace atteint environ 1.5 m d’épaisseur. Il est alors possible de considérer que l’hiver, il n’y a pas de biodégradation.
Idéalement, il serait possible d’utiliser l’épaisseur de la glace comme forcing, en considérant que la biodégradation à lieu uniqumement dans la portion qui n’est pas encore gelée et donc multipliant la biodégradtion par la proportion qui n’est pas gelée. Toutefois, déterminer l’épaisseur de la glace demande bien souvent une étude sur le terrain puisque le processus de congélation est complexe et fait intervenir trop de paramètre pour aisément le paramétriser.
Ainsi, je considère que dès que la température de l’eau atteind 0°C (Temp<273), kbio et kadsorp = 0
f3 : Dégradation des produits de phototransformation
Les produits de dégradations primaires issue de la photodégradation sont presque entièrement composée de 8-chlorocarbazole-1-acetic acid (Poiger, 2001) Le taux de dégradation est donc celui de ce composés.
Cette photodégratation dépend de la température et de la lumière comme pour DFC. Cependant, selon Poiger (2001), le taux de dégradation de ce composés se situe entre 4.55 et 4.75 fois supérieur à celui du DFC, variant avec lui selon la lumière et la température,
Ainsi, il est possible pour l’instant de le décrire comme
f3 = kphoto2 * TPphoto
f3 = f1/DFC * 4.6 * TPphoto
f4 : Dégradation des produits de biotransformations
Pour ce qui est de la biodégradation, le principal produit de décomposition par les bactéries retrouvé dans les sédiments est le 5-hydroxydiclofenac (5HDQI). Celui-ci semble être décomposé à un taux constant, dépendant peu de l’activité microbienne subséquente, probablement par adsorption des sédiments. L’effet de la température est mitigié. Normalement, l’adsorption est augmenté lorsque la température diminue, mais la baisse de température augmente la viscosité de l’eau, ce qui compense pour ce processus (Smith, 1996). Il serait possible d’ajouter un terme éventuellement pour déterminer effet de la température et ajouter un forcings.
f4= kadsorp * TPbio
où :
Paramètres :
Kadsorp = constante de réaction d’adsorption du 5-hydroxydiclofenac = 0.051 d-1 selon groning
Kadsorp de log(2)/(13/24)= 1.279656 d-1 selon Nurmi
En résumé :
Équations de flux
f0 <- ((Cinfluent-DFC)* Vinfluent * POP)/Vtotal
f1 <- (1- exp(-QY*I0*4.57*3600)) * exp((- E/(R * Temp))) * DFC *(1 * 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas))
f2 <- Kbio * DFC * exp((Temp-Tref)/10)*log(q10)
f3 <- f1 /DFC * 4.6 * TPphoto
f4 <- kadsorp * TPbio
f5 <- ((0 -TPphoto)* Vinfluent * POP)/Vtotal
f6 <- ((0 -TPbio)* Vinfluent * POP)/VtotalParamètres
| Paramètre | Valeur | Unité |
|---|---|---|
| Cinfluent | 340 | ng/L |
| Vinfluent | 140 | L/p*d |
| Vbassin | 75 000 000 | L |
| POP | 1000 | p |
| A | 1.707 | min−1 |
| E | 5302 | J mol−1 |
| R | 8.314462 | J mol−1 K−1 |
| Iref | 450 | Watt m-2 |
| q10 | 2 | |
| Kbio | 0.064 et 0.18 | d-1 |
| Tref | 298 | K (kelvin) |
| kadsorp | 0.051 | d-1 |
Forcings
| Variable | unités | Source | Pas de temps |
|---|---|---|---|
| Temp | Kelvin | NOAA | |
| I0 | Watt m-2 | CWEC | Heure |
Tester les équations
Vérification des unités
f0 = \(\frac{ng L-1 * L * p-1 * h-1 * p}{ L + L}\) = ng L-1 h-1
f1 = (1-e^(Eintein * EInstein−1) * \(\frac{J mol−1}{J mol−1 K−1 * K}\) * 60 * ng L-1 * 0.13 = ng L-1 h-1
f2 = h-1 * ng L-1 * e^(\(\frac{K-K}{10}\) * K-1) = ng L-1 h-1
f3 = e^(ln(min−1) - \(\frac{J mol−1}{J mol−1 K−1 * K}\) * 60 * ng L-1 * Watt m-2/Watt m-2 * 4. 6 * 0.13 = ng L-1 h-1
f4= h-1 * ng L-1 = ng L-1 h-1
Conservation de la masse
Dans ce modèle, il y a une source, soit les eaux usées (f0, f5 et f6 ) et deux puits, soit la dégradation des sous produits de photodégradation et de biodégradation, f3 et f4. Ainsi, le taux de changement en concentrationdans le système, considérant une stoechiométrie de 1:1 pour les produits de dégradation, doit être égale à f0 - f3 - f4 + f5 + f6 et 2)
Somme des flux = f0 - f1 - f2 + f1 -f3 + f2 - f4 + f5 + f6 = f0 - f3 - f4 + f5 + f6
Implantation du modèle
# load package with the integration routine:
library(deSolve)
library(insol)
#----------------------#
# the model equations: #
#----------------------#
setwd("/Users/jeremieboucher-fontaine/Desktop/Maitrise/Cours /modelisation/Projet Session//")
dataDFC<-read.csv2("Données/NU_CWEC/Coral_harbour_lumière_heure.csv", dec = ".")
dataTemp <- read.csv2("Données/NU_CWEC/Coral_harbour_temperature de l'eau.csv", dec = ".")
data_suisse <- read.csv2("Données/données_lumière_lac_suisse.csv", dec = ".")
data_temp_suisse <- read.csv2("Données/donnees_température_suisse.csv",dec = ".")
#vient de https://en.climate-data.org/north-america/canada/nunavut/iqaluit-1251/
Diclofenac <- function( t, state, parameters ) {
with( as.list( c( state, parameters ) ), { # unpack the state variables & parameters
##### Forcing #####
#TEMPÉRATURE
#pour cycle annuel avec glace
#Temp = Température
#Temp<-approx(dataTemp$Mois,dataTemp$Temp.kelvin, xout = t/731)$y
#Kbio est inclu dans les forcings pour spécifier qu'il n'y a pas de biodégradation lorsque le bassin est gelé.
#if(Temp<273){Kbio<-0}
#else {Kbio<-(0.064/24)}
#if(Temp<273){kadsorp<-0}
#else {kadsorp<-(0.051/24)}
#Pour cycle annuel sans glace
Temp<-approx(dataTemp$Mois,dataTemp$Temp.kelvin, xout = t/731)$y
#Pour cycle suisse
#Temp<-approx(data_temp_suisse$Mois,(data_temp_suisse$Température+273), xout = t/731)$y
#LUMIÈRE
# POur cycle annuel sans glace
I0<-approx(dataDFC$X,dataDFC$X.Global.Horizontal.Radiation..Wh.m2., xout = t)$y
#I0p<-I0*1.10
#Pour cycle avec glace
#if(Temp<273){I0<-0}
#else {I0<-approx(dataDFC$X,dataDFC$X.Global.Horizontal.Radiation..Wh.m2., xout = t)$y}
#pour test avec lac suisse. https://re.jrc.ec.europa.eu/pvg_tools/en/tools.html
#I0<-approx(data_suisse$heure,data_suisse$global.irradiance, xout = t)$y
#Vbassin
Vbassin = sup_ha * 10000 * prof_bas *1000
#Vtotal
Vtotal <- Vbassin + Vinfluent*POP*t
# Flux
f0 <- ((Cinfluent-DFC)* Vinfluent * POP)/Vtotal
f1 <- (1- exp(-QY*I0*4.57*3600)) * exp((- E/(R * Temp))) * DFC *(1 * 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas))
f2 <- Kbio * DFC * exp((Temp-Tref)/10)*log(q10)
f3 <- f1 /DFC * 4.6 * TPphoto
f4 <- kadsorp * TPbio
f5 <- ((0 -TPphoto)* Vinfluent * POP)/Vtotal
f6 <- ((0 -TPbio)* Vinfluent * POP)/Vtotal
# Differential (balance) equations
dDFC <- f0 - f1 - f2
dTPphoto <- f1 - f3 + f5
dTPbio <- f2 - f4 + f6
# Output, packed as a list
list( c( dDFC,
dTPphoto,
dTPbio
),
c( DFC = DFC, # the output variables
TPphoto = TPphoto,
TPbio = TPbio,
I0 = I0,
f0 = f0,
f1 = f1,
f2 = f2,
f3 = f3,
f4 = f4,
f5 = f5,
Vtotal=Vtotal,
Temp = Temp,
Vbassin = Vbassin
)
)
})
} # end of model#-----------------------#
# Paramètre du modèle: #
#-----------------------#
parameters <- c( Cinfluent = 340, #µg/L
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Kbio = 0.064/24, #d-1
Tref = 298, #K (kelvin)
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)Solution Mathématique
#-------------------------#
# Conditions initales: #
#-------------------------#
state <- c( DFC = 7, #
TPphoto = 1,
TPbio = 3
)
#----------------------#
# RUNNING the model: #
#----------------------#
times <-seq(0:8760)
out <-as.data.frame( ode( state, times, Diclofenac, parameters ) )Illustrastion des résultats
#------------------------#
# PLOTTING model output: #
#------------------------#
out$DFC[times = 6500]## [1] 1.748803plot (times,out$DFC,type="l",lwd=2,main="Diclofenac" ,xlab="Mois", ylab="ng/L", xaxt = "n", las = 1)
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
plot (times,out$TPphoto,type="l",lwd=2,main="TPphoto",xlab="Mois", ylab="ng/L", xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
plot (times,out$TPbio,type="l",lwd=2,main="TPbio",xlab="Mois", ylab="ng/L", xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
plot (times,out$I0,type="l",lwd=2,main="I0 " ,xlab="Mois", ylab="W m-2", xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
plot (times,out$Temp,type="l",lwd=2,main="Température " ,xlab="Mois", ylab="Kelvin", xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
abline(h=273)
par( mfrow=c(2,2))
plot (times,out$Vtotal,type="l",lwd=2,main="Vtotal " ,xlab="time, days",ylab="mmolN/m3", ylim = c(0,1.3e+08))
plot (times,out$Vbassin,type="l",lwd=2,main="Vbassin " ,xlab="time, days",ylab="mmolN/m3", )
plot (times,out$f0 ,type="l",lwd=2,main="f0 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f1 ,type="l",lwd=2,main="f1 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f2 ,type="l",lwd=2,main="f2 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f3 ,type="l",lwd=2,main="f3 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f4 ,type="l",lwd=2,main="f4 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f5 ,type="l",lwd=2,main="f5 " ,xlab="time, days",y2lab="ng/l")## Warning in plot.window(...): "y2lab" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "y2lab" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "y2lab" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "y2lab" is not a
## graphical parameter## Warning in box(...): "y2lab" is not a graphical parameter## Warning in title(...): "y2lab" is not a graphical parameter
plot (times,out$f6 ,type="l",lwd=2,main="f6 " ,xlab="time, days",ylab="ng/l")
#Validation
#Validation
data_valid<-read.csv2("Données/NU_CWEC/Donnees_validation.csv", dec = ".")
#Diclofenac
plot (times,out$DFC,type="l",lwd=2,main="Diclofenac" ,xlab="Mois", ylab="ng/L", xaxt = "n",
ylim = c(0,13), )
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
points(data_valid$Heure,data_valid$concentration, col = "forestgreen", cex = 1.2, pch = 16)
lines(data_valid$Jours,data_valid$Pt_Riley, col="firebrick", lwd=2)
lines(data_valid$Jours,data_valid$Pt_A.G, col="blue",lwd=2)
legend( "topright",
c("Mon modèle","Données Poiger et al."),
lwd = 2,
lty = c(1,NA),
pch = c(NA,16),
col = c("black","forestgreen"),
bty = "n" )
#Température et I0
library(scales)
plot (times,out$I0,type="l",lwd=2, axes = F, xlab = "", ylab = "", col = alpha("blue", 0.4))
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
par(new=TRUE)
plot (times,out$Temp,type="l",lwd=2, axes = F, xlab = "", ylab = "", xaxt = "n", col="firebrick")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
par(new=TRUE)
plot (times,out$DFC,type="l",lwd=2,main="Diclofenac" ,xlab="Mois", ylab="ng/L", xaxt = "n",
ylim = c(0,10), las = 1)
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
legend( "topright",
c("Diclofenac prédit","Éclairement", "Température"),
lwd = 2,
lty = c(1,1,1),
col = c("black",alpha("blue", 0.4),"firebrick"),
bty = "n" )
## Analyses de sensitivité
Cinfluent = 340 #µg/L
Vinfluent = 140/24 #L/p*d
E = 5302
POP = 1000 #p
q10 = 2 #aucune
Kbio = 0.064/24 #d-1
kadsorp = 0.051/24 #d-1
Tref = 298 #K (kelvin)
QY = 0.0375 # mol E-1
sup_ha = 5 #ha
prof_bas = 1.5 #m
#--- Vary parameters one-at-a-time (OAT) by 10%
ss <- 0.3
Cinfluentp <- Cinfluent * (1+ss) # taux de produxction primaire
Vinfluentp <- Vinfluent * (1+ss) # taux de broutage
POPp <- POP * (1+ss) # taux de respiration
q10p <- q10 * (1+ss) #Producteurs primaires initiale
Kbiop <- Kbio * (1+ss)
kadsorpp <- kadsorp * (1+ss)
QYp <- QY * (1+ss)
sup_hap <- sup_ha * (1+ss)
prof_basp <- prof_bas * (1+ss)
Ep <- E * (1+ss)
#--- Run simulations
# référence
parametersref <- c( Cinfluent = 340, #µg/L
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outref <-as.data.frame( ode( state, times, Diclofenac, parametersref ) )
outref <- outref$DFC
# Cinfluent
parametersCinfluent <- c( Cinfluent = Cinfluentp, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
Tref = 298, #K (kelvin)
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outCinfluent <-as.data.frame( ode( state, times, Diclofenac, parametersCinfluent) )
outCinfluent <- outCinfluent$DFC
# Vinfluent
parametersVinfluent <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = Vinfluentp, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outVinfluent <-as.data.frame( ode( state, times, Diclofenac, parametersVinfluent) )
outVinfluent <- outVinfluent$DFC
# POP
parametersPOP <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = POPp, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outPOP <-as.data.frame( ode( state, times, Diclofenac, parametersPOP) )
outPOP <- outPOP$DFC
# q10
parametersq10 <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = q10p, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outq10 <-as.data.frame( ode( state, times, Diclofenac, parametersq10) )
outq10 <- outq10$DFC
# Kbio
parametersKbio <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = Kbiop, #d-1
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outKbio <-as.data.frame( ode( state, times, Diclofenac, parametersKbio) )
outKbio <- outKbio$DFC
# kadsorp
parameterskadsorp <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = kadsorpp, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outkadsorp <-as.data.frame( ode( state, times, Diclofenac, parameterskadsorp) )
outkadsorp <- outkadsorp$DFC
#QY
parametersQY <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
QY = QYp, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outQY <-as.data.frame( ode( state, times, Diclofenac, parametersQY) )
outQY <- outQY$DFC
#E
parametersE <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = Ep, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outE <-as.data.frame( ode( state, times, Diclofenac, parametersE) )
outE <- outE$DFC
#sup_ha
parameterssup_ha <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = sup_hap, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outsup_ha <-as.data.frame( ode( state, times, Diclofenac, parameterssup_ha) )
outsup_ha <- outsup_ha$DFC
#prof_bas
parametersprof_bas <- c( Cinfluent = 340, #g C (g C)-1 d-1 (g cal cm-2 min-1)-1 , # day^-1
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Tref = 298, #K (kelvin)
Kbio = 0.064/24, #d-1
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = prof_basp #m
)
state <- c( DFC = 7, TPphoto = 1,TPbio = 6)
times <-seq(0:8760)
outprof_bas <-as.data.frame( ode( state, times, Diclofenac, parametersprof_bas) )
outprof_bas <- outprof_bas$DFC#Variation d'arrivée de DFC
plot( times,
outref,
type = "l",
lwd = 1,
ylim = c(0,10),
xlab = "Temps (Mois)",
ylab = "Concentration de diclofenac, ng/l",
main = "Sensibilité locale\nVariation de l'apport de DFC",
xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
lines( times,
outVinfluent,
ylim = c(0,10),
lwd = 2,
lty = 2,
col="red")
lines( times,
outCinfluent,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="blue")
lines( times,
outPOP,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="green")
legend( "topright",
c("Référence","Vinfluent + 30%", "Cinfluent + 30%", "POP + 30%"),
lwd = 2,
lty = c(1,1,1),
col = c("black","red","blue","green"),
bty = "n" )
# Variation de la photodégradation
plot( times,
outref,
type = "l",
lwd = 1,
ylim = c(0,10),
xlab = "Temps (Mois)",
ylab = "Concentration de diclofenac, ng/l",
main = "Sensibilité locale\nVariation de la photodégradation",
xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
lines( times,
outQY,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="red")
lines( times,
outE,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="blue")
legend( "topright",
c("Référence","Quantum yield + 30%", "Énergie d'Activation (E) + 30%"),
lwd = 2,
lty = c(1,1,1),
col = c("black","red","blue"),
bty = "n" )
# Variation de la Biodégradation
plot( times,
outref,
type = "l",
lwd = 1,
ylim = c(0,10),
xlab = "Temps (Mois)",
ylab = "Concentration de diclofenac, ng/l",
main = "Sensibilité locale\nVariation de la biodégradation ",
xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
lines( times,
outKbio,
ylim = c(0,10),
lwd = 2,
lty = 2,
col="#FF9900")
lines( times,
outq10,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="#33FFFF")
legend( "topright",
c("Référence","kbio + 30%", "q10 + 30%"),
lwd = 2,
lty = c(1,1,1),
col = c("black","#FF9900","#33FFFF"),
bty = "n" )
# Variation de la taille du bassin
plot( times,
outref,
type = "l",
lwd = 1,
ylim = c(0,10),
xlab = "Temps (Mois)",
ylab = "Concentration de diclofenac, ng/l",
main = "Sensibilité locale\nVariation de la taille du bassin ",
xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
lines( times,
outsup_ha,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="#336666")
lines( times,
outprof_bas,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="#CC0033")
legend( "topright",
c("Référence","Sup_ha + 30%", "prof_ha + 30%"),
lwd = 2,
lty = c(1,1,1),
col = c("black","#336666","#CC0033"),
bty = "n" )
################# Effet de la lumière #################
DiclofenacI0add <- function( t, state, parameters ) {
with( as.list( c( state, parameters ) ), { # unpack the state variables & parameters
Temp<-approx(dataTemp$Mois,dataTemp$Temp.kelvin, xout = t/731)$y
I0add<-approx(dataDFC$X,dataDFC$X.Global.Horizontal.Radiation..Wh.m2., xout = t)$y
I0addp<-I0add+10
Vbassin = sup_ha * 10000 * prof_bas *1000
Vtotal <- Vbassin + Vinfluent*POP*t
# Flux
f0 <- ((Cinfluent-DFC)* Vinfluent * POP)/Vtotal
f1 <- (1- exp(-QY*I0addp*4.57*3600)) * exp((- E/(R * Temp))) * DFC *(1 * 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas))
f2 <- Kbio * DFC * exp((Temp-Tref)/10)*log(q10)
f3 <- f1 /DFC * 4.6 * TPphoto
f4 <- kadsorp * TPbio
f5 <- ((0 -TPphoto)* Vinfluent * POP)/Vtotal
f6 <- ((0 -TPbio)* Vinfluent * POP)/Vtotal
dDFC <- f0 - f1 - f2
dTPphoto <- f1 - f3 + f5
dTPbio <- f2 - f4 + f6
list( c( dDFC,
dTPphoto,
dTPbio
),
c( DFCI0add = DFC # the output variables
)
)
})
}
parameters <- c( Cinfluent = 340, #µg/L
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Kbio = 0.064/24, #d-1
Tref = 298, #K (kelvin)
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1, TPbio = 6 )
times <-seq(0:8760)
outI0add <-as.data.frame( ode( state, times, DiclofenacI0add, parameters ) )# end of model
outI0add <- outI0add$DFCI0add
DiclofenacI0 <- function( t, state, parameters ) {
with( as.list( c( state, parameters ) ), { # unpack the state variables & parameters
Temp<-approx(dataTemp$Mois,dataTemp$Temp.kelvin, xout = t/731)$y
I0<-approx(dataDFC$X,dataDFC$X.Global.Horizontal.Radiation..Wh.m2., xout = t)$y
I0p<-I0*10
Vbassin = sup_ha * 10000 * prof_bas *1000
Vtotal <- Vbassin + Vinfluent*POP*t
# Flux
f0 <- ((Cinfluent-DFC)* Vinfluent * POP)/Vtotal
f1 <- (1- exp(-QY*I0p*4.57*3600)) * exp((- E/(R * Temp))) * DFC *(1 * 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas))
f2 <- Kbio * DFC * exp((Temp-Tref)/10)*log(q10)
f3 <- f1 /DFC * 4.6 * TPphoto
f4 <- kadsorp * TPbio
f5 <- ((0 -TPphoto)* Vinfluent * POP)/Vtotal
f6 <- ((0 -TPbio)* Vinfluent * POP)/Vtotal
dDFC <- f0 - f1 - f2
dTPphoto <- f1 - f3 + f5
dTPbio <- f2 - f4 + f6
list( c( dDFC,
dTPphoto,
dTPbio
),
c( DFCI0 = DFC # the output variables
)
)
})
}
parameters <- c( Cinfluent = 340, #µg/L
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Kbio = 0.064/24, #d-1
Tref = 298, #K (kelvin)
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1, TPbio = 6 )
times <-seq(0:8760)
outI0 <- as.data.frame( ode( state, times, DiclofenacI0, parameters ) )# end of model
outI0 <- outI0$DFCI0
########## Effet de la température ############33
DiclofenacTemp <- function( t, state, parameters ) {
with( as.list( c( state, parameters ) ), { # unpack the state variables & parameters
Temp<-approx(dataTemp$Mois,dataTemp$Temp.kelvin, xout = t/731)$y
Tempp<- Temp * 1.10
I0<-approx(dataDFC$X,dataDFC$X.Global.Horizontal.Radiation..Wh.m2., xout = t)$y
Vbassin = sup_ha * 10000 * prof_bas *1000
Vtotal <- Vbassin + Vinfluent*POP*t
# Flux
f0 <- ((Cinfluent-DFC)* Vinfluent * POP)/Vtotal
f1 <- (1- exp(-QY*I0*4.57*3600)) * exp((- E/(R * Tempp))) * DFC *(1 * 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas))
f2 <- Kbio * DFC * exp((Tempp-Tref)/10)*log(q10)
f3 <- f1 /DFC * 4.6 * TPphoto
f4 <- kadsorp * TPbio
f5 <- ((0 -TPphoto)* Vinfluent * POP)/Vtotal
f6 <- ((0 -TPbio)* Vinfluent * POP)/Vtotal
dDFC <- f0 - f1 - f2
dTPphoto <- f1 - f3 + f5
dTPbio <- f2 - f4 + f6
list( c( dDFC,
dTPphoto,
dTPbio
),
c( DFCTemp = DFC # the output variables
)
)
})
}
parameters <- c( Cinfluent = 340, #µg/L
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Kbio = 0.064/24, #d-1
Tref = 298, #K (kelvin)
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1, TPbio = 6 )
times <-seq(0:8760)
outTemp <-as.data.frame( ode( state, times, DiclofenacTemp, parameters ) )# end of model
outTemp <-outTemp$DFCTemp
DiclofenacTempmoins <- function( t, state, parameters ) {
with( as.list( c( state, parameters ) ), { # unpack the state variables & parameters
Tempmoins<-approx(dataTemp$Mois,dataTemp$Temp.kelvin, xout = t/731)$y
Tempmoinsp<- Tempmoins * 0.90
I0<-approx(dataDFC$X,dataDFC$X.Global.Horizontal.Radiation..Wh.m2., xout = t)$y
Vbassin = sup_ha * 10000 * prof_bas *1000
Vtotal <- Vbassin + Vinfluent*POP*t
# Flux
f0 <- ((Cinfluent-DFC)* Vinfluent * POP)/Vtotal
f1 <- (1- exp(-QY*I0*4.57*3600)) * exp((- E/(R * Tempmoinsp))) * DFC *(1 * 0.10 + 0.1 * 0.90 + 0.01 * (1 - prof_bas))
f2 <- Kbio * DFC * exp((Tempmoinsp-Tref)/10)*log(q10)
f3 <- f1 /DFC * 4.6 * TPphoto
f4 <- kadsorp * TPbio
f5 <- ((0 -TPphoto)* Vinfluent * POP)/Vtotal
f6 <- ((0 -TPbio)* Vinfluent * POP)/Vtotal
dDFC <- f0 - f1 - f2
dTPphoto <- f1 - f3 + f5
dTPbio <- f2 - f4 + f6
list( c( dDFC,
dTPphoto,
dTPbio
),
c( DFCTempmoins = DFC # the output variables
)
)
})
}
parameters <- c( Cinfluent = 340, #µg/L
Vinfluent = 140/24, #L/p*d
POP = 1000, #p
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
q10 = 2, #aucune
Kbio = 0.064/24, #d-1
Tref = 298, #K (kelvin)
kadsorp = 0.051/24, #d-1
QY = 0.0375, # mol E-1
sup_ha = 5, #ha
prof_bas = 1.5 #m
)
state <- c( DFC = 7, TPphoto = 1, TPbio = 6 )
times <-seq(0:8760)
outTempmoins <-as.data.frame( ode( state, times, DiclofenacTempmoins, parameters ) )# end of model
outTempmoins <- outTempmoins$DFCTempmoins#Variation de la lumière
plot (times,outref,type="l",lwd=2,main="Analyses de sensibilité\nVariation de l'éclairement" ,xlab="Mois", ylab="ng/L", xaxt = "n",
ylim = c(0,10), )
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
lines( times,
outI0,
ylim = c(0,10),
lwd = 1,
col = "firebrick",
lty = 2
)
lines( times,
outI0add,
ylim = c(0,10),
lwd = 1,
col = "forestgreen",
lty = 3,
)
legend( "topright",
c("Référence","Éclairement + 10%","Éclairement + 10 watt/m2 "),
lwd = 2,
lty = c(1,2,3),
col = c("black","firebrick","forestgreen"),
bty = "n" )
#Variation de la température
plot (times,outref,type="l",lwd=2,main="Analyses de sensibilité\nVariation de la température" ,xlab="Mois", ylab="ng/L", xaxt = "n",
ylim = c(0,10), )
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
lines( times,
outTemp,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="red")
lines( times,
outTempmoins,
ylim = c(0,10),
lwd = 1,
lty = 2,
col="blue")
legend( "topright",
c("Référence","Température + 10%", "Température - 10%"),
lwd = 2,
lty = c(1,1,1),
col = c("black","red","blue"),
bty = "n" )
# Average response from reference simulation
outm <- mean(outref, na.rm = TRUE)
str(outref)## num [1:8761] 7 7.02 7.05 7.07 7.1 ... Cinfluent = 340 #µg/L
Vinfluent = 140/24 #L/p*d
E = 5302
POP = 1000 #p
q10 = 2 #aucune
Kbio = 0.064/24 #d-1
kadsorp = 0.051/24 #d-1
Tref = 298 #K (kelvin)
QY = 0.0375 # mol E-1
sup_ha = 5 #ha
prof_bas = 1.5 #m
# Local sensitivity S.oat (array) & Sij (scalar)
S.oat <- cbind( Vinfluent = (outVinfluent-outref)/outm/ss,
POP = (outPOP-outref)/outm/ss,
Cinfluent = (outCinfluent-outref)/outm/ss,
q10 = (outq10-outref)/outm/ss,
Kbio = (outKbio-outref)/outm/ss,
kadsorp = (outkadsorp-outref)/outm/ss,
QY = (outQY-outref)/outm/ss,
E = (outE-outref)/outm/ss,
sup_ha = (outsup_ha-outref)/outm/ss,
prof_bas = (outprof_bas-outref)/outm/ss,
Temperature = (outTemp-outref)/outm/ss,
Lumiere = (outI0-outref)/outm/ss
)
colnames(S.oat) <- c( "Volume d'influent/personne", "Population", "Concentration dans influent","q10","k de biodégradation","k d'adsorption","Quantum Yield","E","superficie du bassin (ha)" ,"Profondeur du bassin","Température", "Éclairement")
n <- length(times)
Sij <- ( apply( S.oat^2, 2, sum, na.rm = TRUE )/n )^0.5
write.csv2(Sij, file = "sensibilite.csv" , dec = ".")## Warning in write.csv2(Sij, file = "sensibilite.csv", dec = "."): attempt to set
## 'dec' ignored#--- Figures
# Taux de production primaire (photosynthèse)Section pour le modèle du lac Greifensee (suisse)
# load package with the integration routine:
library(deSolve)
library(insol)
#----------------------#
# the model equations: #
#----------------------#
setwd("/Users/jeremieboucher-fontaine/Desktop/Maitrise/Cours /modelisation/Projet Session//")
data_suisse <- read.csv2("Données/données_lumière_lac_suisse.csv", dec = ".")
data_temp_suisse <- read.csv2("Données/donnees_température_suisse.csv",dec = ".")
#vient de https://en.climate-data.org/north-america/canada/nunavut/iqaluit-1251/
Diclofenac <- function( t, state, parameters ) {
with( as.list( c( state, parameters ) ), { # unpack the state variables & parameters
##### Forcing #####
#TEMPÉRATURE
Temp<-approx(data_temp_suisse$Mois,(data_temp_suisse$Température+273), xout = t/731)$y
#LUMIÈRE
#https://re.jrc.ec.europa.eu/pvg_tools/en/tools.html
I0<-approx(data_suisse$heure,data_suisse$global.irradiance, xout = t)$y
#Vbassin du lac
Q = Vbassin/t_res
Ain = (Intrant)/Q
#arruvé de DFC p
#2.5e+10 /24 = 1041666667 h-1
#ATTÉNUATION PAR LA PROFONDEUR
# profondeur du lac de 17.8 m
#(1 * 0.10 + 0.1 * 0.90 + 0.01 * 16.8 ) / 17.8 = 0.02011236
#si prof de 33 m (max)
# (1 * 0.10 + 0.1 * 0.90 + 0.01 * 32 ) / 33= 0.01545455
# Flux
f0 <- 1/(t_res)*(Ain-DFC)
f1 <- (1- exp(-QY*I0*4.57*3600)) * exp((- E/(R * Temp))) * DFC * 0.02011236
f2 <- Kbio * DFC * exp((Temp-Tref)/10)*log(q10)
f3 <- f1 /DFC * 4.6 * TPphoto
f4 <- kadsorp * TPbio
# Differential (balance) equations
dDFC <- f0 - f1 - f2
dTPphoto <- f1 - f3
dTPbio <- f2 - f4
# Output, packed as a list
list( c( dDFC,
dTPphoto,
dTPbio
),
c( DFC = DFC, # the output variables
TPphoto = TPphoto,
TPbio = TPbio,
I0 = I0,
f0 = f0,
f1 = f1,
f2 = f2,
f3 = f3,
f4 = f4,
Temp = Temp,
Vbassin = Vbassin
)
)
})
} # end of model#-----------------------#
# Paramètre du modèle: #
#-----------------------#
parameters <- c( Intrant = 3.5e+10 /24, #ng/h
t_res = 408*24, #L/p*h
Vbassin = 1.51e+11, #L
E = 5302, #J mol−1
R = 8.314462, #J mol−1 K−1
Iref = 450, #Watt m-2
q10 = 2, #aucune
Kbio = 0, #d-1
#Kbio = 0.064/24, #d-1
Tref = 298, #K (kelvin)
kadsorp = 0.051/24, #d-1
kphoto1 = 2.6, #h-1, test selon poiger
QY = 0.0375, # mol E-1
sup_ha = 5 #ha
)Solution Mathématique
#-------------------------#
# Conditions initales: #
#-------------------------#
state <- c( DFC = 11, #
TPphoto = 2,
TPbio = 6
)
#----------------------#
# RUNNING the model: #
#----------------------#
times <-seq(0:8760)
out <-as.data.frame( ode( state, times, Diclofenac, parameters ) )Illustrastion des résultats
#------------------------#
# PLOTTING model output: #
#------------------------#
out$DFC[times = 6500]## [1] 6.757797plot (times,out$DFC,type="l",lwd=2,main="Diclofenac" ,xlab="Mois", ylab="ng/L", xaxt = "n",
)
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
plot (times,out$TPphoto,type="l",lwd=2,main="8-chlorocarbazole-1-acetic acid",xlab="Mois", ylab="ng/L", xaxt = "n", ylim = c(0,4), col="forestgreen")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
legend( "topright",
c("Mon modèle","Poiger et al. (épi)","Poiger et al. (hypo)"),
lwd = 2,
lty = c(1,1,2),
col = c("forestgreen","black","black"),
bty = "n" )
plot (times,out$TPbio,type="l",lwd=2,main="TPbio",xlab="Mois", ylab="ng/L", xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
plot (times,out$I0,type="l",lwd=2,main="I0 " ,xlab="Mois", ylab="W m-2", xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
plot (times,out$Temp,type="l",lwd=2,main="Température " ,xlab="Mois", ylab="Kelvin", xaxt = "n")
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
abline(h=273)
par( mfrow=c(2,2))
plot (times,out$Vtotal,type="l",lwd=2,main="Vtotal " ,xlab="time, days",ylab="mmolN/m3")
plot (times,out$Vbassin,type="l",lwd=2,main="Vbassin " ,xlab="time, days",ylab="mmolN/m3")
plot (times,out$f0 ,type="l",lwd=2,main="f0 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f1 ,type="l",lwd=2,main="f1 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f2 ,type="l",lwd=2,main="f2 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f3 ,type="l",lwd=2,main="f3 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f4 ,type="l",lwd=2,main="f4 " ,xlab="time, days",ylab="ng/l")
plot (times,out$f5 ,type="l",lwd=2,main="f5 " ,xlab="time, days",y2lab="ng/l")## Warning in plot.window(...): "y2lab" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "y2lab" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "y2lab" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "y2lab" is not a
## graphical parameter## Warning in box(...): "y2lab" is not a graphical parameter## Warning in title(...): "y2lab" is not a graphical parameter
plot (times,out$f6 ,type="l",lwd=2,main="f6 " ,xlab="time, days",ylab="ng/l")
#Validation
#Validation
data_valid<-read.csv2("Données/NU_CWEC/Donnees_validation.csv", dec = ".")
#Diclofenac
plot (times,out$DFC,type="l",lwd=2,main="Validation" ,xlab="Mois", ylab="Concentration de diclofenac (ng/L)", xaxt = "n",
ylim = c(0,14), las = 1, col = "forestgreen" )
axis (1,
at = c(0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016,8760),
labels = c("J", "F", "M", "A", "M","J","J","A","S","O","N","D","J"))
points(data_valid$Heure,data_valid$concentration, col = "black", cex = 1.2, pch = 1)
legend( "topright",
c("Mon modèle","Modèle de Poiger et al.","Données de Poiger et al."),
lwd = 2,
lty = c(1,1,NA),
pch = c(NA,NA,1),
col = c("forestgreen","black","black"),
bty = "n" )