[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