# the extended diagenetic model for shallow water #============================================== # Authigenic carbonate and phosphorus model 0.1 #============================================== # A multicomponent diagenetic model for the formation of authigenic carbonate and phosphorus cycle # pH calculated using Equation (7) in Hofmann et al., 2009 # ===================load required packages================= require(ReacTran) require(seacarb) require(deSolve) require(marelac) require(NORMT3) # ------------------------------------------------------ # # --------------------- Parameters --------------------- # # ------------------------------------------------------ # # ==================== Bottom water properties ========================================== BW_Temp <- 12 # bottom water temperature deg C BW_Sal <- 28.4 # bottom water Salinity BW_Pr <- 2 # bottom water pressure (bar) D_sw <- rho(BW_Sal,BW_Temp,BW_Pr)[[1]] # bottom water density at {T,S,P} (kg/m^3) BW_pH <- 7.64 # bottom water pH BW_H <- 10^(-BW_pH) # mmol/cm^3 BW_SumCO2 <- 2/1e3 # bottom water DIC, mmol/cm^3 BW_SumH2S <- 0 # bottom water H2S, mmol/cm^3 BW_Ca <- 8.6/1e3 # bottom water Ca, mmol/cm^3 BW_Mg <- 46/1e3 # bottom water Mg, mmol/cm^3 BW_Na <- 380.56/1e3 # bottom water Na, mmol/cm^3 BW_F <- 0.07/1e3 # bottom water F, mmol/cm^3 BW_SumSO4 <- 22/1e3 # bottom water sulfate, mmol/cm^3 BW_O2 <- 0.15/1e3 # bottom water oxygen, mmol/cm^3 BW_NO3 <- 11.8/1e6 # bottom water NO3, mmol/cm^3 BW_MnII <- 0/1e6 # bottom water Mn, mmol/cm^3 BW_FeII <- 0/1e6 # bottom water Fe, mmol/cm^3 BW_SumNH4 <- 0/1e6 # bottom water ammonia, mmol/cm^3 BW_CH4 <- 0 # bottom water methane, mmol/cm^3 BW_SumH2PO4 <- 1/1e6 # bottom water phosphate, mmol/cm^3 test # ==================== The composition of organic matter and Fe oxyhydroxides====================== x <- 10/106 # moler ratio of N:C in highly reactive organic matter xx <- 10/106 # moler ratio of N:C in the other organic matter y <- 1.3/106 # moler ratio of P:C in highly reactive organic matter yy <- 0.27/106 # moler ratio of P:C in the other organic matter r0 <- 0.24 # moler ratio of P:Fe for Fe bound P # ==================== The 12 G model from continuum theory ====================== a=0.15 # Average lifetime of more reactive orgC, yr v=0.12 # Shape of orgC distribution integrand <- function(x) {(x^(v-1))*exp(-x)} Total = integrate(integrand, lower = 0, upper = Inf, rel.tol =1e-8)$value k=1:12 f=1:12 k[1]= 1 k[12]=1e-10 for (i in 2:11) {k[i]=10^(-i+1+0.5)} for (i in 2:11) { f[i]=(integrate(integrand, lower = 0, upper = (a*10^(-i+2)), rel.tol =1e-8)$value-integrate(integrand, lower = 0, upper = (a*10^(-i+1)), rel.tol =1e-8)$value)/Total } f[1] = (Total-integrate(integrand, lower = 0, upper = (a*1), rel.tol =1e-8)$value)/Total f[12]= integrate(integrand, lower = 0, upper = (a*1e-10), rel.tol =1e-8)$value/Total # ==================== Solid flux at the sediment-seawater interface ==================================== J.orgC <- 0.3 # mmol/(cm*cm*yr) organic C J.orgC1 <- J.orgC*f[1] # mmol/(cm*cm*yr) rapidly oxidized OC J.orgC2 <- J.orgC*f[2] # mmol/(cm*cm*yr) rapidly oxidized OC J.orgC3 <- J.orgC*f[3] # mmol/(cm*cm*yr) less reactive OC J.orgC4 <- J.orgC*f[4] # mmol/(cm*cm*yr) less reactive OC J.orgC5 <- J.orgC*f[5] # mmol/(cm*cm*yr) less reactive OC J.orgC6 <- J.orgC*f[6] # mmol/(cm*cm*yr) less reactive OC J.orgC7 <- J.orgC*f[7] # mmol/(cm*cm*yr) less reactive OC J.orgC8 <- J.orgC*f[8] # mmol/(cm*cm*yr) less reactive OC J.orgC9 <- J.orgC*f[9] # mmol/(cm*cm*yr) less reactive OC J.orgC10 <- J.orgC*f[10] # mmol/(cm*cm*yr) less reactive OC J.orgC11 <- J.orgC*f[11] # mmol/(cm*cm*yr) less reactive OC J.orgC12 <- J.orgC*f[12] # mmol/(cm*cm*yr) less reactive OC J.orgCc <- J.orgC7+J.orgC8+J.orgC9+J.orgC10+J.orgC11+J.orgC12 J.FeOH3a <- 0.213/100*2.5/56*1000*0.2*0.36 # mmol/(cm*cm*yr) amorphous Fe J.FeOH3b <- 0 # mmol/(cm*cm*yr) more crystalline Fe J.MnO2a <- 1000/1e6*2.5/55*1000*0.2*0.36 # mmol/(cm*cm*yr) amorphous Mn J.MnO2b <- 0 # mmol/(cm*cm*yr) more crystalline Mn J.irPa <- J.FeOH3a*r0 # mmol/(cm*cm*yr) P bound by amorphous Fe J.irPb <- J.FeOH3b*r0 # mmol/(cm*cm*yr) P bound by more crystalline Fe J.S0 <- 0 # mmol/(cm*cm*yr) S0 J.FeS <- 0 # mmol/(cm*cm*yr) FeS J.FeS2 <- 0 # mmol/(cm*cm*yr) FeS2 J.CFA <- 0.1*436/1e6*2.5/31*1000/4.8*0.2*0.36 # mmol/(cm*cm*yr) CFA J.calc <- 2.88/100*2.5/12*1000*0.2*0.36 # mmol/(cm*cm*yr) Calcite J.arag <- 0 # mmol/(cm*cm*yr) aragonite J.biot <- 0.004 # mmol/(cm*cm*yr) Biotite J.magn <- 0.062/100*2.5/56*1000*0.2*0.36 # mmol/(cm*cm*yr) magnetite J.SurfFe <- 0 # mmol/(cm*cm*yr) adsorbed iron # ==================== Equilibrium constants within porewater ==================================================== K1CO2 <- K1(BW_Sal, BW_Temp)*(D_sw/1000) # mmol/(cm*cm*cm) K2CO2 <- K2(BW_Sal, BW_Temp)*(D_sw/1000) # mmol/(cm*cm*cm) K1P <- K1p(BW_Sal, BW_Temp)*(D_sw/1000) # mmol/(cm*cm*cm) K2P <- K2p(BW_Sal, BW_Temp)*(D_sw/1000) # mmol/(cm*cm*cm) K3P <- K3p(BW_Sal, BW_Temp)*(D_sw/1000) # mmol/(cm*cm*cm) KS <- Ks(BW_Sal, BW_Temp)*(D_sw/1000) # mmol/(cm*cm*cm) KN <- Kn(BW_Sal, BW_Temp)*(D_sw/1000) # mmol/(cm*cm*cm) KW <- Kw(BW_Sal, BW_Temp)*(D_sw/1000)*(D_sw/1000) # mmol*mmol/cm^6 KHS <- Khs(BW_Sal,BW_Temp)[[1]]*(D_sw/1000) # mmol/(cm*cm*cm) Kspa <- Kspa(BW_Sal, BW_Temp)[[1]]*(D_sw/1000)*(D_sw/1000) # mmol*mmol/cm^6 Kspc <- Kspc(BW_Sal, BW_Temp)[[1]]*(D_sw/1000)*(D_sw/1000) # mmol*mmol/cm^6 TK <- 273.15 + BW_Temp Kspr <- 10^(-9) # mmol*mmol/cm^6 KspCFA <- 10^(-83.231) # mmol*mmol/cm^6 # ===================== Derived concentrations ========================================== BW_CO2 <- BW_H*BW_H/(BW_H*K1CO2 + BW_H*BW_H + K1CO2*K2CO2)*BW_SumCO2 BW_HCO3 <- BW_H*K1CO2/(BW_H*K1CO2 + BW_H*BW_H + K1CO2*K2CO2)*BW_SumCO2 BW_CO3 <- K1CO2*K2CO2/(BW_H*K1CO2 + BW_H*BW_H + K1CO2*K2CO2)*BW_SumCO2 BW_HS <- BW_SumH2S/(1+BW_H/KHS) BW_H2S <- BW_SumH2S - BW_HS BW_PO4 <- 1/(1+BW_H/K3P+BW_H*BW_H/K2P/K3P+BW_H*BW_H*BW_H/K1P/K2P/K3P)*BW_SumH2PO4 BW_HPO4 <- (BW_H/K3P)/(1+BW_H/K3P+BW_H*BW_H/K2P/K3P+BW_H*BW_H*BW_H/K1P/K2P/K3P)*BW_SumH2PO4 BW_H2PO4 <- (BW_H*BW_H/K2P/K3P)/(1+BW_H/K3P+BW_H*BW_H/K2P/K3P+BW_H*BW_H*BW_H/K1P/K2P/K3P)*BW_SumH2PO4 BW_H3PO4 <- (BW_H*BW_H*BW_H/K1P/K2P/K3P)/(1+BW_H/K3P+BW_H*BW_H/K2P/K3P+BW_H*BW_H*BW_H/K1P/K2P/K3P)*BW_SumH2PO4 BW_SO4 <- BW_SumSO4/(1+BW_H/KS) BW_HSO4 <- BW_SumSO4-BW_SO4 BW_NH3 <- BW_SumNH4/(1+BW_H/KN) BW_NH4 <- BW_SumNH4-BW_NH3 # ===================== activity coefficient ========================================== rCa <- 0.21866266 rNa <- 0.65687384 rMg <- 0.2819323 rPO4 <- 0.000037 rCO3 <- 0.028861 rF <- 0.31540788 # ===================== Caculation of diffusion coeffiecients ================================================================================================== select <- c("SO4", "Ca", "O2", "HCO3", "CO2", "CO3", "HS","H2S", "H", "NO3", "NH4", "H3PO4", "H2PO4", "HPO4", "PO4", "Mn", "Fe", "CH4", "Na", "Mg", "F","HSO4","NH3") diffsal <- diffcoeff(S = BW_Sal, t= BW_Temp, P=BW_Pr, species = select) # convert diffusion coefficients into cm*cm/yr diffsal <- diffsal *100*100*31556926 # extract diffusion coefficients from diffsal DSO4 <- diffsal[1,1] DCa <- diffsal[1,2] DO2 <- diffsal[1,3] DHCO3 <- diffsal[1,4] DCO2 <- diffsal[1,5] DCO3 <- diffsal[1,6] DHS <- diffsal[1,7] DH2S <- diffsal[1,8] DH <- diffsal[1,9] DNO3 <- diffsal[1,10] DNH4 <- diffsal[1,11] DH3PO4 <- diffsal[1,12] DH2PO4 <- diffsal[1,13] DHPO4 <- diffsal[1,14] DPO4 <- diffsal[1,15] DMn <- diffsal[1,16] DFe <- diffsal[1,17] DCH4 <- diffsal[1,18] DNa <- diffsal[1,19] DMg <- diffsal[1,20] DF <- diffsal[1,21] DHSO4 <- diffsal[1,22] DNH3 <- diffsal[1,23] #=======================parameters for bioturbation============================= Db <- 1e-7*31556926 # cm*cm/yr Biodiffusion coefficient at surface xbt <- 3 # cm Biodiffusion attenuation coefficient a0 <- 100 # yr-1 Bioirrigation coefficient at surface xbi <- 0.8 # cm Bioirrigation attenuation coefficient #====================== Reaction parameters =============================================== Ks_O2 <- 0.02/1e3 # mmol/cm^3 limitting concentration of O2 Ks_NO3 <- 0.004/1e3 # mmol/cm^3 limitting concentration of NO3 Ks_MnO2 <- 0.032*2.5 # mmol/cm^3 limitting concentration of MnO2 Ks_FeOH3 <- 0.065*2.5 # mmol/cm^3 limitting concentration of FeOH3 Ks_SO4 <- 1.6/1e3 # mmol/cm^3 limitting concentration of SO4 a_SO4 <- 0.2 # attenuation factor for SO4 reduction KPO4 <- 10/1e6 k7 <- 1e7 # cm^3/(mmol*yr) rate constant for E7 k8 <- 6.94*1e5 # cm^3/(mmol*yr) rate constant for E8 k9 <- 1.4*1e8 # cm^3/(mmol*yr) rate constant for E9 k10 <- 3*1e5 # cm^3/(mmol*yr) rate constant for E10 k11 <- 1.89*1e4 # cm^3/(mmol*yr) rate constant for E11 k12 <- 1.6*1e5 # cm^3/(mmol*yr) rate constant for E12 k13 <- 1e10 # cm^3/(mmol*yr) rate constant for E13 k14 <- 3*1e6 # cm^3/(mmol*yr) rate constant for E14 k15 <- 2*1e4 # cm^3/(mmol*yr) rate constant for E15 k16 <- 8*1e3 # cm^3/(mmol*yr) rate constant for E16 k17 <- 1.48e6 # cm^3/(mmol*yr) rate constant for E17 k18 <- 1e4 # cm^3/(mmol*yr) rate constant for E18 k19 <- 3.16 # /yr rate constant for E19 k20 <- 7.26*1e3 # cm^3/(mmol*yr) rate constant for E20 k21 <- 0.57 # /yr rate constant for E21 k22 <- 1.7 # /yr rate constant for E22 k28 <- 3.25*1e3 # cm^3/(mmol*yr) rate constant for E28 k.biotite <- 3e-4 # yr-1 apparent rate constants for biotite dissolution k.magn <- 0.95e-5*570 # mM-0.5yr-1 apparent rate constants for magnetite dissolution kFe <- 500 # adsorption coefficient of Fe k32 <- 5e6 # cm^3/(mmol*yr) rate constant for E32 #====================== Carbonate precipitation kinetics ================================================================== Kpmaxa <- 3e-6 # mmol/cm*cm*cm/yr reaction constant for aragonite precipitation Kdissa <- 0.5 # /yr kinetic constant of aragonite dissolution Kpmaxc <- 3e-7 # mmol/cm*cm*cm/yr reaction constant for calcite precipitation Kdissc <- 0.5 # /yr kinetic constant of calcite dissolution Kpmaxr <- 3e-6 # mmol/cm*cm*cm/yr reaction constant for rhodochrosite precipitation Kdissr <- 0.25 # /yr kinetic constant of rhodochrosite dissolution KpmaxCFA <- 2.7e-8 # mmol/cm*cm*cm/yr reaction constant for CFA precipitation # ----------------------------------------------------------------------------- # # --------------------- setup model domain and properties --------------------- # # ----------------------------------------------------------------------------- # L <- 300 # size of sediment column (cm) N <- 20 # number of grid layers dbl <- 0.1 # thickness of diffusive boundary layer (cm) por.0 <- 0.64 # porosity at sediment seawater interface por.inf <- 0.64 # porosity at infinite depth por_coeff <- 5 # attenuation coefficient for porosity profile cm rho_sed <- 2.5 # sediment density (g/cm^3) MAR <- rho_sed*0.36*0.2 # g/cm3/yr v.0 <- MAR/(1-por.0)/rho_sed # solid-phase advection velocity at SWI (cm/y) v.inf <- NULL # solid-phase advection velocity at infinite depth (cm/y) grid <- setup.grid.1D(x.up = 0, L = L, N = N, dx.1 = 0.3) # create finite element grid por.grid <- setup.prop.1D(func = p.exp, grid = grid, y.0 = por.0, y.inf = por.inf, x.att = por_coeff) # setup porosity profile svf.grid <- setup.prop.1D(func = p.exp, grid = grid, y.0 = 1-por.0, y.inf = 1-por.inf, x.att = por_coeff) # setup solid volume fraction tort <- 1 - 2*log(por.grid$int) # calculate tortuosity from porosity adv <- setup.compaction.1D(v.0 = v.0, v.inf = v.inf, por.0 = por.0, por.inf = por.inf, por.grid = por.grid) # setup advection grid u.grid <- adv$u # porewater advection velocities (cm/y) v.grid <- adv$v # solid-phase advection velocities (cm/y) f.grid <- (1-por.grid$mid)/por.grid$mid #====================== setup of bioturbation grid ============================================== Dbfunc <- function (x) return(Db*exp(-(x/xbt)^2)) Db.grid <- setup.prop.1D(func = Dbfunc, grid = grid)$int # biodiffusion afunc <- function (x) return(a0*exp(-x/xbi)) a.grid <- setup.prop.1D(func = afunc, grid = grid)$mid # bioirrigation # ------------------------------------------------------------- # # --------------------- model formulation --------------------- # # ------------------------------------------------------------- # Pmodel <- function (t = 0, Conc, pars = NULL) { #------------------------ pass in intial concentrations ------------------------ orgC1 <- Conc[1:N] orgC2 <- Conc[(N+1):(2*N)] orgC3 <- Conc[(2*N+1):(3*N)] orgC4 <- Conc[(3*N+1):(4*N)] orgC5 <- Conc[(4*N+1):(5*N)] orgC6 <- Conc[(5*N+1):(6*N)] CFA <- Conc[(6*N+1):(7*N)] O2 <- Conc[(7*N+1):(8*N)] NO3 <- Conc[(8*N+1):(9*N)] MnO2a <- Conc[(9*N+1):(10*N)] FeOH3a <- Conc[(10*N+1):(11*N)] SO4 <- Conc[(11*N+1):(12*N)] CH4 <- Conc[(12*N+1):(13*N)] SumNH4 <- Conc[(13*N+1):(14*N)] MnII <- Conc[(14*N+1):(15*N)] FeII <- Conc[(15*N+1):(16*N)] FeS <- Conc[(16*N+1):(17*N)] FeS2 <- Conc[(17*N+1):(18*N)] SumH2S <- Conc[(18*N+1):(19*N)] MnO2b <- Conc[(19*N+1):(20*N)] SumH2PO4 <- Conc[(20*N+1):(21*N)] S0 <- Conc[(21*N+1):(22*N)] FeOH3b <- Conc[(22*N+1):(23*N)] F <- Conc[(23*N+1):(24*N)] SumCO2 <- Conc[(24*N+1):(25*N)] H <- Conc[(25*N+1):(26*N)] Ca <- Conc[(26*N+1):(27*N)] arag <- Conc[(27*N+1):(28*N)] calc <- Conc[(28*N+1):(29*N)] rhod <- Conc[(29*N+1):(30*N)] Mg <- Conc[(30*N+1):(31*N)] magn <- Conc[(31*N+1):(32*N)] irPa <- Conc[(32*N+1):(33*N)] irPb <- Conc[(33*N+1):(34*N)] #------------------------ Dissociation terms ----------------------------- HS <- SumH2S/(1+H/KHS) H2S <- SumH2S - HS PO4 <- 1/(1+H/K3P+H*H/K2P/K3P+H*H*H/K1P/K2P/K3P)*SumH2PO4 HPO4 <- (H/K3P)/(1+H/K3P+H*H/K2P/K3P+H*H*H/K1P/K2P/K3P)*SumH2PO4 H2PO4 <- (H*H/K2P/K3P)/(1+H/K3P+H*H/K2P/K3P+H*H*H/K1P/K2P/K3P)*SumH2PO4 H3PO4 <- (H*H*H/K1P/K2P/K3P)/(1+H/K3P+H*H/K2P/K3P+H*H*H/K1P/K2P/K3P)*SumH2PO4 NH3 <- SumNH4/(1+H/KN) NH4 <- SumNH4-NH3 CO2 <- H*H/(H*K1CO2 + H*H + K1CO2*K2CO2)*SumCO2 HCO3 <- H*K1CO2/(H*K1CO2 + H*H + K1CO2*K2CO2)*SumCO2 CO3 <- K1CO2*K2CO2/(H*K1CO2 + H*H + K1CO2*K2CO2)*SumCO2 SurfFe <- kFe*FeII/f.grid #------------------------ differential equations to keep track of pH -------------------------- deltaTP.deltaSumCO2 <- 2*H*H/(H*K1CO2 + H*H + K1CO2*K2CO2) + H*K1CO2/(H*K1CO2 + H*H + K1CO2*K2CO2) deltaTP.deltaSumH2S <- H/(H+KHS) deltaTP.deltaSumH2PO4 <- 2*K1P*H*H/(K1P*K2P*K3P+K1P*K2P*H+K1P*H*H+H*H*H)+K1P*K2P*H/(K1P*K2P*K3P+K1P*K2P*H+K1P*H*H+H*H*H) deltaTP.deltaSumNH4 <- H/(H+KN) deltaCO2.deltaH <- (2*H*(H*K1CO2 + H*H + K1CO2*K2CO2)-(2*H + K1CO2)*H*H)/(H*K1CO2 + H*H + K1CO2*K2CO2)^2 * SumCO2 deltaHCO3.deltaH <- (K1CO2*(H*K1CO2 + H*H + K1CO2*K2CO2)-(2*H + K1CO2)*H*K1CO2)/(H*K1CO2 + H*H + K1CO2*K2CO2)^2 * SumCO2 deltaH2PO4.deltaH <- (2*K1P*H*(K1P*K2P*K3P+K1P*K2P*H+K1P*H*H+H*H*H)-(3*H*H+2*H*K1P+K1P*K2P)*K1P*H*H)/(K1P*K2P*K3P+K1P*K2P*H+K1P*H*H+H*H*H)^2*SumH2PO4 deltaHPO4.deltaH <- (K1P*K2P*(K1P*K2P*K3P+K1P*K2P*H+K1P*H*H+H*H*H)-(3*H*H+2*H*K1P+K1P*K2P)*H*K1P*K2P)/(K1P*K2P*K3P+K1P*K2P*H+K1P*H*H+H*H*H)^2*SumH2PO4 deltaNH4.deltaH <- KN/(H+KN)^2*SumNH4 deltaH2S.deltaH <- (KHS /(H*H + 2*H*KHS + KHS*KHS)) * SumH2S deltaTP.deltaH <- 2*deltaCO2.deltaH + deltaHCO3.deltaH + 2*deltaH2PO4.deltaH + deltaHPO4.deltaH + deltaH2S.deltaH + deltaNH4.deltaH + 1 #------------------------ molecule diffusion, biodiffusion, irrigation and advection functions ----------------------------------------------- tranorgC1 <- tran.1D(C = orgC1, flux.up = J.orgC1, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranorgC2 <- tran.1D(C = orgC2, flux.up = J.orgC2, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranorgC3 <- tran.1D(C = orgC3, flux.up = J.orgC3, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranorgC4 <- tran.1D(C = orgC4, flux.up = J.orgC4, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranorgC5 <- tran.1D(C = orgC5, flux.up = J.orgC5, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranorgC6 <- tran.1D(C = orgC6, flux.up = J.orgC6, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranCFA <- tran.1D(C = CFA, flux.up = J.CFA, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranMnO2a <- tran.1D(C = MnO2a, flux.up = J.MnO2a, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranFeOH3a <- tran.1D(C = FeOH3a, flux.up = J.FeOH3a, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranFeS <- tran.1D(C = FeS, flux.up = J.FeS, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranFeS2 <- tran.1D(C = FeS2, flux.up = J.FeS2, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranMnO2b <- tran.1D(C = MnO2b, flux.up = J.MnO2b, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranS0 <- tran.1D(C = S0, flux.up = J.S0, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranFeOH3b <- tran.1D(C = FeOH3b, flux.up = J.FeOH3b, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranarag <- tran.1D(C = arag, flux.up = J.arag, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC trancalc <- tran.1D(C = calc, flux.up = J.calc, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranrhod <- tran.1D(C = rhod, flux.up = 0, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranmagn <- tran.1D(C = magn, flux.up = J.magn, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranSurfFe <- tran.1D(C = SurfFe, flux.up = J.SurfFe, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranirPa <- tran.1D(C = irPa, flux.up = J.irPa, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranirPb <- tran.1D(C = irPb, flux.up = J.irPb, D = Db.grid, v = v.grid, VF = svf.grid, dx = grid)$dC tranO2 <- tran.1D(C = O2, C.up = BW_O2, D = DO2/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(O2-BW_O2) tranNO3 <- tran.1D(C = NO3, C.up = BW_NO3, D = DNO3/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(NO3-BW_NO3) tranSO4 <- tran.1D(C = SO4, C.up = BW_SO4, D = DSO4/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(SO4-BW_SO4) tranCH4 <- tran.1D(C = CH4, C.up = BW_CH4, D = DCH4/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(CH4-BW_CH4) tranNH3 <- tran.1D(C = NH3, C.up = BW_NH3, D = DNH3/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(NH3-BW_NH3) tranNH4 <- tran.1D(C = NH4, C.up = BW_NH4, D = DNH4/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(NH4-BW_NH4) tranMnII <- tran.1D(C = MnII, C.up = BW_MnII, D = DMn/tort, v = u.grid, VF = por.grid, dx = grid)$dC - 0.2*a.grid*(MnII-BW_MnII) tranFeII <- tran.1D(C = FeII, C.up = BW_FeII, D = DFe/tort, v = u.grid, VF = por.grid, dx = grid)$dC tranHS <- tran.1D(C = HS, C.up = BW_HS, D = DHS/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(HS-BW_HS) tranH2S <- tran.1D(C = H2S, C.up = BW_H2S, D = DH2S/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(H2S-BW_H2S) tranPO4 <- tran.1D(C = PO4, C.up = BW_PO4, D = DPO4/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(PO4-BW_PO4) tranHPO4 <- tran.1D(C = HPO4, C.up = BW_HPO4, D = DHPO4/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(HPO4-BW_HPO4) tranH2PO4 <- tran.1D(C = H2PO4, C.up = BW_H2PO4, D = DH2PO4/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(H2PO4-BW_H2PO4) tranH3PO4 <- tran.1D(C = H3PO4, C.up = BW_H3PO4, D = DH3PO4/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(H3PO4-BW_H3PO4) tranF <- tran.1D(C = F, C.up = BW_F, D = DF/tort, v = u.grid, VF = por.grid, dx = grid)$dC - 2*a.grid*(F-BW_F) tranHCO3 <- tran.1D(C = HCO3, C.up = BW_HCO3, D = DHCO3/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(HCO3-BW_HCO3) tranCO2 <- tran.1D(C = CO2, C.up = BW_CO2, D = DCO2/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(CO2-BW_CO2) tranCO3 <- tran.1D(C = CO3, C.up = BW_CO3, D = DCO3/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(CO3-BW_CO3) tranH <- tran.1D(C = H, C.up = BW_H, D = DH/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(H-BW_H) tranCa <- tran.1D(C = Ca, C.up = BW_Ca, D = DCa/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(Ca-BW_Ca) tranMg <- tran.1D(C = Mg, C.up = BW_Mg, D = DMg/tort, v = u.grid, VF = por.grid, dx = grid)$dC - a.grid*(Mg-BW_Mg) tranSumNH4 <- tranNH3 + tranNH4 tranSumH2S <- tranHS + tranH2S tranSumH2PO4 <- tranPO4 + tranHPO4 + tranH2PO4 + tranH3PO4 tranSumCO2 <- tranCO3 + tranHCO3 + tranCO2 tranTP <- 2*tranCO2 + tranHCO3 + 2*tranH2PO4 + tranHPO4 + tranH2S + tranNH4 + tranH #------------------------ reaction terms -------------------------------------------------------------------- reacOC1i <- k[1]*orgC1 reacOC2i <- k[2]*orgC2 reacOC3i <- k[3]*orgC3 reacOC4i <- k[4]*orgC4 reacOC5i <- k[5]*orgC5 reacOC6i <- k[6]*orgC6 reacorgC11 <- reacOC1i*O2/(Ks_O2+O2) # R1 reacorgC21 <- reacOC2i*O2/(Ks_O2+O2) # R1 reacorgC31 <- reacOC3i*O2/(Ks_O2+O2) # R1 reacorgC41 <- reacOC4i*O2/(Ks_O2+O2) # R1 reacorgC51 <- reacOC5i*O2/(Ks_O2+O2) # R1 reacorgC61 <- reacOC6i*O2/(Ks_O2+O2) # R1 reacorgC12 <- reacOC1i*(Ks_O2/(Ks_O2+O2))*(NO3/(Ks_NO3+NO3)) # R2 reacorgC22 <- reacOC2i*(Ks_O2/(Ks_O2+O2))*(NO3/(Ks_NO3+NO3)) # R2 reacorgC32 <- reacOC3i*(Ks_O2/(Ks_O2+O2))*(NO3/(Ks_NO3+NO3)) # R2 reacorgC42 <- reacOC4i*(Ks_O2/(Ks_O2+O2))*(NO3/(Ks_NO3+NO3)) # R2 reacorgC52 <- reacOC5i*(Ks_O2/(Ks_O2+O2))*(NO3/(Ks_NO3+NO3)) # R2 reacorgC62 <- reacOC6i*(Ks_O2/(Ks_O2+O2))*(NO3/(Ks_NO3+NO3)) # R2 reacorgC13 <- reacOC1i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(MnO2a/(Ks_MnO2+MnO2a)) # R3 reacorgC23 <- reacOC2i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(MnO2a/(Ks_MnO2+MnO2a)) # R3 reacorgC33 <- reacOC3i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(MnO2a/(Ks_MnO2+MnO2a)) # R3 reacorgC43 <- reacOC4i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(MnO2a/(Ks_MnO2+MnO2a)) # R3 reacorgC53 <- reacOC5i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(MnO2a/(Ks_MnO2+MnO2a)) # R3 reacorgC63 <- reacOC6i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(MnO2a/(Ks_MnO2+MnO2a)) # R3 reacorgC14 <- reacOC1i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(FeOH3a/(Ks_FeOH3+FeOH3a)) # R4 reacorgC24 <- reacOC2i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(FeOH3a/(Ks_FeOH3+FeOH3a)) # R4 reacorgC34 <- reacOC3i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(FeOH3a/(Ks_FeOH3+FeOH3a)) # R4 reacorgC44 <- reacOC4i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(FeOH3a/(Ks_FeOH3+FeOH3a)) # R4 reacorgC54 <- reacOC5i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(FeOH3a/(Ks_FeOH3+FeOH3a)) # R4 reacorgC64 <- reacOC6i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(FeOH3a/(Ks_FeOH3+FeOH3a)) # R4 reacorgC15 <- a_SO4*reacOC1i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (SO4/(Ks_SO4+SO4)) # R5 reacorgC25 <- a_SO4*reacOC2i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (SO4/(Ks_SO4+SO4)) # R5 reacorgC35 <- a_SO4*reacOC3i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (SO4/(Ks_SO4+SO4)) # R5 reacorgC45 <- a_SO4*reacOC4i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (SO4/(Ks_SO4+SO4)) # R5 reacorgC55 <- a_SO4*reacOC5i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (SO4/(Ks_SO4+SO4)) # R5 reacorgC65 <- a_SO4*reacOC6i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (SO4/(Ks_SO4+SO4)) # R5 reacorgC16 <- a_SO4*reacOC1i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (Ks_SO4/(Ks_SO4+SO4)) # R6 reacorgC26 <- a_SO4*reacOC2i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (Ks_SO4/(Ks_SO4+SO4)) # R6 reacorgC36 <- a_SO4*reacOC3i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (Ks_SO4/(Ks_SO4+SO4)) # R6 reacorgC46 <- a_SO4*reacOC4i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (Ks_SO4/(Ks_SO4+SO4)) # R6 reacorgC56 <- a_SO4*reacOC5i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (Ks_SO4/(Ks_SO4+SO4)) # R6 reacorgC66 <- a_SO4*reacOC6i*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(Ks_FeOH3/(Ks_FeOH3+FeOH3a))* (Ks_SO4/(Ks_SO4+SO4)) # R6 reacOC1 <- reacorgC11+reacorgC12+reacorgC13+reacorgC14+reacorgC15+reacorgC16 reacOC2 <- reacorgC21+reacorgC22+reacorgC23+reacorgC24+reacorgC25+reacorgC26 reacOC3 <- reacorgC31+reacorgC32+reacorgC33+reacorgC34+reacorgC35+reacorgC36 reacOC4 <- reacorgC41+reacorgC42+reacorgC43+reacorgC44+reacorgC45+reacorgC46 reacOC5 <- reacorgC51+reacorgC52+reacorgC53+reacorgC54+reacorgC55+reacorgC56 reacOC6 <- reacorgC61+reacorgC62+reacorgC63+reacorgC64+reacorgC65+reacorgC66 reacorgCa1 <- reacorgC11+reacorgC21 reacorgCa2 <- reacorgC12+reacorgC22 reacorgCa3 <- reacorgC13+reacorgC23 reacorgCa4 <- reacorgC14+reacorgC24 reacorgCa5 <- reacorgC15+reacorgC25 reacorgCa6 <- reacorgC16+reacorgC26 reacorgCb1 <- reacorgC31+reacorgC41+reacorgC51+reacorgC61 reacorgCb2 <- reacorgC32+reacorgC42+reacorgC52+reacorgC62 reacorgCb3 <- reacorgC33+reacorgC43+reacorgC53+reacorgC63 reacorgCb4 <- reacorgC34+reacorgC44+reacorgC54+reacorgC64 reacorgCb5 <- reacorgC35+reacorgC45+reacorgC55+reacorgC65 reacorgCb6 <- reacorgC36+reacorgC46+reacorgC56+reacorgC66 reacorgCa <- reacOC1+reacOC2 reacorgCb <- reacOC3+reacOC4+reacOC5+reacOC6 reacorgC4P <- (reacOC1i+reacOC2i+reacOC3i+reacOC4i+reacOC5i+reacOC6i)*(Ks_O2/(Ks_O2+O2))*(Ks_NO3/(Ks_NO3+NO3))*(Ks_MnO2/(Ks_MnO2+MnO2a))*(irPa/(Ks_FeOH3+FeOH3a)) # R4 reac7 <- k7*O2*SumNH4 # R7 reac8 <- k8*O2*MnII # R8 reac9 <- k9*O2*FeII # R9 reac10 <- k10*O2*FeS # R10 reac11 <- k11*O2*FeS2 # R11 reac12 <- k12*O2*SumH2S # R12 reac13 <- k13*O2*CH4 # R13 reac14a <- k14*MnO2a*FeII # R14a reac14b <- k14*MnO2b*FeII # R14b reac14 <- reac14a+reac14b # R14 reac15a <- k15*MnO2a*SumH2S # R15a reac15b <- k15*MnO2b*SumH2S # R15a reac15 <- reac15a+reac15b # R15 reac16a <- k16*FeOH3a*SumH2S # R16a reac16b <- k16*FeOH3b*SumH2S # R16b reac16aP <- k16*irPa*SumH2S # R16a reac16bP <- k16*irPb*SumH2S # R16b reac16 <- reac16a+reac16b # R16 reac17 <- k17*FeII*SumH2S # R17 reac18 <- k18*SO4*CH4 # R18 reac19 <- k19*S0 # R19 reac20 <- k20*FeS*S0 # R20 reac21 <- k21*FeOH3a # R21 reac21P <- k21*irPa # R21 reac22 <- k22*MnO2a # R22 reac28 <- k28*FeS*SumH2S # R28 reac29 <- k.biotite*J.biot/(1-por.0)/v.0 # R29 reac30 <- magn*k.magn*((SumH2S*1000)^0.5) # R30 reac32 <- k32*SurfFe*O2 # R32 #------------------------ authigenic carbonates ------------------------------------------------------------------- omegaa <- (CO3 * Ca)/Kspa reac23 <- ifelse (omegaa>1,Kpmaxa* (omegaa-1),-Kdissa*(1-omegaa)*arag ) # R23 omegac <- (CO3 * Ca)/Kspc reac24 <- ifelse (omegac>1,Kpmaxc* (omegac-1),-Kdissc*(1-omegac)*calc) # R24 omegar <- (CO3 * MnII)/Kspr reac26 <- ifelse (omegar>1,Kpmaxr* (omegar-1),-Kdissr*(1-omegar)*rhod) # R26 omegaCFA <- (Ca*rCa)^9.54*(BW_Na*rNa)^0.33*(Mg*rMg)^0.13*(PO4*rPO4)^4.8*(CO3*rCO3)^(1.2-2.3307)*(F*rF)^2.48/KspCFA reac27 <- ifelse(omegaCFA>1,KpmaxCFA* (omegaCFA-1),0) # R27 #------------------------ differential equations ------------------------------------------------------------------- dorgC1 <- tranorgC1-reacOC1 dorgC2 <- tranorgC2-reacOC2 dorgC3 <- tranorgC3-reacOC3 dorgC4 <- tranorgC4-reacOC4 dorgC5 <- tranorgC5-reacOC5 dorgC6 <- tranorgC6-reacOC6 dCFA <- tranCFA + reac27 dO2 <- tranO2 - (1+2*x)*f.grid*reacorgCa1 - (1+2*xx)*f.grid*reacorgCb1 -2*reac7 - 0.5*reac8 - 0.25*reac9 - 2*f.grid*reac10 - 3.5*f.grid*reac11 - 2*reac12 - 2*reac13 - 0.25*f.grid*reac32 dNO3 <- tranNO3 + x*f.grid*reacorgCa1 + xx*f.grid*reacorgCb1 - (0.8+0.6*x)*f.grid*reacorgCa2 - (0.8+0.6*xx)*f.grid*reacorgCb2 + reac7 dMnO2a <- tranMnO2a - 2*reacorgCa3 - 2*reacorgCb3 + 1/f.grid*reac8 - reac14a - reac15a - reac22 dFeOH3a <- tranFeOH3a - 4*reacorgCa4 - 4*reacorgCb4 + 1/f.grid*reac9 + 2*reac14 - 2*reac16a - reac21 + reac32 dSO4 <- tranSO4 - 0.5*f.grid*(reacorgCa5+reacorgCb5) + f.grid*reac10 + 2*f.grid*reac11 + reac12 - reac18 + f.grid*reac19 dCH4 <- tranCH4 + 0.5*f.grid*(reacorgCa6+reacorgCb6)- reac13 - reac18 dSumNH4 <- tranSumNH4 + f.grid*(reacorgCa3*x+reacorgCb3*xx+reacorgCa4*x+reacorgCb4*xx+reacorgCa5*x+reacorgCb5*xx+reacorgCa6*x+reacorgCb6*xx) - reac7 dMnII <- tranMnII + 2*f.grid*(reacorgCa3+reacorgCb3) - reac8 + f.grid*reac14 + f.grid*reac15 - f.grid*reac26 RTFeII <- tranFeII + 4*f.grid*(reacorgCa4+reacorgCb4) - reac9 + f.grid*reac10 + f.grid*reac11 - 2*f.grid*reac14 + 2*f.grid*reac16 - reac17 + 2*f.grid*reac29 + 3*f.grid*reac30 RTSurfFe <- tranSurfFe - reac32 dFeII <- 1/(1+kFe)*RTFeII + 1/(1+kFe)*f.grid*RTSurfFe dFeS <- tranFeS - reac10 + 1/f.grid*reac17 - reac20 - reac28 dFeS2 <- tranFeS2 - reac11 + reac20 + reac28 dSumH2S <- tranSumH2S + 0.5*f.grid*(reacorgCa5+reacorgCb5) - reac12 - f.grid*reac15 - f.grid*reac16 - reac17 + reac18 + 3*f.grid*reac19 - f.grid*reac28 - f.grid*reac30 dMnO2b <- tranMnO2b - reac14b - reac15b + reac22 dSumH2PO4 <- tranSumH2PO4 + f.grid*(reacorgCa*y+reacorgCb*yy) - 4.8*f.grid*reac27 + 4*f.grid*reacorgC4P - r0*reac9*SumH2PO4/(SumH2PO4+KPO4) - 2*r0*f.grid*reac14*SumH2PO4/(SumH2PO4+KPO4) + 2*f.grid*(reac16aP+reac16bP) - r0*f.grid*reac32*SumH2PO4/(SumH2PO4+KPO4) dS0 <- tranS0 + reac15 + reac16 - 4*reac19 - reac20 + reac30 dFeOH3b <- tranFeOH3b - 2*reac16b + reac21 dF <- tranF - 2.48*f.grid*reac27 dSumCO2 <- tranSumCO2 + f.grid*(reacorgCa + reacorgCb - 0.5*reacorgCa6 - 0.5*reacorgCb6) + reac13 + reac18 - f.grid*reac23 - f.grid*reac24 - f.grid*reac26 - 1.2*f.grid*reac27 dTP <- tranTP + (2+x)*f.grid*reacorgCa1 +(2+xx)*f.grid*reacorgCb1 + 1.2*f.grid*(reacorgCa2 + reacorgCb2) - 0.6*f.grid*(reacorgCa2*x + reacorgCb2*xx) - 2*f.grid*(reacorgCa3 + reacorgCb3) - 6*f.grid*(reacorgCa4 + reacorgCb4) + 1.5*f.grid*(reacorgCa5 + reacorgCb5) + f.grid*(reacorgCa6 + reacorgCb6) + 3*f.grid*(reacorgCa*y+reacorgCb*yy) + 12*f.grid*reacorgC4P + reac7 + 2*reac8 + (2-3*r0*SumH2PO4/(SumH2PO4+KPO4))*reac9 + 2*f.grid*reac11 + reac12 + 2*reac13 + (2-6*r0*SumH2PO4/(SumH2PO4+KPO4))*f.grid*reac14 - 3*f.grid*reac15 - 5*f.grid*(reac16a+reac16b) + 6*f.grid*(reac16aP+reac16bP) + reac17 + reac18 + 5*f.grid*reac19 - f.grid*reac28 - 7*f.grid*reac29 - 7*f.grid*reac30 + (1-3*r0*SumH2PO4/(SumH2PO4+KPO4))*f.grid*reac32 + kFe/(1+kFe)*RTFeII - 1/(1+kFe)*f.grid*RTSurfFe dH <- dTP/deltaTP.deltaH - (dSumCO2*deltaTP.deltaSumCO2 + dSumH2S*deltaTP.deltaSumH2S + dSumH2PO4*deltaTP.deltaSumH2PO4 + dSumNH4*deltaTP.deltaSumNH4)/deltaTP.deltaH dCa <- tranCa - f.grid*reac23 - f.grid*reac24 - 9.54*f.grid*reac27 darga <- tranarag + reac23 dcalc <- trancalc + reac24 drhod <- tranrhod + reac26 dMg <- tranMg - 0.13*f.grid*reac27 dmagn <- tranmagn -reac30 dirPa <- tranirPa - 4*reacorgC4P + r0*1/f.grid*reac9*SumH2PO4/(SumH2PO4+KPO4) + 2*r0*reac14*SumH2PO4/(SumH2PO4+KPO4) - 2*reac16aP - reac21P + r0*reac32*SumH2PO4/(SumH2PO4+KPO4) dirPb <- tranirPb - 2*reac16bP + reac21P #------------------------ outputs ------------------------------------------------------------------------------------------------------------------------------ return(list(c(dorgC1=dorgC1,dorgC2=dorgC2,dorgC3=dorgC3,dorgC4=dorgC4,dorgC5=dorgC5,dorgC6=dorgC6, dCFA=dCFA,dO2=dO2, dNO3=dNO3, dMnO2a=dMnO2a,dFeOH3a=dFeOH3a, dSO4=dSO4,dCH4=dCH4,dSumNH4=dSumNH4,dMnII=dMnII,dFeII=dFeII,dFeS=dFeS,dFeS2=dFeS2,dSumH2S=dSumH2S,dMnO2b=dMnO2b,dSumH2PO4=dSumH2PO4,dS0=dS0,dFeOH3b=dFeOH3b,dF=dF,dSumCO2=dSumCO2,dH=dH,dCa=dCa,darga=darga, dcalc=dcalc,drhod=drhod,dMg=dMg,dmagn=dmagn,dirPa=dirPa,dirPb=dirPb),omegaCFA=omegaCFA,PO4=PO4,HPO4=HPO4,H2PO4=H2PO4,H3PO4=H3PO4)) } # ---------------------------------------------------------- # # --------------------- model solution --------------------- # # ---------------------------------------------------------- # # ---------------------initial guess of state variables--------------------- orgC1 <- rep(0,N) orgC2 <- rep(0,N) orgC3 <- rep(0,N) orgC4 <- rep(0,N) orgC5 <- rep(0,N) orgC6 <- rep(0,N) CFA <- rep(0,N) O2 <- rep(0,N) NO3 <- rep(0,N) MnO2a <- rep(0,N) FeOH3a <- rep(0,N) SO4 <- rep(0,N) CH4 <- rep(0,N) SumNH4 <- rep(0,N) MnII <- rep(0,N) FeII <- rep(0,N) FeS <- rep(0,N) FeS2 <- rep(0,N) SumH2S <- rep(0,N) MnO2b <- rep(0,N) SumH2PO4 <- rep(0,N) S0 <- rep(0,N) FeOH3b <- rep(0,N) F <- rep(0,N) SumCO2 <- rep(BW_SumCO2,N) H <- rep(BW_H,N) Ca <- rep(0,N) arga <- rep(0,N) calc <- rep(0,N) rhod <- rep(0,N) Mg <- rep(0,N) magn <- rep(0,N) irPa <- rep(0,N) irPb <- rep(0,N) Conc <- c(orgC1,orgC2,orgC3,orgC4,orgC5,orgC6,CFA,O2,NO3,MnO2a,FeOH3a,SO4,CH4,SumNH4,MnII,FeII,FeS,FeS2,SumH2S,MnO2b,SumH2PO4,S0,FeOH3b,F,SumCO2,H,Ca,arga,calc,rhod,Mg,magn,irPa,irPb) #----------------------- run model --------------------------------------------------------------------------------------------------------------- options(max.print=20000) times <- seq(0, 20000, by = 1000) system.time (sol <- ode.1D (y=Conc, time=times, func=Pmodel, maxstep = 10000, parms = NULL,nspec=34,method="vode",rtol = 1e-14, atol = 1e-14,verbose =TRUE,mf=-22)) # show run time, sol contains the solution #---------------------- collecting results ------------------------------------------------------------------------------------------- orgC1 <- sol[21, (2):(N+1)] orgC2 <- sol[21, (N + 2):(2 * N + 1)] orgC3 <- sol[21, (2*N + 2):(3 * N + 1)] orgC4 <- sol[21, (3*N + 2):(4 * N + 1)] orgC5 <- sol[21, (4*N + 2):(5 * N + 1)] orgC6 <- sol[21, (5*N + 2):(6 * N + 1)] CFA <- sol[21, (6*N + 2):(7 * N + 1)] O2 <- sol[21, (7*N + 2):(8 * N + 1)] NO3 <- sol[21, (8*N + 2):(9 * N + 1)] MnO2a <- sol[21, (9*N + 2):(10 * N + 1)] FeOH3a <- sol[21, (10*N + 2):(11 * N + 1)] SO4 <- sol[21, (11*N + 2):(12 * N + 1)] CH4 <- sol[21, (12*N + 2):(13 * N + 1)] SumNH4 <- sol[21, (13*N + 2):(14 * N + 1)] MnII <- sol[21, (14*N + 2):(15 * N + 1)] FeII <- sol[21, (15*N + 2):(16 * N + 1)] FeS <- sol[21, (16*N + 2):(17 * N + 1)] FeS2 <- sol[21, (17*N + 2):(18 * N + 1)] SumH2S <- sol[21, (18*N + 2):(19 * N + 1)] MnO2b <- sol[21, (19*N + 2):(20 * N + 1)] SumH2PO4 <- sol[21, (20*N + 2):(21 * N + 1)] S0 <- sol[21, (21*N + 2):(22 * N + 1)] FeOH3b <- sol[21, (22*N + 2):(23 * N + 1)] F <- sol[21, (23*N + 2):(24 * N + 1)] SumCO2 <- sol[21, (24*N + 2):(25 * N + 1)] H <- sol[21, (25*N + 2):(26 * N + 1)] Ca <- sol[21, (26*N + 2):(27 * N + 1)] arga <- sol[21, (27*N + 2):(28 * N + 1)] calc <- sol[21, (28*N + 2):(29 * N + 1)] rhod <- sol[21, (29*N + 2):(30 * N + 1)] Mg <- sol[21, (30*N + 2):(31 * N + 1)] magn <- sol[21, (31*N + 2):(32 * N + 1)] irPa <- sol[21, (32*N + 2):(33 * N + 1)] irPb <- sol[21, (33*N + 2):(34 * N + 1)] omegaCFA <- sol[21, (34*N + 2):(35 * N + 1)] PO4 <- sol[21, (35*N + 2):(36 * N + 1)] HPO4 <- sol[21, (36*N + 2):(37 * N + 1)] H2PO4 <- sol[21, (37*N + 2):(38 * N + 1)] H3PO4 <- sol[21, (38*N + 2):(39 * N + 1)] # --------------------- convert results ------------------------------ orgCat <- (orgC1+orgC2)*12/1000/rho_sed*100 orgCbt <- (orgC3+orgC4+orgC5+orgC6)*12/1000/rho_sed*100 orgCct <- (J.orgCc/0.2/0.36)*12/1000/rho_sed*100 orgCt <- orgCat + orgCbt + orgCct orgNt <- orgCat*x + orgCbt*xx + orgCct*xx CFA_Pt <- CFA*4.8*1000/rho_sed orgPt <- ((orgC1+orgC2)*y + (orgC3+orgC4+orgC5+orgC6+J.orgCc/0.2/0.36)*yy)*1000/rho_sed O2t <- O2*1e6 NO3t <- NO3*1e6 MnO2t <- (MnO2a + MnO2b)*1000/rho_sed + 88/55 FeOH3t <- (FeOH3a + FeOH3b)*1000/rho_sed SO4t <- SO4*1e3 CH4t <- CH4*1e6 SumNH4t <- SumNH4*1e3 MnIIt <- MnII*1e6 FeIIt <- FeII*1e6 SumH2St <- SumH2S*1e6 SumH2PO4t <- SumH2PO4*1e6 S0t <- S0*1000/rho_sed SumCO2t <- SumCO2*1e3 pH <- -log10(H) Cat <- Ca*1e3 magnt <- magn*56*100/1000/2.5 FeS.Fe <- FeS*100/2.5*56/1000 FeS2.Fe <- FeS2*100/2.5*56/1000 Ft <- F*1e6 irPat <- irPa*1000/rho_sed irPbt <- irPb*1000/rho_sed irPt <- irPat + irPbt calct <- calc*100/2.5*12/1000 C_Pratio <- (orgC1+orgC2+orgC3+orgC4+orgC5+orgC6+J.orgCc/0.2/0.36)/((orgC1+orgC2)*y + (orgC3+orgC4+orgC5+orgC6+J.orgCc/0.2/0.36)*yy) # --- output of data --- # mydf<-as.data.frame(cbind(pH,orgCt,FeOH3t,SO4t,CH4t,SumNH4t,FeIIt,SumH2St,orgPt,SumH2PO4t,SumCO2t,Cat,MnO2t,magnt,Ft,CFA_Pt,MnIIt,FeS.Fe,FeS2.Fe,calct,irPt,irPat,irPbt,omegaCFA)) mydf