### 0
#####
A <- aquaenv(35, 15, SumCO2=0.003, TA=seq(0.001,0.004, 0.0001))
plot(A, xval=A$TA, xlab="[TA]/(mol/kg-soln)")
plot(A, what=c("CO2", "HCO3", "CO3"), bjerrum=TRUE, log=TRUE)
plot(A, xval=A$TA, xlab="[TA]/(mol/kg-soln)", what=c("CO2", "HCO3", "CO3"),
cumulative=TRUE, ylab="mol/kg-soln", ylim=c(0,0.0031))
### 1
#####
SumCO2 <- 0.0020
pH <- 8
S <- 30
t <- 1:15
p <- 10
ae <- aquaenv(S, t, p, SumCO2=SumCO2, pH=pH, revelle=TRUE, dsa=TRUE)
plot(ae, xval=t, xlab="T/(deg C)", newdevice=FALSE)
### 2
#####
S <- 35
t <- 15
SumCO2 <- 0.003500
SumNH4 <- 0.000020
mass_sample <- 0.01 # the mass of the sample solution in kg
mass_titrant <- 0.02 # the total mass of the added titrant solution in
# kg
conc_titrant <- 0.01 # the concentration of the titrant solution in
# mol/kg-soln
S_titrant <- 0.5 # the salinity of the titrant solution (the
# salinity of a solution with a ionic strength of
# 0.01 according to: I = (19.924 S) / (1000 - 1.005S)
steps <- 50 # the amount of steps the mass of titrant is added
# in
type <- "HCl"
pHstart <- 11.3
ae <- titration(aquaenv(S=S, t=t, SumCO2=SumCO2, SumNH4=SumNH4,
pH=pHstart), mass_sample, mass_titrant, conc_titrant,
S_titrant, steps, type)
# plotting everything
plot(ae, xval=ae$delta_mass_titrant, xlab="HCl solution added [kg]",
mfrow=c(10,10))
# plotting selectively
size <- c(12,8) #inches
mfrow <- c(4,4)
what <- c("TA", "pH", "CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH",
"NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F", "pCO2")
plot(ae, xval=ae$delta_mass_titrant, xlab="HCl solution added [kg]",
what=what, size=size, mfrow=mfrow)
plot(ae, xval=ae$pH, xlab="free scale pH", what=what, size=size,
mfrow=mfrow)
# different x values
plot(ae, xval=ae$delta_conc_titrant, xlab="[HCl] offset added
[mol/kg-soln]", what=what, size=size, mfrow=mfrow)
plot(ae, xval=ae$delta_moles_titrant, xlab="HCl added [mol]", what=what,
size=size, mfrow=mfrow, newdevice=FALSE)
# bjerrum plots
plot(ae, bjerrum=TRUE)
what <- c("CO2", "HCO3", "CO3")
plot(ae, what=what, bjerrum=TRUE)
plot(ae, what=what, bjerrum=TRUE, lwd=4, palette=c("cyan", "magenta",
"yellow"), bg="gray", legendinset=0.1, legendposition="topleft")
what <- c("CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH", "NH4", "NH3",
"H2SO4", "HSO4", "SO4", "HF", "F")
plot(ae, what=what, bjerrum=TRUE, log=TRUE, newdevice=FALSE)
plot(ae, what=what, bjerrum=TRUE, log=TRUE, ylim=c(-6,-1),
legendinset=0, lwd=3, palette=c(1,3,4,5,6,colors()[seq(100,250,6)]))
### 3
#####
parameters <- list(
t = 15 , # degrees C
S = 35 , # psu
SumCO2_t0 = 0.002 , # mol/kg-soln (comparable to Wang2005)
TA_t0 = 0.0022 , # mol/kg-soln (comparable to Millero1998)
kc = 0.5 , # 1/d proportionality factor
# for air-water exchange
kp = 0.000001 , # mol/(kg-soln*d) max rate of calcium
# carbonate precipitation
n = 2.0 , # - exponent for kinetic
# rate law of precipitation
modeltime = 20 , # d duration of the model
outputsteps = 100 # number of outputsteps
)
boxmodel <- function(timestep, currentstate, parameters)
{
with (
as.list(c(currentstate,parameters)),
{
ae <- aquaenv(S=S, t=t, SumCO2=SumCO2, pH=-log10(H), SumSiOH4=0,
SumBOH3=0, SumH2SO4=0, SumHF=0, dsa=TRUE)
Rc <- kc * ((ae$CO2_sat) - (ae$CO2))
Rp <- kp * (1-ae$omega_calcite)^n
dSumCO2 <- Rc - Rp
dHRc <- ( -(ae$dTAdSumCO2*Rc ))/ae$dTAdH
dHRp <- (-2*Rp -(ae$dTAdSumCO2*(-Rp)))/ae$dTAdH
dH <- dHRc + dHRp
ratesofchanges <- c(dSumCO2, dH)
processrates <- c(Rc=Rc, Rp=Rp)
outputvars <- c(dHRc=dHRc, dHRp=dHRp)
return(list(ratesofchanges, list(processrates, outputvars, ae)))
}
)
}
with (as.list(parameters),
{
aetmp <- aquaenv(S=S, t=t, SumCO2=SumCO2_t0,
TA=TA_t0, SumSiOH4=0, SumBOH3=0,
SumH2SO4=0, SumHF=0)
H_t0 <- 10^(-aetmp$pH)
initialstate <<- c(SumCO2=SumCO2_t0, H=H_t0)
times <<- seq(0,modeltime,(modeltime/outputsteps))
output <<- as.data.frame(vode(initialstate,times,
boxmodel,parameters, hmax=1))
})
what <- c("SumCO2", "TA", "Rc", "Rp",
"omega_calcite", "pH", "dHRc", "dHRp")
plot(aquaenv(ae=output, from.data.frame=TRUE), xval=output$time,
xlab="time/d", mfrow=c(3,3), size=c(15,10), what=what)
what <- c("dHRc", "dHRp")
plot(aquaenv(ae=output, from.data.frame=TRUE), xval=output$time,
xlab="time/d", what=what, ylab="mol-H/(kg-soln*d)",
legendposition="topright", cumulative=TRUE)Run the code above in your browser using DataLab