[R-sig-dyn-mod] ode integration troubles
Romain PETE
Romain.Pete at ifremer.fr
Wed Feb 3 17:24:54 CET 2016
Dear all,
I'm really struggling to get this biogeochemical model to work, main
issue, I never get a full integration (over a year period) despite
trying many different time steps (or variable time step solvers). I've
checked conservativity, steady state conditions and parameters... any
help would be much appreciated
This biogeochemical model is to be included in a LOICZ type of model
however in the following code, no inputs are yet considered
Many thanks
romain
#-----------------------#
# Model equations: #
#-----------------------#
BIOGEO6 <- function(t,state,parameters) {
with(as.list(c(state,parameters)), { # unpack the state variables
## LIGHT & TEMPERATURE
# Computed environmental variables
I0 <- 0.48*(540+440*sin(2*pi*t/365-1.4))
watertemp <- function(Time, A, C, W, O){A * cos(pi * (Time - C) /
W) + O}
Temp <- watertemp(Time=t, A=6, C=-131, W=182, O=16)
# Safe state variables: function max used to reset low values to zero
phytoS <- max(PhytoS, 0)
phytoL <- max(PhytoL, 0)
zs <- max(ZS, 0)
zl <- max(ZL, 0)
ulva <- max(ULVA, 0)
graci <- max(GRACI, 0)
zost <- max(ZB, 0)
qN_ulva <- max(Q_N_Ulva, 0)
qN_graci <- max(Q_N_Graci, 0)
qN_zost <- max(Q_N_Zost, 0)
qP_ulva <- max(Q_P_Ulva, 0)
qP_graci <- max(Q_P_Graci, 0)
qP_zost <- max(Q_P_Zost, 0)
l_pon <- max(l_PON, 0)
l_pop <- max(l_POP, 0)
r_pon <- max(r_PON, 0)
r_pop <- max(r_POP, 0)
l_don <- max(l_DON, 0)
l_dop <- max(l_DOP, 0)
r_don <- max(r_DON, 0)
r_dop <- max(r_DOP, 0)
no3 <- max(0, NO3)
nh4 <- max(0, NH4)
po4 <- max(0, PO4)
sedN <- max(N_sed, 0)
sedP <- max(P_sed, 0)
burialN <- max(N_burial, 0)
burialP <- max(P_burial, 0)
## PHYTOPLANKTON (d'aprés Chapelle 2000)
# Phytoplankton computed parameters: Light, temperature
Chla <- (phytoS + phytoL) / CHLaN_ratio
I_phyto <- I0*exp(-kI * z_phyto)
f_I_PS <- 1- exp(-(I_phyto/Iopt_PS))
f_I_PL <- 1 - exp(-(I_phyto/Iopt_PL))
f_Temp_PS <- 1/(1+exp(-zeta_PS*(Temp-Temp_PS)))
f_Temp_PL <- 1/(1+exp(-zeta_PL*(Temp-Temp_PL)))
# Uptake limitation functions
f_N_PS <- ((no3/(no3+K_N_PS)) * exp(-phi_PS * nh4)) +
(nh4/(nh4+K_N_PS)) # Ammonia preference introduced with correction
factor Phi
f_P_PS <- (po4/(po4+K_P_PS))
f_N_PL <- ((no3+nh4)/(no3+nh4+K_N_PL))
f_P_PL <- (po4/(po4+K_P_PL))
f_SN_PS <- min(f_N_PS, f_P_PS)
f_SN_PL <- min(f_N_PL, f_P_PL)
# growth rate of phytoplankton
mu_PS <- mu_max_PS * f_Temp_PS * f_I_PS * f_SN_PS
mu_PL <- mu_max_PL * f_Temp_PL * f_I_PL * f_SN_PL
# Uptakes: N & P with NH4 preference feedback control for PS
NO3_uptake_phyto <- mu_PS * (((no3/(no3+K_N_PS)) * exp(-phi_PS *
nh4))/f_N_PS) * phytoS + mu_PL * (no3/(nh4 + no3)) * phytoL
NH4_uptake_phyto <- mu_PS * ((nh4/(nh4+K_N_PS))/f_N_PS) * phytoS
+ mu_PL * (nh4/(nh4 + no3)) * phytoL
PO4_uptake_phyto <- mu_PS * phytoS * (po4/(po4+K_P_PS)) + mu_PL *
phytoL * (po4/(po4+K_P_PL))
# Growth
Growth_PS <- mu_PS * phytoS
Growth_PL <- mu_PL * phytoL
# Lysis
Lysis_PS <- K_lys_PS * phytoS
Lysis_PL <- K_lys_PL * phytoL
# Exudation
Exudation_PS <- K_e_PS * phytoS
Exudation_PL <- K_e_PL * phytoL
# Sedimentation
Sed_PL <- vsed_PL*phytoL
# IN-OUT: ADD Vx * Phyto for open sea-lagoon exchange
## ZOOPLANKTON (d'aprés Zaldivar 2003)
# Zooplankton computed parameters
gmax_ZS <- g_prim_max_ZS*exp(-((Temp-Temp_opt_ZS)/Temp_width_ZS)^2)
gmax_ZL <- g_prim_max_ZL*exp(-((Temp-Temp_opt_ZL)/Temp_width_ZL)^2)
#K_app_B <- K_B*((1+PS)/K_PS)
K_app_PS <- K_PS #*(1+(B/K_B)) #
K_app_PL <- K_PL*(1+(zs/K_ZS))
K_app_ZS <- K_ZS*(1+(phytoL/K_PL))
# Zoo Grazing
#grazing_B <- gmax_ZS*ZS*((B)/(B+K_app_B)) # No Bacteria in this
version of the model
grazing_PS <- gmax_ZS*zs*(phytoS/(phytoS+K_app_PS))
grazing_ZS_tot <- grazing_PS #+grazing_B
grazing_PL <- gmax_ZL*zl*(phytoL/(phytoL+K_app_PL))
grazing_ZS <- gmax_ZL*zl*(zs/(zs+K_app_ZS))
grazing_ZL_tot <- grazing_PL + grazing_ZS
# Zoo Growth
growth_ZS <- Y_ZS*grazing_ZS_tot
growth_ZL <- Y_ZL*grazing_ZL_tot
# Zoo egestion
egestion_ZL <- E_ZL*grazing_ZL_tot
# Zoo lysis (mortality)
Lysis_ZS <- K_d_ZS*zs
Lysis_ZL <- K_d_ZL*zl
## MACROPHYTES
## Ulva: ULVA (d'aprés Zaldivar 2003 et Coffaro 1997)
# Ulva accesory functions
I_Ulva <- I0*exp(-kI_Ulva*z_fond)
f_I_Ulva <- ((I_Ulva/Iopt_Ulva)*exp(1-(I_Ulva/Iopt_Ulva)))
f_Temp_Ulva <- 1/(1+exp(-zeta_Ulva*(Temp-Temp_Ulva)))
# Ulva nutrient quota
f_QN_Ulva <- (qN_ulva - Qmin_N_Ulva)/(Qmax_N_Ulva-Qmin_N_Ulva)
f_QP_Ulva <- (qP_ulva - Qmin_P_Ulva)/(Qmax_P_Ulva-Qmin_P_Ulva)
# Ulva uptake of N & P
NO3_uptake_Ulva <- V_NO3_Ulva*(no3/(no3+K_NO3_Ulva)) * (1 - f_QN_Ulva)
NH4_uptake_Ulva <- V_NH4_Ulva*(nh4/(nh4+K_NH4_Ulva)) * (1 - f_QN_Ulva)
PO4_uptake_Ulva <- V_PO4_Ulva*(po4/(po4+K_PO4_Ulva)) * (1 - f_QP_Ulva)
# Lim function for nutrients
f_SN_Ulva <- min(f_QN_Ulva, f_QP_Ulva)
mu_Ulva <- mu_max_Ulva * f_I_Ulva * f_Temp_Ulva * f_SN_Ulva
# Ulva growth
growth_Ulva <- mu_Ulva * ulva
# Ulva death
death_Ulva <- MaxMort_Ulva*ulva^beta_Ulva
#death_Ulva <- MaxMort_Ulva * exp(Temp - 27)
## Gracilaria: GRACI (calqué sur le modéle Ulve de Zaldivar 2003)
# Gracilaria accesory functions
I_Graci <- I0*exp(-kI_Graci*z_fond)
f_I_Graci <- (I_Graci/Iopt_Graci) * exp(1-(I_Graci/Iopt_Graci))
f_Temp_Graci <- 1/(1+exp(-zeta_Graci*(Temp-Temp_Graci)))
f_QN_Graci <- (qN_graci - Qmin_N_Graci)/(Qmax_N_Graci-Qmin_N_Graci)
f_QP_Graci <- (qP_graci - Qmin_P_Graci)/(Qmax_P_Graci-Qmin_P_Graci)
# Gracilaria uptake of N & P
NO3_uptake_Graci <- V_NO3_Graci*(no3/(no3+K_NO3_Graci)) * (1 -
f_QN_Graci)
NH4_uptake_Graci <- V_NH4_Graci*(nh4/(nh4+K_NH4_Graci)) * (1 -
f_QN_Graci)
PO4_uptake_Graci <- V_PO4_Graci*(po4/(po4+K_PO4_Graci)) * (1 -
f_QP_Graci)
# Lim function for nutrients
f_SN_Graci <- min(f_QN_Graci, f_QP_Graci)
mu_Graci <- mu_max_Graci * f_I_Graci * f_Temp_Graci * f_SN_Graci
# Gracilaria growth
growth_Graci <- mu_Graci * graci
# Gracilaria death
death_Graci <- MaxMort_Graci*graci^beta_Graci
#death_Graci <- MaxMort_Graci * exp(Temp - 27)
## Zostera: ZB (d'aprés Plus 2003 et discussion Plus/Pete Nov 2015)
# Photosynthesis of Zostera
Zost_Pmax <- theta_Pmax * Temp - Pmax_0C
Qcan <- I0 * exp(-K1*(z_fond-0.3))
f_can <- exp(-K3 * z_zost)
PAR_k1 <- 0.5*(Ikmax + Ikmin) + 0.5*(Ikmax - Ikmin) *
cos(2*pi*((t-110)/365))
f_I_Zost <- tanh((Qcan*f_can)/PAR_k1)
Zost_Ptot <- Zost_Pmax * f_I_Zost #* z_zost
# Zostera growth
f_QN_Zost <- (qN_zost - LN_min) / (LN_max - LN_min)
f_QP_Zost <- (qP_zost - LP_min) / (LP_max - LP_min)
f_LN <- f_QN_Zost^epsilon1_zost
f_LP <- f_QP_Zost^epsilon2_zost
f_SN_zost <- min(f_LN, f_LP)
# Zostera uptake of N & P
LABS_NH4 <- LVmaxN * (nh4/(nh4 + KL_zost_NH4)) * (1 -
f_QN_Zost)^delta2_zost
LABS_NO3 <- LVmaxN * (no3/(no3 + KL_zost_NO3)) * (1 -
f_QN_Zost)^delta1_zost
LABS_PO4 <- LVmaxP * (po4/(po4 + KL_zost_PO4)) * (1 -
f_QP_Zost)^delta3_zost
ZGR <- max(Zost_Ptot * f_SN_zost, 0)
# Zostera mortality
ZM_fT <- theta_LM^(Temp-20)
#ZM_fV <- LMR_V * VitVent/10 * exp(-K4*z_zost)
ZM <- LMR_20C * ZM_fT #+ ZM_fV
## ORGANIC MATTER (simplifié d'aprés Lancelot 2002)
# Particulate organic matter
# Hydrolysis of POM
# Hydrolysis constante of labile POM as a function of temperature
K_lP <- K_lP_max * exp(-((Temp-Temp_opt_OM)/Temp_width_OM)^2)
# Hydrolysis constante of refractory POM as a function of temperature
K_rP <- K_rP_max * exp(-((Temp-Temp_opt_OM)/Temp_width_OM)^2)
# Biological lysis
Lys_Bio_N <- Lysis_PS + Lysis_PL + Lysis_ZS + Lysis_ZL #
expressed in N
Lys_Bio_P <- (Lysis_PS + Lysis_PL)/NP_phyto + (Lysis_ZS +
Lysis_ZL)/NP_zoo # expressed in P (mmolP/m3)
# Dissolved Organic Matter
# Hydrolysis of DOM, temperature dependant
# Max hydrolysis rate of labile DOM
V_lD <- Vmax_lD * exp(-((Temp-Temp_opt_OM)/Temp_width_OM)^2)
# Max hydrolysis rate of refractory DOM
V_rD <- Vmax_rD * exp(-((Temp-Temp_opt_OM)/Temp_width_OM)^2)
# Hydrolysis of labile DON to NH4 following M-M curve
elys_l_DON <- V_lD * (l_don/(l_don + K_l))
# Hydrolysis of refractory DON to NH4 (slower)
elys_r_DON <- V_rD * (r_don/(r_don + K_r))
# Hydrolysis of labile DOP to PO4 (follows M-M expression)
elys_l_DOP <- V_lD * (l_dop/(l_dop + K_l))
# Hydrolysis of refractory DOP to PO4
elys_r_DOP <- V_rD * (r_dop/(r_dop + K_r))
## SEDIMENT (determiner a partir de donnée terrain Ifremer, V.
Ouisse, projet DEPART & RESTOLAG)
# Empirical relationships between Temperature, stocks (N, P) and
fluxes of NH4, NO3, PO4, DON & DOP
# See Sediment_data.R with LM modelling
FxNH4 <- max(0, (0.43785*sedN - 54.99)*1e-3/14) # je passe de mg en
g puis de g en mol
FxNO3 <- max(0, (-0.037205*sedN + 4.39)*1e-3/14)
#
FxPO4 <- max(0, (0.306*Temp + 0.0318*sedP - 6.75)*1e-3/31)
# FxPO4_t <- 0.0256 * Temp - 0.1674
#
FxDON <- max(((0.63*Temp + 0.29)*1e-3)/14, 0) ## relation non
significative !!!!
FxDOP <- max(((0.62467*Temp - 17.53)*1e-3)/31, 0)
#Qt?s N (mmol /m?/jr) = (-0.0102*Temp?rature?) +
(0.5172*Temp?rature) - 3.264
#Qt?s P (mmol /m?/jr) = 0.0256 * Temp?rature - 0.1674
## NUTRIENTS
# Nitrification in the water column (from Chapelle 2000)
f_T_nit <- exp(kT*Temp)
Nitrif <- K_nit * f_T_nit * nh4
## Rates of change
#PHYTO <- PhytoS + PhytoL
dPhytoS <- + Growth_PS - Lysis_PS - grazing_PS - Exudation_PS
dPhytoL <- + Growth_PL - Lysis_PL - grazing_PL - Exudation_PL - Sed_PL
#ZOO <- ZS + ZL
dZS <- growth_ZS - Lysis_ZS - grazing_ZS
dZL <- growth_ZL - Lysis_ZL - egestion_ZL
#MACROPHYTES
dULVA <- growth_Ulva - death_Ulva * ulva
dQ_N_Ulva <- (NO3_uptake_Ulva + NH4_uptake_Ulva)*ulva -
death_Ulva*qN_ulva
dQ_P_Ulva <- PO4_uptake_Ulva*ulva - death_Ulva*qP_ulva
dGRACI <- growth_Graci - death_Graci*graci
dQ_N_Graci <- (NO3_uptake_Graci + NH4_uptake_Graci)*graci -
death_Graci*qN_graci
dQ_P_Graci <- PO4_uptake_Graci*graci - death_Graci*qP_graci
dZB <- + ZGR * zost - ZM * zost
dQ_N_Zost <- (LABS_NO3 + LABS_NH4)*zost - ZM * qN_zost
dQ_P_Zost <- (LABS_PO4) * zost - ZM * qP_zost
#ORGANIC MATTER
dl_PON <- (+ gammaP*egestion_ZL + epsilonP*Lys_Bio_N +
etaP*death_Ulva*ulva + omegaP*death_Graci*graci - Vsed_OM*l_pon -
K_lP*l_pon)
dl_POP <- (+ gammaP*(egestion_ZL/NP_zoo) + epsilonP*Lys_Bio_P +
phiP*death_Ulva*ulva + thetaP*death_Graci*graci - Vsed_OM*l_pop -
K_lP*l_pop)
dr_PON <- (+ ZM*qN_zost/z_fond - K_rP*r_pon - Vsed_OM*r_pon)
dr_POP <- (+ ZM*qP_zost/z_fond - K_rP*r_pop - Vsed_OM*r_pop)
dl_DON <- (+ Exudation_PS + Exudation_PL + K_lP*l_pon +
etaD*death_Ulva*ulva + omegaD*death_Graci*graci +
labi*FxDON*ratio_lagune + gammaD*egestion_ZL + epsilonD*Lys_Bio_N -
elys_l_DON*l_don)
dl_DOP <- (+ Exudation_PS/NP_phyto + Exudation_PL/NP_phyto +
K_lP*l_pop + phiD*death_Ulva*ulva + thetaD*death_Graci*graci +
labi*FxDOP*ratio_lagune + gammaD*(egestion_ZL/NP_zoo) +
epsilonD*Lys_Bio_P - elys_l_DOP*l_dop)
dr_DON <- (+ K_rP*r_pon + (1 - labi)*FxDON*ratio_lagune -
elys_r_DON*r_don)
dr_DOP <- (+ K_rP*r_pop + (1 - labi)*FxDOP*ratio_lagune -
elys_r_DOP*r_dop)
#SEDIMENT
dN_sed <- + Vsed_OM*(l_pon + r_pon) + Sed_PL - (FxNO3 + FxNH4 +
FxDON)*ratio_lagune - burial*sedN
dP_sed <- + Vsed_OM*(l_pop + r_pop) + Sed_PL/NP_phyto - (FxPO4 +
FxDOP)*ratio_lagune - burial*sedP
dN_burial <- + burial*sedN
dP_burial <- + burial*sedP
#NUTRIENTS
dNO3 <- (+ FxNO3*ratio_lagune + Nitrif - NO3_uptake_phyto -
NO3_uptake_Ulva*ulva - NO3_uptake_Graci*graci - LABS_NO3*ZB/z_zost)
dNH4 <- (+ FxNH4*ratio_lagune + elys_r_DON*r_don + elys_l_DON*l_don
- Nitrif - NH4_uptake_phyto - NH4_uptake_Ulva*ulva -
NH4_uptake_Graci*graci - LABS_NH4*(ZB/z_zost))
dPO4 <- (+ FxPO4*ratio_lagune + elys_r_DOP*r_dop + elys_l_DOP*l_dop
- PO4_uptake_phyto - PO4_uptake_Ulva*ulva - PO4_uptake_Graci*graci -
LABS_PO4*(ZB/z_zost))
# the output, packed as a list
res <- c(dPhytoS, dPhytoL, dZS, dZL, dULVA, dGRACI, dZB,
dQ_N_Ulva, dQ_N_Graci, dQ_N_Zost, dQ_P_Ulva, dQ_P_Graci,
dQ_P_Zost,
dl_PON, dl_POP, dr_PON, dr_POP, dl_DON, dl_DOP, dr_DON,
dr_DOP,
dNO3, dNH4, dPO4, dN_sed, dP_sed, dN_burial, dP_burial) #
the rate of change
accessory_all <- c(Chla=Chla, I0=I0, Temp=Temp,
I_phyto=I_phyto, f_I_PS=f_I_PS, f_I_PL=f_I_PL,
f_Temp_PS=f_Temp_PS, f_Temp_PL=f_Temp_PL,
NO3_uptake_phyto=NO3_uptake_phyto,
NH4_uptake_phyto=NH4_uptake_phyto, PO4_uptake_phyto=PO4_uptake_phyto,
f_N_PS=f_N_PS, f_P_PS=f_P_PS, f_N_PL=f_N_PL,
f_P_PL=f_P_PL, f_SN_PS=f_SN_PS, f_SN_PL=f_SN_PL,
mu_PS=mu_PS, mu_PL=mu_PL, Growth_PS=Growth_PS,
Growth_PL=Growth_PL, Lysis_PS=Lysis_PS, Lysis_PL=Lysis_PL,
Exudation_PS=Exudation_PS,
Exudation_PL=Exudation_PL, Sed_PL=Sed_PL,
gmax_ZS=gmax_ZS, gmax_ZL=gmax_ZL,
K_app_PS=K_app_PS, K_app_PL=K_app_PL, K_app_ZS=K_app_ZS,
grazing_PS=grazing_PS,
grazing_ZS_tot=grazing_ZS_tot, grazing_PL=grazing_PL,
grazing_ZS=grazing_ZS,
grazing_ZL_tot=grazing_ZL_tot, growth_ZS=growth_ZS, growth_ZL=growth_ZL,
egestion_ZL=egestion_ZL, Lysis_ZS=Lysis_ZS,
Lysis_ZL=Lysis_ZL,
I_Ulva=I_Ulva, f_I_Ulva=f_I_Ulva,
f_Temp_Ulva=f_Temp_Ulva,
f_QN_Ulva=f_QN_Ulva, f_QP_Ulva=f_QP_Ulva,
NO3_uptake_Ulva=NO3_uptake_Ulva,
NH4_uptake_Ulva=NH4_uptake_Ulva,
PO4_uptake_Ulva=PO4_uptake_Ulva, mu_Ulva=mu_Ulva,
growth_Ulva=growth_Ulva, death_Ulva=death_Ulva,
I_Graci=I_Graci, f_I_Graci=f_I_Graci,
f_Temp_Graci=f_Temp_Graci,
f_QN_Graci=f_QN_Graci, f_QP_Graci=f_QP_Graci,
NO3_uptake_Graci=NO3_uptake_Graci,
NH4_uptake_Graci=NH4_uptake_Graci,
PO4_uptake_Graci=PO4_uptake_Graci, mu_Graci=mu_Graci,
growth_Graci=growth_Graci, death_Graci=death_Graci,
Zost_Pmax=Zost_Pmax, Qcan=Qcan, f_can=f_can,
PAR_k1=PAR_k1, f_I_Zost=f_I_Zost, Zost_Ptot=Zost_Ptot,
f_QN_Zost=f_QN_Zost,f_QP_Zost=f_QP_Zost,
f_LN=f_LN, f_LP=f_LP, f_SN_zost=f_SN_zost, ZGR=ZGR, LABS_NO3=LABS_NO3,
LABS_NH4=LABS_NH4,
LABS_PO4=LABS_PO4, ZM_fT=ZM_fT, ZM=ZM,
K_lP=K_lP, K_rP=K_rP, Lys_Bio_N=Lys_Bio_N,
Lys_Bio_P=Lys_Bio_P,
V_lD=V_lD, V_rD=V_rD, elys_l_DON=elys_l_DON,
elys_r_DON=elys_r_DON,
elys_l_DOP=elys_l_DOP, elys_r_DOP=elys_r_DOP,
FxNO3=FxNO3, FxNH4=FxNH4, FxPO4=FxPO4,
FxDON=FxDON, FxDOP=FxDOP,
f_T_nit=f_T_nit, Nitrif=Nitrif)
list(res, accessory_all)
})
} # end of model
#-----------------------#
# the model parameters: #
#-----------------------#
parameters6 <- c(
# Water column
z_fond = 1.4, # average
depth of the lagoon - m
z_zost = 1.2, # average
depth of ZOstera canopy - m
z_phyto = 0.7, #
mid-average depth of the lagoon - m
# Phytoplanktno parameters
kI = 0.4 , # Phyto light
extinction coeff
#KT_phyto = 0.07, # Temperature
coefficient for phytoplankton (C°-1)
zeta_PS = 0.4,
zeta_PL = 0.4,
Temp_PS = 15,
Temp_PL = 12,
mu_max_PS = 0.07 , # Max growth rate PS
(d-1)
mu_max_PL = 0.05 , # Max growth rate PL
(d-1)
Iopt_PS = 200 , # Optimal light
intensity for PS (W.m-2.d-1)
Iopt_PL = 150 , # Optimal light
intensity for PL (W.m-2.d-1)
phi_PS = 1.5,
K_N_PS = 0.8, # PIco-nano
Half-saturation constante for nitrate (mmolN.m-3)
K_N_PL = 1.2,
K_P_PS = 0.02, #
Half-saturation constante of pico-nano for phosphate (mmolP.m-3)
K_P_PL = 0.1, #
Half-saturation constante of diatoms for phosphate (mmolP.m-3)
K_lys_PS = 0.000083 , # Lysis rate (d-1)
K_lys_PL = 0.00002 , #
K_e_PS = 0.0002 , # Exudation
rate of PS (d-1)
K_e_PL = 0.0002 , # Exudation
rate of PL (d-1)
vsed_PL = 0.00034 , #
Sedimentation rate of large phyto (m / d)
CHLaN_ratio = 1, #
Chlorophyll a to Nitrogen ratio in phytoplankton
NP_phyto = 16,
# Zooplankton parameters
g_prim_max_ZS = 0.0041 , # Grazing
rate of ZS at optimal temperature (d-1)
Temp_opt_ZS = 23, # ZS optimal
temperature (?C)
Temp_width_ZS = 12, # ZS optimal
temperature window (?C)
g_prim_max_ZL = 0.005 , # Grazing rate
of ZL at optimal temperature (d-1)
Temp_opt_ZL = 23, # ZL optimal
temperature (?C)
Temp_width_ZL = 8, # ZL
temperature window (?C)
K_PS = 1, #
Half-saturation constante on small phyto grazing (mmol.m-3)
K_PL = 0.8, #
K_ZS = 0.4, #
Y_ZS = 0.3, # Growth
efficiency for ZS (dimensionless)
Y_ZL = 0.33, #
E_ZL = 0.25, # Egestion
constante for ZL (dimensionless)
K_d_ZS = 0.000125 , # Martoality
rate for ZS (d-1)
K_d_ZL = 0.000042 , #
NP_zoo = 19.5, # Nitrogen
to Phosphorus ratio for zooplankton (average from Golz et al 2015)
# Macrophytes parameters
# Ulva
kI_Ulva = 0.4,
Iopt_Ulva = 180 , # Photosynthetic efficiency
(moles.m-2.d-1)
Temp_Ulva = 15 , # Temperature
reference for Ulva (?C)
zeta_Ulva = 0.2 , #
Temperature coefficient (?C-1)
Qmin_N_Ulva = 0.01 , # Min N quota
(mol N / gdw)
Qmin_P_Ulva = 0.001 , # Min P quota
(mol P / gdw)
mu_max_Ulva = 0.05 , # ulva max specific
growth rate (d-1)
Qmax_N_Ulva = 0.06 , # Max N quota
(mol N / gdw)
Qmax_P_Ulva = 0.0035 , # Max P quota
(mol P / gdw)
V_NO3_Ulva = 0.009 , # Max specific
NO3 uptake (d-1)
K_NO3_Ulva = 3.57 , # Half-saturation
constante for nitrate (mmol.m-3)
V_NH4_Ulva = 0.012 , # Max specific
NH4 uptake (d-1)
K_NH4_Ulva = 7.14 , # Half-saturation
constante for ammonium (mmol.m-3)
V_PO4_Ulva = 0.0003 , # Max specific
PO4 uptake (d-1)
K_PO4_Ulva = 5 , # Half-saturation
constante for phosphate (mmol.m-3)
MaxMort_Ulva = 0.0035 , # Max
death rate (d-1)
beta_Ulva = 0.84 , #Coeff of
intrinsec mortality
# Gracilaria
kI_Graci = 0.4,
Iopt_Graci = 150 , # Photosynthetic
efficiency (moles.m-2.d-1)
Temp_Graci = 13.5 , #
Temperature reference for Gracilaria (?C)
zeta_Graci = 0.2 , # Temperature
coefficient (?C-1)
Qmin_N_Graci = 0.01 , # Min N quota
(mol N / gdw)
Qmin_P_Graci = 0.001 , # Min P quota
(mol P / gdw)
mu_max_Graci = 0.03 , # Graci max
specific growth rate (d-1)
Qmax_N_Graci = 0.05 , # Max N quota
(mmol N / gdw)
Qmax_P_Graci = 0.0035 , # Max P quota
(mmol P / gdw)
V_NO3_Graci = 0.007 , # Max specific
NO3 uptake (d-1)
K_NO3_Graci = 27 , # Half-saturation constante
for nitrate (mmol.m-3)
V_NH4_Graci = 0.012 , # Max specific
NH4 uptake (d-1)
K_NH4_Graci = 51 , # Half-saturation constante
for ammonium (mmol.m-3)
V_PO4_Graci = 0.0004 , # Max specific
PO4 uptake (d-1)
K_PO4_Graci = 5 , # Half-saturation
constante for phosphate (mmol.m-3)
MaxMort_Graci = 0.0035 , # Max
death rate (d-1)
beta_Graci = 0.84 , #Coeff of
intrinsec mortality
# Zostera
theta_Pmax = 0.0265 , # Production
increasing rate with temperature
Pmax_0C = -0.0467 , # Theoritical
max production at 0?C
K1 = 0.4 , # Light
extinction coefficient (m-1)
K3 = 0.6 , # Part of
incoming light intercepted by the canopy (dimensionless)
Ikmax = 80 , # Max saturation
light intensity (moles.m-2.d-1) 80 *60*60*24/1000*217.5
Ikmin = 35 , # Min saturation
light intensity (moles.m-2.d-1) 35 *60*60*24/1000*217.5
LVmaxN = 0.007 , # Max Zostera
nitrogen uptake rate (d-1)
LVmaxP = 0.00054 , #Max
zostera P uptake (molP/molC/d) Zostera CP leaf: 230
KL_zost_NO3 = 9.2 , #
Half-saturation for Zostera uptake (mmol.m-3)
KL_zost_NH4 = 3 , #
Half-saturation for Zostera uptake (mmol.m-3)
KL_zost_PO4 = 12 , #
Half-saturation for Zostera uptake (mmol.m-3)
delta1_zost = 0.6 , # Limitation
coeff for zostera nitrate uptake
delta2_zost = 0.6 , # Limitation
coeff for zostera ammonium uptake
delta3_zost = 0.6 , # Limitation
coeff for zostera phosphate uptake
theta_LM = 1.1 , # Zostera
mortality increasing rate with temperature
#K4 = 1.2 , # Wind effect
attenuation with depth (m-1)
#LMR_V = 0.08 , # Zostera
leaf slooghing coeff due to wind (d-1)
LMR_20C = 0.025 , # Max zsotera
(leaf) mortality rate at 20?C (d-1)
epsilon1_zost = 0.4 , # Zostera
(leaf) nitrogen content limitation coeff (dimensionless)
epsilon2_zost = 0.4 , # Zostera
(leaf) phosphorus content limitation coeff (dimensionless)
LN_max = 0.07 , # Max N quota
for Zostera (mol:mol)
LN_min = 0.03 , # Min N quota
for Zostera
LP_max = 0.0035 , # Max P quota
LP_min = 0.0015 , # Min P quota
# Organic matter parameters
K_lP_max = 0.0002 , # Max.
hydrolysis rate of labile POM (d-1)
K_rP_max = 0.00001 , # Max. hydrolysis
rate of refractory POM (d-1)
Temp_opt_OM = 30 , # Optimal
temperature for POM hydrolysis in ?C
Temp_width_OM = 18 , # Temperature
window for POM hydrolysis in ?C
Vmax_lD = 0.032 , # Max. rate of
labile DON hydrolysis (d-1)
Vmax_rD = 0.010 , # Max. rate of
refractcory DON hydrolysis (d-1)
K_l = 1.1,
K_r = 11.1,
# Zaldivar gives a ratio of 10 between labile and refractory
half-saturation cste which is used here to differenaciate K_l & K_r DOM
gammaP = 0.7 , # Fraction of
Zooplankton egestion going to labile POM
gammaD = 0.3 , # Fraction of
Zooplankton egestion going to labile DOM
epsilonP = 0.8 , # Fraction of
biological lysis going to labile POM
epsilonD = 0.2 , # Fraction of
biological lysis going to labile DOM
etaP = 0.2 , # Fraction of
Ulva N to labile POM
etaD = 0.3 , # Fraction of
Ulva N to labile DOM
phiP = 0.1 , # Fraction of
Ulva P to labile POM
phiD = 0.4 , # Fraction of
Ulva P to labile DOM
omegaP = 0.6 , # Fraction of
Gracilaria Nitrogen going to labile POM
omegaD = 0.3 , # Fraction of
Gracilaria Nitrogen going to labile DOM
thetaP = 0.2 , # Fraction of
Gracilaria Phosphorus going to labile POM
thetaD = 0.1 , # Fraction of
Gracilaria Phosphorus going to labile DOM
Vsed_OM = 0.000542 , # Sedimentation
rate of particulate organic matter (d-1)
# Sediment
burial = 0.00005, # Fraction
of organic matter buried in deeper sediment
labi = 0.2, # Fraction of
labile DOM released from the sediment
ratio_lagune = 0.6, # ratio
surface/volume (m2/m3) to transforme Fx in mmol/m3/d
# Nutrient
kT = 0.07 , # Temperature
increasing rate (?C-1)
K_nit = 0.00035 # Nitrification
rate at 0?C (d-1)
)
#-------------------------#
# the initial conditions: #
#-------------------------#
state6 <- c(
PhytoS = 1.5,
PhytoL = 1.5,
ZS = 0.1,
ZL = 0.1,
ULVA = 0.1,
GRACI = 0.1,
ZB = 0.1,
Q_N_Ulva = 0.02,
Q_N_Graci = 0.02,
Q_N_Zost = 0.04,
Q_P_Ulva = 0.003,
Q_P_Graci= 0.003,
Q_P_Zost = 0.002,
l_PON = 0.01,
l_POP = 0.01,
r_PON = 0.01,
r_POP = 0.01,
l_DON = 0.01,
l_DOP = 0.01,
r_DON = 0.01,
r_DOP = 0.01,
NO3 = 15, # mmol/m3 (µmol/L)
NH4 = 5, # mmol/m3 (?mol/L)
PO4 = 1, # mmol/m3 (?mol/L)
N_sed = 0.01,
P_sed = 0.01,
N_burial = 0.01,
P_burial = 0.01
)
--
Romain Pete
UMR MARBEC
Ifremer - Laboratoire LER/LR
Avenue Jean Monnet
CS 30171 - 34203 Sète cedex France
Tél : (+33)4 99 57 32 93 - Fax : 04 99 57 32 96
More information about the R-sig-dynamic-models
mailing list