\documentclass[article,nojss,shortnames]{jss}
\DeclareGraphicsExtensions{.pdf,.eps}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{longtable}
\usepackage{array}
\setlength{\extrarowheight}{0.1cm}
\newcommand{\noun}[1]{\textsc{#1}}
\newcommand{\aq}{\textbf{\textsf{AquaEnv}}}
\newcommand{\lm}{\textbf{\textsf{minpack.lm}}}
\newcommand{\R}{\proglang{R }}
\newcommand{\ds}{\textbf{\textsf{deSolve}}}
\setkeys{Gin}{width=\textwidth} %width of graphics produced by Sweave
% enumerate and add to table down to subparagraph
\addtocounter{tocdepth}{+3}
\addtocounter{secnumdepth}{+3}
\title{Package \aq: an \noun{Aqua}tic modelling \noun{Env}ironment in \proglang{R}}
\Plaintitle{Package AquaEnv an Aquatic modelling Environment in R}
\Keywords{aquatic modelling, pH, pH scales, dissolved inorganic carbon, total alkalinity, total alkalinity curve fitting, theoretical titration, revelle factor, omega, solubility products, $\rm CO_2$, ocean acidification, estuaries, carbonate system, seawater, buffer factors, \proglang{R}}
\Plainkeywords{aquatic modelling, pH, pH scales, dissolved inorganic carbon, total alkalinity, total alkalinity curve fitting, theoretical titration, revelle factor, omega, solubility products, CO2, ocean acidification, estuaries, carbonate system, seawater, buffer factors, R}
\author{Andreas F. Hofmann, Karline Soetaert and Mathilde Hagens}
\Plainauthor{Andreas F. Hofmann, Karline Soetaert and Mathilde Hagens}
\Abstract{
\noindent
\code{AquaEnv} \citep{AquaEnv} is an integrated development toolbox for aquatic chemical model generation focused on (ocean) acidification and $\rm CO_2$ air-water exchange. It contains
\begin{itemize}
\item All elements necessary to model the pH, the related $\rm CO_2$ air-water exchange, as well as aquatic acid-base chemistry in general for an arbitrary marine, estuarine or freshwater system.
\item A suite of tools to visualize this information.
\item It can be used to build dynamic models of aquatic systems that include acid-base chemistry.
\item The sensitivity of the system to variations in the input variables can be visualized.
\item A number of example ``applications'' that make use of \code{AquaEnv} are:
\begin{itemize}
\item a theoretical titration simulator
\item a routine to determine total alkalinity ([TA]), the total dissolved inorganic carbon concentration ($[\rm \sum CO_2]$), as well as additionally the electrode standard potential ($\rm E_0$) and the first dissociation constant of the carbonate system ($\rm K^*_{CO_2}$) from titration data.
\end{itemize}
\end{itemize}
}
\Address{
Karline Soetaert\\
Royal Netherlands Institute\\
of Sea Research (NIOZ)\\
Yerseke, The Netherlands\\
and\\
Mathilde Hagens\\
Utrecht University
}
%% need no \usepackage{Sweave.sty}
%\VignetteIndexEntry{AquaEnv}
\begin{document}
%\SweaveOpts{concordance=TRUE}
\SweaveOpts{engine=R,eps=FALSE}
%\SweaveOpts{keep.source=TRUE}
<<Preliminaries,echo=FALSE,results=hide>>=
library("AquaEnv")
library("deSolve")
options(width=80)
Plotit <- AquaEnv:::plot.aquaenv
plot.aquaenv <- function(...) {
if ("newdevice" %in% list(...))
Plotit(...)
else
Plotit(newdevice = FALSE, ...) # else Sweave does not work...
}
@
\maketitle
\newpage
\tableofcontents
\newpage
\section{Introduction}
\code{AquaEnv} can be used in three ways
\begin{itemize}
\item It provides separate functions to calculate the stoichiometric equilibrium constants ($\rm K^*$) for key acid base systems in natural seawater, the Henry's constants ($\rm K_0$), as well as the solubility products ($\rm K_{sp}$) for calcite and aragonite.
This functionality is provided via the functions \code{K\_CO2, K\_HCO3, K\_BOH3, K\_W},\\
\code{K\_HSO4, K\_HF, K\_NH4, K\_H2S, K\_H3PO4, K\_H2PO4, K\_HPO4, K\_SiOH4} \\,
\code{K\_SiOOH3, K0\_CO2, K0\_O2, Ksp\_aragonite} and \code{Ksp\_calcite}.
\item It is also possible to estimate all the acid-base chemistry using one single function: \code{aquaenv}. This function returs a list of class \code{aquaenv} that contains, in addition to the input parameters:
\begin{itemize}
\item the clorinity, the ionic strength, $[\rm \sum B(OH)_3]$, $[\rm \sum H_2SO_4]$, $[\rm \sum HF]$, $[\rm Cl^-]$, $[\rm Cl^-]$, \\
$[\rm \sum Br]$, $[\rm Na^+]$, $[\rm Mg^{2+}]$, $[\rm Ca^{2+}]$, $[\rm K^{+}]$, $[\rm Sr^{2+}]$ calculated from salinity as given in \cite{DOE1994}
\footnote{Please note that if values for $[\rm \sum B(OH)_3]$, $[\rm \sum H_2SO_4]$, and/or $[\rm \sum HF]$ are given as input parameters, these parameters are used and not the ones calculated from salinity.}
\item the gauge pressure p \citep[total pressure minus atmospheric pressure,][]{Feistel2008} either given as input variable, or calculated from depth \citep[according to][]{Fofonoff1983}, or calculated from the total pressure P and the atmospheric pressure Pa,
both of which can be given as input variables and are also stored in an object of class \code{aquaenv}
\item the seawater density calculated from temperature and salinity as given by \cite{Millero1981}
\item a set of conversion factors to convert between different pH scales \citep{Dickson1984, Zeebe2001} and between mol/kg-$\rm H_2O$ and mol/kg-solution (inferred from \citet{Roy1993b} and \citet{DOE1994})
\item the Henry's constants for $\rm CO_2$ \citep{Weiss1974} and for $\rm O_2$ \citep[inferred from][]{Weiss1970} calculated from temperature and salinity as well as the associated saturation concentrations of $\rm CO_2$ and $\rm O_2$.
\item the ion product of water \citep{Millero1995}, the stoichiometric equilibrium constants of $\rm HSO_4^-$ \citep{Dickson1990}, $\rm HF$ \citep{Dickson1979a} ,
$\rm CO_2$ \citep{Roy1993b}, $\rm HCO_3^-$ \citep{Roy1993b}, $\rm B(OH)_3$ \citep{Dickson1990}, $\rm NH_4^+$\citep{Millero1995a}, $\rm H_2S$ \citep{Millero1995},
$\rm H_3PO_4$\citep{Millero1995}, $\rm H_2PO_4^-$ \citep{Millero1995}, $\rm HPO_4^{2-}$ \citep{Millero1995}, $\rm Si(OH)_4$ \citep{Millero1988},
$\rm SiO(OH)_3^-$ \citep{Wischmeyer2003}, $\rm HNO_2$ \citep{Riordan2005}, $\rm HNO_3$, $\rm H_2SO_4$ \citep{Atkins1996}, $\rm HS$ \citep{Atkins1996} mostly calculated
as functions of temperature and salinity and pressure corrected according to \cite{Millero1995}.
\item the solubility products of calcite and aragonite \citep{Mucci1983} as well as the associated $\Omega$'s if a full speciation is calculated (see below)
\item the fugacity of $\rm CO_2$ - if a full speciation is calculated (see below)
\item if $[\rm \sum CO_2]$ and pH are given $[\rm TA]$ is calculated, if $[\rm \sum CO_2]$ and $[\rm TA]$ are given pH is calculated, if $[\rm \sum CO_2]$ and $\rm [CO_2]$ or f$\rm CO_2$
are given, pH and $[\rm TA]$ are calculated.
\item if either one of the pairs pH and $\rm [CO_2]$ or f$\rm CO_2$, pH and $[\rm TA]$, or $[\rm TA]$ and $\rm [CO_2]$ or f$\rm CO_2$ is given, $[\rm \sum CO_2]$ is calculated
\item if sufficient information is given and the flag \code{speciation=TRUE} is set, a full speciation of $[\rm \sum CO2]$, $[\rm \sum NH4]$, $[\rm \sum H_2S]$, $[\rm \sum HNO3]$, $[\rm \sum HNO2]$,
$[\rm \sum H_3PO4]$, \\
$[\rm \sum Si(OH)_4]$, $[\rm \sum B(OH)_3]$, $[\rm \sum H_2SO_4]$, $[\rm \sum HF]$, as well as water itself is calculated
\item if the flag \code{dsa = TRUE} is set, all necessary quantities for the explicit ``direct substitution approach'' (DSA) to pH modelling as given in \cite{Hofmann2008} are
calculated. These are the buffer factor (the partial derivative of $\rm [TA]$ with respect to $\rm [H^+]$) and the partial derivatives of $\rm [TA]$ with respect to the other
total quantities. Furthermore, the partial derivatives of $\rm [TA]$ with respect to changes in the equilibrium constants ($K^*$), multiplied with the partial derivatives
of the equilibrium constants with respect to their variables needed for the DSA with time variable equilibrium constants as described in \citet{Hofmann2009} are calculated.
Finally, the ionization fractions as defined by \cite{Stumm1996} and used in \cite{Hofmann2010a} are calculated for the full speciation.
\item Thirdly, a generic function \code{BufferFactors} is provided, based on \cite{Hagens2016}. This function internally calls the function \code{aquaenv} and uses its output to analytically calculate the sensitivity of pH and concentrations of $\rm CO_2$ and other acid-base species to a change in ocean chemistry, as well as the Revelle factor.
\end{itemize}
\item Input for \code{aquaenv} and \code{BufferFactors} has to be supplied in standard SI units, the free proton pH scale and in molinity (mol/kg-solution)\footnote{Note that it is not sufficient to give a gravimetric concentration in
mol/kg since there is a substantial difference between mol/kg-$\rm H_2O$ (molality) and mol/kg-solution (molinity).}.
Conversion of input parameters to these necessary units and pH scale can be done with the generic function \code{convert}.
\item The information created with \code{aquaenv} is also supplied in standard SI units and in molinity. All elements of an object of class \code{aquaenv} of a certain unit or pH scale
can be converted into other units or pH scales with the function \code{convert} as well.
\item One can use input vectors of salinity S, temperature t, or gauge pressure p (as well as total pressure P and depth d) for \code{aquaenv} to obtain vectors of all calculated information as function of
the input vector. This can be visualized in a large variety of ways using the \code{plot} function specially defined for objects of class \code{aquaenv}.
\item Objects of class \code{aquaenv} can be used in dynamic models to define the state of the system in each timestep of the numerical integration (done with e.g. \ds).
With the function \code{aquaenv} and the flag \code{from.data.frame=TRUE} it is possible to convert output of those dynamic models into objects of class
\code{aquaenv} which allows the user to use the whole suite of visualisation tools that is provided by the function \code{plot} in \aq.
\item \cite{Hofmann2008}, \cite{Hofmann2009}, and \cite{Hofmann2010a} describe methods for an ``explicit'' pH modelling which allows
for the quantification of the influences of kinetically modelled processes on pH. Objects of class \code{aquaenv} provide all needed quantities
(partial derivatives of $\rm[TA]$, ionization fractions, etc.) to employ both of those methods in dynamic models.
\aq$\,$ also provides the functionality to cumulatively plot the obtained influences on pH.
\item As an example of how to use the aquatic chemical toolbox that is provided by \aq, two applications are provided:
\begin{itemize}
\item The function \code{titration}: creates theoretical titrations which can be used to e.g. create Bjerrum plots with the function \code{plot.aquaenv} in \aq.
\item The function \code{TAfit}: a routine based on a method in \cite{DOE1994} that makes use of that theoretical titration function and allows for
determining total alkalinity ($\rm [TA]$), the total dissolved inorganic carbon concentration ([$\sum$CO2]), as well as additionally the electrode standard potential ($\rm E_0$) and the first dissociation
constant of the carbonate system ($\rm K^*_{CO_2}$) using the Levenberg-Marquart algorithm (least squares optimization procedure) as provided in the \R package \lm.
\end{itemize}
\end{itemize}
\section[The elements of an object of class aquaenv]{The elements of an object of class \code{aquaenv}}
The function \code{aquaenv}, the central function of \aq, returns an object of class \code{aquaenv}. This object is a list of different elements which can be accesses with the \$ character or with the [[]] operator
<<AquaenvElements, fig=FALSE, echo=TRUE>>=
test <- aquaenv(S = 35, t = 10)
test$t
@
If enough input data is supplied to define the pH of the system and the flags \code{speciation} and \code{dsa} are \code{TRUE} while the flag
\code{skeleton} is \code{FALSE}, an object of class \code{aquaenv} contains the following elements\footnote{Literature references are given in Appendix \ref{app: references}.
}
\begin{footnotesize}
\begin{longtable}{l|l|p{7cm}}
\textbf{element}& \textbf{unit} & \textbf{explanation} \\ \hline
\code{S} & ``psu'' (no unit) & salinity \\
\code{t} & \textdegree C & temperature \\
\code{p} & bar & gauge pressure \citep[total pressure minus atmospheric pressure,][]{Feistel2008} \\
\code{T} & K & absolute temperature \\
\code{Cl} & \textperthousand & chlorinity \\
\code{I} & mol/kg-$\rm H_2O$ & ionic strength \\
\code{P} & bar & total pressure \\
\code{Pa} & bar & atmospheric pressure \\
\code{d} & m & depth \\
\code{density} & kg/$\rm m^3$ & (seawater) density \\
\code{SumCO2} & mol/kg-soln & $[\rm \sum CO_2]$, total dissolved inorganic carbon concentration \\
\code{SumNH4} & mol/kg-soln & $[\rm \sum NH_4^+]$, total ammonium concentration\\
\code{SumH2S} & mol/kg-soln & $[\rm \sum H_2S]$, total sulfide concentration\\
\code{SumHNO3} & mol/kg-soln & $[\rm \sum HNO_3]$, total nitrate concentration\\
\code{SumHNO2} & mol/kg-soln & $[\rm \sum HNO_2]$, total nitrite concentration\\
\code{SumH3PO4} & mol/kg-soln & $[\rm \sum H_3PO_4]$, total phosphate concentration\\
\code{SumSiOH4} & mol/kg-soln & $[\rm \sum Si(OH)_4]$, total silicate concentration\\
\code{SumBOH3} & mol/kg-soln & $[\rm \sum B(OH)_3]$, total borates concentration\\
\code{SumH2SO4} & mol/kg-soln & $[\rm \sum H_2SO_4]$, total sulfate concentration\\
\code{SumHF} & mol/kg-soln & $[\rm \sum HF]$, total fluoride concentration\\
\code{Br} & mol/kg-soln & $[\rm Br^-]$, bromide concentration\\
\code{ClConc} & mol/kg-soln & $[\rm Cl^-]$, chloride concentration\\
\code{Na} & mol/kg-soln & $[\rm Na^{+}]$, sodium concentration\\
\code{Mg} & mol/kg-soln & $[\rm Mg^{2+}]$, magnesium concentration\\
\code{Ca} & mol/kg-soln & $[\rm Ca^{2+}]$, calcium concentration\\
\code{K} & mol/kg-soln & $[\rm K^+]$, potassium concentration\\
\code{Sr} & mol/kg-soln & $[\rm Sr^{2+}]$, strontium concentration\\
\code{molal2molin} & (mol/kg-soln)/(mol/kg-H2O) & concentration conversion factor: from molality to molinity\\
\code{free2tot} & - & pH conversion factor: free scale to total scale\\
\code{free2sws} & - & pH conversion factor: free scale to sawater scale\\
\code{tot2free} & - & pH conversion factor: total scale to free scale\\
\code{tot2sws} & - & pH conversion factor: total scale to seawater scale\\
\code{sws2free} & - & pH conversion factor: seawater scale to free scale\\
\code{sws2tot} & - & pH conversion factor: seawater scale to total scale\\
\code{K0\_CO2} & mol/(kg-soln*atm) & Henry's constant for $\rm CO_2$\\
\code{K0\_O2} & mol/(kg-soln*atm) & Henry's constant for $\rm O_2$\\
\code{fCO2atm} & atm & atmospheric fugacity of CO$_2$\\
\code{fO2atm} & atm & atmospheric fugacity of O$_2$\\
\code{CO2\_sat} & mol/kg-soln & $\rm CO_2$ saturation concentration at an atmospheric fugacity of \code{fCO2atm}\\
\code{O2\_sat} & mol/kg-soln & $\rm O_2$ saturation concentration at an atmospheric fugacity of \code{fO2atm}\\
\code{K\_W} &(mol/kg-soln)$^2$, free pH scale & stoichiometric equilibrium ion product of $\rm H_2O$:\\
& & $\rm K^*_W = [H^+][OH^-]$\\
\code{K\_HSO4} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant \\
& & $\rm K^*_{HSO_4^-} = [H^+][SO_4^{2-}] / [HSO_4^-]$\\
\code{K\_HF} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{HF} = [H^+][F^{-}] / [HF]$\\
\code{K\_CO2} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{CO_2} = [H^+][HCO_3^{-}] / [CO_2]$\\
\code{K\_HCO3} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant \\
& & $\rm K^*_{HCO_3^{-}} = [H^+][CO_3^{2-}] / [HCO_3^{-}]$\\
\code{K\_BOH3} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{B(OH)_3} = [H^+][B(OH)_4^-] / [B(OH)_3]$\\
\code{K\_NH4} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{NH_4^+} = [H^+][NH_3] / [NH_4^+]$\\
\code{K\_H2S} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant \\
& & $\rm K^*_{H_2S} = [H^+][HS^-] / [H_2S]$\\
\code{K\_H3PO4} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{H_3PO_4} = [H^+][H_2PO_4^-] / [H_3PO_4]$\\
\code{K\_H2PO4} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{H_2PO_4^-} = [H^+][HPO_4^{2-}] / [H_2PO_4^-]$\\
\code{K\_HPO4} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{HPO_4^{2-}} = [H^+][PO_4^{3-}] / [HPO_4^{2-}]$\\
\code{K\_SiOH4} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{Si(OH)_4} = [H^+][SiO(OH)_3^-] / [Si(OH)_4]$\\
\code{K\_SiOOH3} &mol/kg-soln, free pH scale & stoichiometric equilibrium constant\\
& & $\rm K^*_{SiO(OH)_3^-} = [H^+][SiO_2(OH)_2^{2-}] / [SiO(OH)_3^-]$\\
\code{K\_HNO2} &mol/kg-soln; mol/kg-H2O; mol/l & approximate value for equilibrium constant \\
& & $\rm K^*_{HNO_2} = [H^+][NO_2^-] / [HNO_2]$\\
\code{K\_HNO3} &mol/kg-soln; mol/kg-H2O; mol/l & approximate value for equilibrium constant \\
& & $\rm K^*_{HNO_3} = [H^+][NO_3^-] / [HNO_3]$\\
\code{K\_H2SO4} &mol/kg-soln; mol/kg-H2O; mol/l & approximate value for equilibrium constant \\
& & $\rm K^*_{H_2SO_4} = [H^+][HSO_4^-] / [H_2SO_4]$\\
\code{K\_HS} & mol/kg-soln; mol/kg-H2O; mol/l & approximate value for equilibrium constant \\
& & $\rm K^*_{HS^-} = [H^+][S^{2-}] / [HS^-]$\\
\code{Ksp\_calcite}& (mol/kg-soln)$^2$ & stoichiometric equilibrium solubility product of\\
& & calcite: $\rm Ksp^*_{cal} = [Ca^{2+}][CO_3^{2-}]$\\
\code{Ksp\_aragonite} &(mol/kg-soln)$^2$ & stoichiometric equilibrium solubility product of\\
& & aragonite: $\rm Ksp^*_{ara} = [Ca^{2+}][CO_3^{2-}]$\\
\code{TA} & mol/kg-soln & [TA], total alkalinity\\
\code{pH} & -, free scale & pH\\
\code{fCO2} & atm & fugacity of $\rm CO_2$ in the water (i.e. in a small\\
& & volume of air equilibrated with the water)\\
\code{CO2} & mol/kg-soln & $[\rm CO_2]$\\
\code{HCO3} & mol/kg-soln & $[\rm HCO_3^-]$\\
\code{CO3} & mol/kg-soln & $[\rm CO_3^{2-}]$\\
\code{BOH3} & mol/kg-soln & $[\rm B(OH)_3]$\\
\code{BOH4} & mol/kg-soln & $[\rm B(OH)_4^-]$\\
\code{OH} & mol/kg-soln & $[\rm OH^-]$\\
\code{H3PO4} & mol/kg-soln & $[\rm H_3PO_4]$\\
\code{H2PO4} & mol/kg-soln & $[\rm H_2PO_4^-]$\\
\code{HPO4} & mol/kg-soln & $[\rm HPO_4^{2-}]$\\
\code{PO4} & mol/kg-soln & $[\rm PO_4^{3-}]$\\
\code{SiOH4} & mol/kg-soln & $[\rm Si(OH)_4]$\\
\code{SiOOH3} & mol/kg-soln & $[\rm SiO(OH)_3^-]$\\
\code{SiO2OH2} & mol/kg-soln & $[\rm SiO_2(OH)_2^{2-}]$\\
\code{H2S} & mol/kg-soln & $[\rm H_2S]$\\
\code{HS} & mol/kg-soln & $[\rm HS^-]$\\
\code{S2min} & mol/kg-soln & $[\rm S^{2-}]$\\
\code{NH4} & mol/kg-soln & $[\rm NH_4^+]$\\
\code{NH3} & mol/kg-soln & $[\rm NH_3]$\\
\code{H2SO4} & mol/kg-soln & $[\rm H_2SO_4]$\\
\code{HSO4} & mol/kg-soln & $[\rm HSO_4^-]$\\
\code{SO4} & mol/kg-soln & $[\rm SO_4^{2-}]$\\
\code{HF} & mol/kg-soln & $[\rm HF]$\\
\code{F} & mol/kg-soln & $[\rm F^-]$\\
\code{HNO3} & mol/kg-soln & $[\rm HNO_3]$\\
\code{NO3} & mol/kg-soln & $[\rm NO_3^-]$\\
\code{HNO2} & mol/kg-soln & $[\rm HNO_2]$\\
\code{NO2} & mol/kg-soln & $[\rm NO_2^-]$\\
\code{omega\_calcite} & - & saturation state $\Omega$ with respect to calcite\\
\code{omega\_aragonite} & - & saturation state $\Omega$ with respect to aragonite\\
\code{revelle} & - & Revelle factor (redundant in current version of \aq)\\
\code{c1} & - & ionization fraction $c_1 = [\rm CO_2]/[\rm \sum CO_2]$\\
\code{c2} & - & ionization fraction $c_2 = [\rm HCO_3^-]/[\rm \sum CO_2]$\\
\code{c3} & - & ionization fraction $c_3 = [\rm CO_3^{2-}]/[\rm \sum CO_2]$\\
\code{dTAdSumCO2} & - & $\rm \frac{\partial [TA]}{[\partial \sum CO_2]}$\\
& & with $\rm[TA] = f([H^+], [\sum CO_2], ...)$\\
\code{b1} & - & ionization fraction $b_1 = [\rm B(OH)_3]/[\rm \sum B(OH)_3]$\\
\code{b2} & - & ionization fraction $b_2 = [\rm B(OH)_4^-]/[\rm \sum B(OH)_3]$\\
\code{dTAdSumBOH3} & - & $\rm \frac{\partial [TA]}{[\partial \sum B(OH)_3]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{so1} & - & ionization fraction $so_1 = [\rm H_2SO_4]/[\rm \sum H_2SO_4]$\\
\code{so2} & - & ionization fraction $so_2 = [\rm HSO_4^-]/[\rm \sum H_2SO_4]$\\
\code{so3} & - & ionization fraction $so_3 = [\rm SO_4^{2-}]/[\rm \sum H_2SO_4]$\\
\code{dTAdSumH2SO4}& - & $\rm \frac{\partial [TA]}{[\partial \sum H_2SO_4]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{f1} & - & ionization fraction $f_1 = [\rm HF]/[\rm \sum HF]$\\
\code{f2} & - & ionization fraction $f_1 = [\rm F^-]/[\rm \sum HF]$\\
\code{dTAdSumHF} & - & $\rm \frac{\partial [TA]}{[\partial \sum HF]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{p1} & - & ionization fraction $p_1 [\rm H_3PO_4]/[\rm \sum H_3PO_4]$\\
\code{p2} & - & ionization fraction $p_2 [\rm H_2PO_4^-]/[\rm \sum H_3PO_4]$\\
\code{p3} & - & ionization fraction $p_3 [\rm HPO_4^{2-}]/[\rm \sum H_3PO_4]$\\
\code{p4} & - & ionization fraction $p_4 [\rm PO_4^{3-}]/[\rm \sum H_3PO_4]$\\
\code{dTAdSumH3PO4} & - & $\rm \frac{\partial [TA]}{[\partial \sum H_3PO_4]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{si1} & - & ionization fraction $si_1 = [\rm Si(OH)_4] / [\rm \sum Si(OH)_4]$\\
\code{si2} & - & ionization fraction $si_2 = [\rm SiO(OH)_3] / [\rm \sum Si(OH)_4]$\\
\code{si3} & - & ionization fraction $si_3 = [\rm SiO_2(OH)_2] / [\rm \sum Si(OH)_4]$\\
\code{dTAdSumSumSiOH4} & - & $\rm \frac{\partial [TA]}{[\partial \sum Si(OH)_4]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{s1} & - & ionization fraction $s_1 = [\rm H_2S] / [\rm \sum H_2S]$\\
\code{s2} & - & ionization fraction $s_2 = [\rm HS^-] / [\rm \sum H_2S]$\\
\code{s3} & - & ionization fraction $s_3 = [\rm S^{2-}] / [\rm \sum H_2S]$ \\
& & Note that we do assume that $\rm S^{2-}$ exists. However, $s_3$ is very small.\\
\code{dTAdSumH2S} & - & $\rm \frac{\partial [TA]}{[\partial \sum H_2S]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{n1} & - & ionization fraction $n_1 = [\rm NH_4^+] / [\rm \sum NH_4^+]$\\
\code{n2} & - & ionization fraction $n_2 = [\rm NH_3] / [\rm \sum NH_4^+]$\\
\code{dTAdSumNH4} & - & $\rm \frac{\partial [TA]}{[\partial \sum NH_4^+]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{na1} & - & ionization fraction $na_1 = [\rm HNO_3] / [\sum HNO_3]$\\
\code{na2} & - & ionization fraction $na_2 = [\rm NO_3^-] / [\sum HNO_3]$\\
\code{dTAdSumHNO3} & - & $\rm \frac{\partial [TA]}{[\partial \sum HNO_3]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{ni1} & - & ionization fraction $ni_1 = [\rm [HNO_2] / [\rm \sum HNO_2]$\\
\code{ni2} & - & ionization fraction $ni_2 = [\rm [NO_2^-] / [\rm \sum HNO_2]$\\
\code{dTAdSumHNO2} & - & $\rm \frac{\partial [TA]}{[\partial \sum HNO_2]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{dTAdH} & - & $\rm \frac{\partial [TA]}{[\partial [H^+]}$: buffer factor\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ...)$\\
\code{dTAdKdKdS} & - & $\rm \sum_i \frac{\partial [TA]}{\partial K^*_i} \frac{\partial K^*_i}{\partial S}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ..., K^*_i)$\\
\code{dTAdKdKdT} & - & $\rm \sum_i \frac{\partial [TA]}{\partial K^*_i} \frac{\partial K^*_i}{\partial T}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ..., K^*_i)$\\
\code{dTAdKdKdp} & - & $\rm \sum_i \frac{\partial [TA]}{\partial K^*_i} \frac{\partial K^*_i}{\partial p}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ..., K^*_i)$\\
\code{dTAdKdKdSumH2SO4} & - & $\rm \sum_i \frac{\partial [TA]}{\partial K^*_i} \frac{\partial K^*_i}{\partial [\sum H_2SO_4]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ..., K^*_i)$\\
\code{dTAdKdKdSumHF} & - & $\rm \sum_i \frac{\partial [TA]}{\partial K^*_i} \frac{\partial K^*_i}{\partial [\sum HF]}$\\
& & with $\rm [TA] = f([H^+], [\sum CO_2], ..., K^*_i)$\\ \hline
\end{longtable}
\end{footnotesize}
\section[Using AquaEnv]{Using \textbf{\textsf{AquaEnv}}}
\subsection[Calling the K functions directly]{Calling the ``K'' functions directly}
The elements \code{K\_CO2, K\_HCO3, K\_BOH3, K\_W, K\_HSO4, K\_HF, K\_NH4, K\_H2S, K\_H3PO4,}\\
\code{K\_H2PO4, K\_HPO4, K\_SiOH4, K\_SiOOH3, K0\_CO2, K0\_O2, Ksp\_aragonite,} and\\
\code{Ksp\_calcite} can be calculated directly. This is done via functions that bear the same name as those elements.
<<CallingKFuncsDirectly, fig=FALSE, echo=TRUE>>=
K_CO2(S = 15, t = 30)
K0_CO2(S = 15, t = 30)
Ksp_calcite(S = 15, t = 30, p = 100)
@
One of the input arguments can also be a vector:
<<Avector, fig=FALSE, echo=TRUE>>=
K_CO2(S = 10:15, t = 30)
@
\subsection[The function aquaenv]{The function \code{aquaenv}}
\subsubsection[Minimal aquaenv definition]{Minimal \code{aquaenv} definition}
Minimally, an object of class \code{aquaenv} can be defined with just a temperature (\code{t}) and salinity (\code{S}) value
<<MinimalAquaenv, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15)
ae$K_CO2
@
Optionally the pressure (here the gauge pressure p) can be given. As in the above case, the returned object of class \code{aquaenv} then contains
a standard set of elements as shown by the \code{names} command.
<<AquaenvWithDepth, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15, p = 10)
names(ae)
ae$Ksp_calcite
@
The pressure can also be given via the total pressure (\code{P}) or the water depth (\code{d}). The atmospheric pressure can be given as well (\code{Pa}, e.g. for the case of a mountain lake).
Furthermore, if the depth is given, the latitude (\code{lat}) can also be specified (default is 0 degrees).
<<AquaenvPressures, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15, p = 10)
ae[c("p", "P")]
ae <- aquaenv(S = 30, t = 15, P = 10)
unlist(ae[c("p", "P")])
ae <- aquaenv(S = 30, t = 15, P = 10, Pa = 0.5)
unlist(ae[c("p", "P")])
ae <- aquaenv(S = 30, t = 15, d = 100)
unlist(ae[c("p", "P")])
ae <- aquaenv(S = 30, t = 15, d = 100, lat = 51)
unlist(ae[c("p", "P")])
@
A minimal set of elements in an object of class \code{aquaenv} can be obtained by setting the flag \code{skeleton} to \code{TRUE}.
<<SkeletonAquaenv, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15, p = 10, skeleton = TRUE)
names(ae)
@
\subsubsection[Defining the complete aquaenv system in different ways]{Defining the complete \code{aquaenv} system in different ways}
If enough information is given to define a complete speciation, i.e. either one of the pairs
$[\rm \sum CO_2]$ and pH, $[\rm \sum CO_2]$ and $\rm [TA]$, $[\rm \sum CO_2]$ and $[\rm CO_2]$, or $[\rm \sum CO_2]$ and $\rm fCO_2$,
a full \code{aquaenv} system can be defined.
<<CompleteAquaenvSystem, fig=FALSE, echo=TRUE>>=
S <- 30
t <- 15
p <- 10
SumCO2 <- 0.0020
pH <- 8
TA <- 0.002142233
fCO2 <- 0.0005272996
CO2 <- 2.031241e-05
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, pH = pH)
ae$TA
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, TA = TA)
ae$pH
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, CO2 = CO2)
ae$pH
names(ae)
@
As seen above, a full speciation is calculated along with the pH or total alkalinity, respectively. If only pH or total alkalinity is needed,
the calculation of the full speciation can be toggled off. Furthermore, the flag \code{skeleton} also works for a full system.
<<SpeciationSkeletonAquaenv, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, pH = pH, speciation = FALSE)
names(ae)
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, pH = pH, speciation = FALSE,
skeleton = TRUE)
names(ae)
@
All the quantities needed for the explicit pH modelling approaches as given in \cite{Hofmann2008} and \cite{Hofmann2010a} can be calculated by
setting the flag \code{dsa} to \code{TRUE}.The flag \code{revelle} is redundant; the Revelle factor can instead be calculated analytically using \code{BufferFactors$RF} following \cite{Hagens2016}. This is further detailed in section XX.
<<DSAAquaenv, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, fCO2 = fCO2, dsa = TRUE)
ae$dTAdH
ae$revelle
@
If an ambivalent situation is created because the user enters too much information, an error message is displayed
\begin{scriptsize}
<<ErrorMessagesAquaenv, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, CO2 = CO2, fCO2 = fCO2)
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, pH = pH, TA = TA)
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, pH = pH, CO2 = CO2)
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, pH = pH, fCO2 = fCO2)
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, TA = TA, CO2 = CO2)
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, TA = TA, fCO2 = fCO2)
@
\end{scriptsize}
\subsubsection[Calculating SumCO2]{Calculating $\rm [\sum CO_2]$}
$\rm [\sum CO_2]$ can be calculated by giving a constant pair of either pH and $[\rm CO_2]$, pH and $\rm fCO_2$, pH and $\rm [TA]$, $\rm [TA]$ and $[\rm CO_2]$, or $\rm [TA]$ and $\rm fCO_2$
<<SumCO2Calculating, fig=FALSE, echo=TRUE>>=
fCO2 <- 0.0006943363
CO2 <- 2.674693e-05
pH <- 7.884892
TA <- 0.0021
S <- 30
t <- 15
p <- 10
ae <- aquaenv(S, t, p, SumCO2 = NULL, pH = pH, CO2 = CO2)
ae$SumCO2
ae <- aquaenv(S, t, p, SumCO2 = NULL, pH = pH, fCO2 = fCO2)
ae$SumCO2
ae <- aquaenv(S, t, p, SumCO2 = NULL, pH = pH, TA = TA)
ae$SumCO2
ae <- aquaenv(S, t, p, SumCO2 = NULL, TA = TA, CO2 = CO2)
ae$SumCO2
ae <- aquaenv(S, t, p, SumCO2 = NULL, TA = TA, fCO2 = fCO2)
ae$SumCO2
@
\subsubsection[Cloning an object of class aquaenv]{Cloning an object of class \code{aquaenv}}
It is possible to clone an object of class \code{aquaenv}, either 1 to 1 or with different pH, $\rm [TA]$, or $\rm K^*_{CO_2}$.
<<CloningAquaenv, fig=FALSE, echo=TRUE>>=
S <- 30
t <- 15
SumCO2 <- 0.0020
TA <- 0.00214
ae <- aquaenv(S, t, SumCO2 = SumCO2, TA = TA)
ae$pH
ae1 <- aquaenv(ae = ae) # this is the same
ae1$pH
ae2 <- aquaenv(ae = ae, pH = 9)
c(ae$TA, ae2$TA)
ae3 <- aquaenv(ae = ae, TA = 0.002)
c(ae$pH, ae3$pH)
K_CO2 <- 1e-6
ae4 <- aquaenv(ae = ae, k_co2 = 1e-6)
c(ae$TA, ae4$TA)
@
Note that \code{k\_co2} as an input argument is in lower case (in constrast to the element \code{K\_CO2} of class \code{aquaenv}!
\subsubsection{Preparing input arguments}
Input arguments for the function \code{aquaenv} need to be in mol/kg-solution and on the free pH scale.
Data in other concentration units or pH scales can be converted using the function \code{convert}. This is a rather complex function, but its typical use is:
\begin{verbatim}
convert(x, vartype, what, S, t, p, SumH2SO4, SumHF, khf)
convert(x, from, to, factor, convattr)
\end{verbatim}
<<PreparingInput, fig=FALSE, echo=TRUE>>=
S <- 10
t <- 15
pH_NBS <- 8.142777
SumCO2molar <- 0.002016803
(pH_free <- convert(pH_NBS, "pHscale", "nbs2free", S = S, t = t))
(SumCO2molin <- convert(SumCO2molar, "conc", "molar2molin", S = S, t = t))
ae <- aquaenv(S, t, SumCO2 = SumCO2molin, pH = pH_free)
ae$pH
ae$SumCO2
@
\subsubsection{Vectors as input variables}
One of the input variables for the function \code{aquaenv} may be a vector. All the other input variables are then assumed to be constant.
The elements of the resulting two dimensional object of class \code{aquaenv} are then vectors containing the elements as a function of the input variable
for which a vector is given.
<<InputVectors, fig=FALSE, echo=TRUE>>=
SumCO2 <- 0.0020
pH <- 8
S <- 30
t <- 1:5
p <- 10
ae <- aquaenv(S, t, p, SumCO2 = SumCO2, pH = pH)
rbind(t = ae$t, TA = ae$TA)
@
A two dimensional object of class \code{aquaenv} can be visualized using the \code{plot} function.
For convenience of the user, the default setting for the \code{plot} function for an object of class \code{aquaenv} results in a
new plotting device being opened. Setting the flag \code{newdevice} to \code{FALSE} prevents that.
The \code{plot} function plots all elements of the respective object of class \code{aquaenv}. This, however, might not be what the user wants,
especially if a larger plotting device cannot properly displayed like in the case above. In this case the parameter \code{what} can be used.
Note, however, that the default setting for calling \code{plot} with the parameter \code{what} is that \code{mfrow=c(1,1)}. So
if one wants to plot several elements, \code{mfrow} needs to be set to a suitable value.
<<Selectplot, fig=TRUE, echo=TRUE, width=10, height=2>>=
plot(ae, xval = t, xlab = "t/(deg C)",
what = c("pH", "CO2", "HCO3", "CO3"),
mfrow = c(1, 4))
@
Other input arguments can also be vectors, e.g. the obvious S, t, or p:
<<MoreInputVectors1, fig=FALSE, echo=TRUE, eval=FALSE>>=
ae <- aquaenv(S=20:24, t=15, p=10, SumCO2 = SumCO2, pH = pH, dsa = TRUE)
rbind(ae$S, ae$TA)
@
But also $\rm [\sum CO_2]$, $\rm [TA]$, pH and $[\rm \sum NH_4^+]$ can be vectors, e.g.
<<MoreInputVectors6, echo=TRUE, width=10, height=2>>=
ae <- aquaenv(20, 10, SumCO2=seq(0.001, 0.002, 0.00025), TA = 0.002)
rbind(ae$SumCO2, ae$pH, ae$HCO3)
@
\subsubsection[Calculating SumCO2 from input vectors]{Calculating $\rm [\sum CO_2]$ from input vectors}
The functionality of calculating $\rm [\sum CO_2]$ can also be used together with vectors as input variables.
<<SumCO2CalcInputVecs1, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 11:15, SumCO2 = NULL, pH = pH, CO2 = CO2,
dsa = TRUE)
ae$SumCO2
@
\subsubsection{Conversion from and to a dataframe}
Objects of class \code{aquaenv} can be converted to an \R \code{data.frame} to further post-process them with standard
\R means. Similarly, \R \code{data.frames} can be converted to objects of class \code{aquaenv} to use the plotting
facilities that exist for objects of class \code{aquaenv}. This can be helpful for plotting output of a dynamic model run, e.g. from \R package \ds,
as will be shown later in this document.
<<Dataframe, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 11:15, SumCO2 = NULL, pH = 8, CO2 = 2e-5)
aedataframe <- as.data.frame(ae)
dim(aedataframe)
aedataframe[, 1:3]
aetest <- aquaenv(ae = aedataframe, from.data.frame = TRUE)
@
\subsubsection[Converting elements in an object of class aquaenv]{Converting elements in an object of class \code{aquaenv}}
Elements of an object of class \code{aquaenv} are calculated in predefined units, e.g.,
the concentrations are in unit mol/kg-solution (molinity).
The function \code{convert} can be used to convert all elements in an object of class \code{aquaenv} that share a common attribute,
e.g. the unit.
<<Elementconversion, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 10)
ae$SumBOH3
ae <- convert(ae, "mol/kg-soln", "umol/kg-H2O", 1e6/ae$molal2molin, "unit")
ae$SumBOH3
@
\subsubsection{Quantities needed for explicit pH modelling}
As mentioned above, the quantities needed for the explicit pH modelling approach (direct substitution approach - DSA) as presented by \cite{Hofmann2008}
can be calculated with the function \code{aquaenv} by setting the flag \code{dsa} to \code{TRUE}.
<<DSAQuantities1, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15, d = 10, SumCO2 = 0.002, pH = 8, dsa = TRUE)
@
This command calculates the buffer factor and the partial derivatives of \code{[TA]} with respect to other summed quantities referred to in \cite{Hofmann2008}. This functionality is expanded with the function \code{BufferFactors}, which is discussed in a separate subsection.
<<DSAQuantities2, fig=FALSE, echo=TRUE>>=
ae$dTAdH
ae$dTAdSumCO2
@
Moreover, setting the flag \code{dsa} to \code{TRUE} calculates the sums of the partial derivatives of [TA] with respect to the equilibrium constants ($K^*$'s) multiplied with the partial derivatives of the respective
equilibrium constant with one of their variables (i.e., S, T, d, $[\rm \sum H_2SO_4]$, or $[\rm \sum HF]$) as introduced in \cite{Hofmann2009}.
<<DSAQuantities3, fig=FALSE, echo=TRUE>>=
ae$dTAdKdKdS
ae$dTAdKdKdSumH2SO4
@
Furthermore, the ionization fractions used for the pH dependent fractional stoichiometric pH modelling approach described in \cite{Hofmann2010a} are calculated as well
<<DSAQuantities4, fig=FALSE, echo=TRUE>>=
ae$c1
@
\subsection[The use of BufferFactors]{The use of \code{BufferFactors}}
\subsubsection[An introduction to BufferFactors]{An introduction to \code{BufferFactors}}
As described above, by setting the flag \code{dsa} to \code{TRUE}, the function \code{aquaenv} calculates the total buffer factor and the partial derivatives of \code{[TA]} with respect to other summed quantities referred to in \cite{Hofmann2008}. The function \code{BufferFactors} allows for a more generic, analytical calculation of partial derivatives related to carbonate system calculations, as presented in \citet{Hagens2016}. This includes the sensitivity of pH and concentrations of $[\rm CO_2]$ and other acid-base species to a change in ocean chemistry, i.e. a change in $[\rm TA]$ or other summed quantities.
\code{BufferFactors} internally calls the \code{aquaenv} function and uses its output to calculate the sensitivities. Its output is a list of elements, the first of which is of class \code{aquaenv} and is the output of the internal call to \code{aquaenv}. \code{BufferFactors} runs without specifying any input variables; in this case, the global contemporary surface ocean values given in Table 4 of \citet{Hagens2016} are used as input.
<<BufferFactors1, fig=FALSE, echo=TRUE>>=
BF <- BufferFactors()
names(BF)
BF$dtotX.dpH
@
\subsubsection[Specifying input parameters for BufferFactors]{Specifying input parameters for \code{BufferFactors}}
An object of class \code{aquaenv} where a complete speciation is calculated can be used as input for \code{BufferFactors}:
<<BufferFactors2, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15, d = 10, SumCO2 = 0.002, pH = 8.1,
skeleton = TRUE)
BF <- BufferFactors(ae = ae)
BF$RF
@
Note that \code{ae$dTAdH} and \code{BF$beta.H} both represent the total buffer factor and are thus the same:
<<BufferFactors3, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15, d = 10, SumCO2 = 0.002, pH = 8.1, dsa = TRUE)
BF <- BufferFactors(ae = ae)
cbind(ae$dTAdH,BF$beta.H)
@
As an alternative to an object of class \code{aquaenv}, input parameters can be given by specifying the argument \code{parameters}. This argument is a vector which can contain the following inputs: \code{DIC}, \code{TotNH3}, \code{TotP}, \code{TotNO3}, \code{TotNO2}, \code{TotS}, \code{TotSi}, \code{TB}, \code{TotF}, \code{TotSO4}, \code{sal}, \code{temp}, \code{pres}, \code{Alk}. All of these should be supplied in the same units as for \code{aquaenv}, i.e. mol/kg-soln.
<<BufferFactors4, fig=FALSE, echo=TRUE>>=
parameters <- c(DIC = 0.002, Alk = 0.0022)
BF <- BufferFactors(parameters = parameters)
BF$RF
@
Note that if temperature, salinity, pressure or a specific total concentration is not given in \code{parameters}, the default value of the contemporary global ocean is used. However, if an object of class \code{aquaenv} is provided, all of its total concentrations are used, even if they are not specified in the calculation of the object, as is the case for e.g. $\rm [\sum NH_4^+]$ in the example above. However, one or more output variables of \code{aquenv} can be overruled by specifying them in the \code{parameters} argument, and specifying both \code{aquaenv} and \code{parameters} as input parameters for \code{BufferFactors}:
<<BufferFactors5, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15, d = 10, SumCO2 = 0.002, pH = 8.1,
skeleton = TRUE)
BF <- BufferFactors(ae = ae)
BF$RF
parameters <- c(Alk = 0.0022)
BF_2 <- BufferFactors(ae = ae, parameters = parameters)
BF_2$RF
@
\subsubsection[Specifying which output is produced by BufferFactors]{Specifying which output is produced by \code{BufferFactors}}
By default, \code{BufferFactors} only displays the sensitivities related to $\rm [\sum CO_2]$. This is indicated by the argument \code{species}, which defaults to \code{species = c("SumCO2")}. However, sensitivities can be calculkated for any total concentration contributing to $\rm [TA]$, or any acid-base species contributing to this total concentration:
<<BufferFactors6, fig=FALSE, echo=TRUE>>=
BF <- BufferFactors(species = c("CO2", "HCO3", "CO3", "SumNH4"))
BF$dtotX.dX
@
In the case when species are defined which corresponding total concentration equals zero, the corresponding output produces \code{NaN}:
<<BufferFactors7, fig=FALSE, echo=TRUE>>=
ae <- aquaenv(S = 30, t = 15, d = 10, SumCO2 = 0.002, pH = 8.1,
skeleton = TRUE)
BF <- BufferFactors(ae = ae, species = c("SumCO2", "SumNH4"))
BF$dtotX.dX
@
\subsubsection[Other arguments of BufferFactors]{Other arguments of \code{BufferFactors}}
\code{BufferFactors} allows using different pre-defined sets of equilibrium constants, as well as specifying specific values for them, similarly to how this is done in the function \code{aquaenv}:
<<BufferFactors7, fig=FALSE, echo=TRUE>>=
BF <- BufferFactors(k1k2 = "roy")
BF$RF
BF <- BufferFactors(k_co2 = 1e-6)
BF$RF
@
\subsection[The plot.aquaenv function]{The \code{plot.aquaenv} function}
In the previous sections, the \code{plot} function has been introduced.
When the first element of the arguments list of
\code{plot} is an object of class \code{aquaenv}, the function \code{plot.aquaenv} is called.
This is a multifunctional tool to visualize information contained in an object of class
\code{aquaenv}.
For the convenience of the users, \code{plot.aquaenv} combines the call of standard \R plotting functions and the previous call of the function \code{par}
to set parameters like \code{mfrow}, \code{mar}, etc. as well as the opening of a plotting device with a certain size. As already shown above, setting the flag \code{newdevice} to \code{FALSE} suppresses the opening of a new plotting device.\\
\noindent
For example
<<Plot1, fig=TRUE, echo=TRUE, width=10, height=2>>=
ae <- aquaenv(20:30, 10)
plot(ae, xval = 20:30, xlab = "S", what = c("K_CO2", "K_HCO3", "K_BOH3"),
size = c(10, 2), mfrow = c(1,3))
@
Furthermore the parameter \code{device} can be specified which allows the user to write the plots to .eps and .pdf files.
The parameter \code{filename} can be used to specify a filename other than the default filename ``aquaenv''.
Whereas these features make the function \code{plot.aquaenv} different from standard \R plotting functions,
when setting the flags \code{newdevice} and \code{setpar} to \code{FALSE}, \code{plot.aquaenv} behaves like a ``normal'' \R plotting function.
Furthermore, the function \code{plot.aquaenv} can be used to create ``cumulative'' plots and ``Bjerrum'' plots. This will be explained in upcoming sections.
\subsection[Using objects of class aquaenv in dynamic models]{Using objects of class \code{aquaenv} in dynamic models}
\subsubsection{Ordinary dynamic equation models}
It is convenient to use objects of class \code{aquaenv} in a dynamic model, e.g. solved using the \R package \ds.
This is illustrated with an example\footnote{For information about how to set up a dynamic model with \ds, consult the documentation of \ds}.\\
\noindent
First, the input parameters are specified:
\begin{scriptsize}
<<OrdDynModParamList, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
parameters <- list(
S = 25 , # psu
t_min = 5 , # degrees C
t_max = 25 , # degrees C
k = 0.4 , # 1/d proportionality factor for air-water exchange
rOx = 0.0000003 , # mol-N/(kg*d) maximal rate of oxic mineralisation
rNitri = 0.0000002 , # mol-N/(kg*d) maximal rate of nitrification
rPP = 0.000006 , # mol-N/(kg*d) maximal rate of primary production
ksDINPP = 0.000001 , # mol-N/kg
ksNH4PP = 0.000001 , # mol-N/kg
D = 0.1 , # 1/d (dispersive) transport coefficient
O2_io = 0.000296 , # mol/kg-soln
NO3_io = 0.000035 , # mol/kg-soln
SumNH4_io = 0.000008 , # mol/kg-soln
SumCO2_io = 0.002320 , # mol/kg-soln
TA_io = 0.002435 , # mol/kg-soln
C_Nratio = 8 , # mol C/mol N C:N ratio of organic matter
a = 30 , # time at which PP begins
b = 50 , # time at which PP stops again
modeltime = 100 # duration of the model
)
@
\end{scriptsize}
Next the model function that calculates the derivatives is defined. In this function, an object of class \code{aquaenv} is created in each timestep, some of its elements are used to calculated kinetic rate expressions and part of the object is returned as output. A function that returns the temperature is defined first.
\begin{scriptsize}
<<OrdDynModFunction, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
temperature <- with (parameters,
approxfun(x = 0:101,
y = c(seq(t_min, t_max, (t_max-t_min)/50),
seq(t_max, t_min, -(t_max-t_min)/50)))
)
boxmodel <- function(time, state, parameters) {
with (
as.list(c(state, parameters)), {
t <- temperature(time)
ae <- aquaenv(S = S, t = t, SumCO2 = SumCO2,
SumNH4 = SumNH4, TA = TA)
ECO2 <- k * (ae$CO2_sat - ae$CO2)
EO2 <- k * (ae$O2_sat - O2)
# dilution
TO2 <- D*(O2_io - O2)
TNO3 <- D*(NO3_io - NO3)
TSumNH4 <- D*(SumNH4_io - SumNH4)
TTA <- D*(TA_io - TA)
TSumCO2 <- D*(SumCO2_io - SumCO2)
RNit <- rNitri * SumNH4/(SumNH4+1e-8)
ROx <- rOx
ROxCarbon <- ROx * C_Nratio
pNH4PP <- 0
RPP <- 0
if ((time > a) && (time < b)) {
RPP <- rPP * ((SumNH4+NO3)/(ksDINPP + (SumNH4+NO3)))
pNH4PP <- SumNH4/(SumNH4+NO3)
}
RPPCarbon <- RPP * C_Nratio
dO2 <- TO2 + EO2 - ROxCarbon - 2*RNit + (2-2*pNH4PP)*RPP + RPPCarbon
dNO3 <- TNO3 + RNit -(1-pNH4PP)*RPP
dSumCO2 <- TSumCO2 + ECO2 + ROxCarbon - RPPCarbon
dSumNH4 <- TSumNH4 + ROx - RNit - pNH4PP*RPP
dTA <- TTA + ROx - 2*RNit -(2*pNH4PP-1)*RPP
ratesofchanges <- c(dO2, dNO3, dSumNH4, dSumCO2, dTA)
return(list(ratesofchanges, ae[c("t", "NH4", "NH3", "pH")]))
} )
}
@
\end{scriptsize}
The model is solved, and the output can be plotted in the same way as any \code{deSolve} object.
The package \pkg{deSolve} \citep{deSolve} has to be loaded first.
\begin{scriptsize}
<<OrdDynModSolution, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
initialstate <- with (parameters,
c(O2=O2_io, NO3=NO3_io, SumNH4=SumNH4_io, SumCO2=SumCO2_io, TA=TA_io))
times <- 1:100
output <- vode(initialstate,times,boxmodel,parameters, hmax = 1)
@
\end{scriptsize}
<<OrdDynModePlotting, fig = TRUE, echo=TRUE, eval=FALSE>>=
plot(output)
@
\subsubsection{Models using the explicit pH modelling approach}
\paragraph{In one single model} $\;$\\
Since an object of class \code{aquaenv} can contain all quantities necessary to employ the explicit pH modelling approaches as introduced by
\cite{Hofmann2008, Hofmann2009, Hofmann2010a}, they can readily be used in an explicit pH model.\\
\noindent
As an example, we give a model that calculates the pH in the ``classical'' way in every timestep using \code{aquaenv}, also employs the explicit pH modelling
approach (direct substitution approach - DSA) given in \cite{Hofmann2008} and additionally employs fractional stoichiometry as given in \cite{Hofmann2010a}.
The pH evolution is thus calculated in three different ways which allows comparing the three values for consistency.\\
\noindent
Again, a list of parameters is defined
\begin{scriptsize}
<<SinglePHModParams, fig=FALSE, echo=TRUE, keep.source=TRUE>>=
parameters <- list(
S = 25 , # psu
t = 15 , # degrees C
k = 0.4 , # 1/d proportionality factor for air-water exchange
rOx = 0.0000003 , # mol-N/(kg*d) maximal rate of oxic mineralisation
rNitri = 0.0000002 , # mol-N/(kg*d) maximal rate of nitrification
rPP = 0.0000006 , # mol-N/(kg*d) maximal rate of primary production
ksSumNH4 = 0.000001 , # mol-N/kg
D = 0.1 , # 1/d (dispersive) transport coefficient
SumNH4_io = 0.000008 , # mol/kg-soln
SumCO2_io = 0.002320 , # mol/kg-soln
TA_io = 0.002435 , # mol/kg-soln
C_Nratio = 8 , # mol C/mol N C:N ratio of organic matter
a = 5 , # time from which PP begins
b = 20 , # time where PP shuts off again
modeltime = 30 # duration of the model
)
@
\end{scriptsize}
and a model function is defined. Again, an object of class \code{aquaenv} is created in each timestep and the respective elements are used.
\begin{scriptsize}
<<SinglePHModFunction, fig=FALSE, echo=TRUE, keep.source=TRUE>>=
boxmodel <- function(timestep, currentstate, parameters)
{
with (
as.list(c(currentstate,parameters)),
{
# the "classical" implicit pH calculation method is applied in aquaenv
ae <- aquaenv(S=S, t=t, SumCO2=sumCO2,
SumNH4=sumNH4, TA=alkalinity, dsa=TRUE)
ECO2 <- k * (ae$CO2_sat - ae$CO2)
RNit <- rNitri
ROx <- rOx
if ((timestep > a) && (timestep < b))
RPP <- rPP * (sumNH4/(ksSumNH4 + sumNH4))
else
RPP <- 0
dsumCO2 <- ECO2 + C_Nratio*ROx - C_Nratio*RPP
dsumNH4 <- ROx - RNit - RPP
dalkalinity <- ROx - 2*RNit - RPP
# The DSA pH
dH <- (dalkalinity - (dsumCO2*ae$dTAdSumCO2 + dsumNH4*ae$dTAdSumNH4))/ae$dTAdH
DSApH <- -log10(H)
# The DSA pH using pH dependent fractional stoichiometry (= using partitioning coefficients)
rhoHECO2 <- ae$c2 + 2*ae$c3
rhoHRNit <- 1 + ae$n1
rhoHROx <- C_Nratio * (ae$c2 + 2*ae$c3) - ae$n1
rhoHRPP <- -(C_Nratio * (ae$c2 + 2*ae$c3)) + ae$n1
dH_ECO2 <- rhoHECO2*ECO2/(-ae$dTAdH)
dH_RNit <- rhoHRNit*RNit/(-ae$dTAdH)
dH_ROx <- rhoHROx*ROx /(-ae$dTAdH)
dH_RPP <- rhoHRPP*RPP /(-ae$dTAdH)
dH_stoich <- dH_ECO2 + dH_RNit + dH_ROx + dH_RPP
DSAstoichpH <- -log10(H_stoich)
ratesofchanges <- c(dsumNH4, dsumCO2, dalkalinity, dH, dH_stoich)
processrates <- c(ECO2=ECO2, RNit=RNit, ROx=ROx, RPP=RPP)
DSA <- c(DSApH=DSApH, rhoHECO2=rhoHECO2, rhoHRNit=rhoHRNit, rhoHROx=rhoHROx,
rhoHRPP=rhoHRPP, dH_ECO2=dH_ECO2, dH_RNit=dH_RNit, dH_ROx=dH_ROx,
dH_RPP=dH_RPP, DSAstoichpH=DSAstoichpH)
return(list(ratesofchanges, processrates, DSA, ae))
}
)
}
@
\end{scriptsize}
The model is solved
\begin{scriptsize}
<<SinglPHModSolution, fig=FALSE, echo=TRUE, keep.source=TRUE>>=
with (as.list(parameters),
{
H_init <<- 10^(-(aquaenv(S=S, t=t, SumCO2=SumCO2_io, SumNH4=SumNH4_io, TA=TA_io,
speciation=FALSE)$pH))
initialstate <<- c(sumNH4=SumNH4_io, sumCO2=SumCO2_io, alkalinity =TA_io, H=H_init,
H_stoich=H_init)
times <<- c(0:modeltime)
})
output <- vode(initialstate, times, boxmodel, parameters, hmax=1)
@
\end{scriptsize}
and output can be plotted, again using \code{plot.aquaenv}. Note that here the parameter \code{what} is used.
<<SinglePHModePlotting1, fig=TRUE, echo=TRUE, width=10, height=7>>=
select <- c("sumCO2", "alkalinity", "sumNH4", "RPP", "dTAdH", "dTAdSumCO2",
"dTAdSumNH4","rhoHECO2", "rhoHRNit", "rhoHROx","dH_ECO2","dH_RNit",
"dH_RPP","pH", "DSApH", "DSAstoichpH")
plot(output, which = select, xlab="time/d", mfrow=c(4,4))
@
Here, the cumulative plotting functionality of \code{plot.aquaenv} can be employed as well to visualize the influences of the different kinetically
modelled processes on $\rm [H^+]$.
<<SinglepHModPlotting2, fig=TRUE, echo=TRUE, width=10, height=5>>=
what <- c("dH_ECO2", "dH_RNit", "dH_ROx", "dH_RPP")
plot(aquaenv(ae=as.data.frame(output), from.data.frame=TRUE),
xval = times, what = what, xlab = "time/d", size = c(7,5),
ylab = "mol-H/(kg-soln*d)", legendposition = "bottomright", cumulative = TRUE)
@
Finally, the pH values calculated with the three different methods can be plotted in one single graph to see that they are identical, i.e. the three methods of pH calculation are consistent with each other
<<SinplePHModConsistencyCheck, fig=TRUE, echo=TRUE, width=10, height=5>>=
matplot(output[, "time"],
output[, c("pH", "DSApH", "DSAstoichpH")], type = "l", lty = 1, lwd = 2)
@
\paragraph{In three separate models} $\;$\\
\subparagraph{The implicit pH modelling approach}$\,$\\
A list of parameters:
\begin{scriptsize}
<<ImplicitPHModParams, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
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
)
@
\end{scriptsize}
The model function:
\begin{scriptsize}
<<ImplicitPHModFunction, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
boxmodel <- function(timestep, currentstate, parameters)
{
with (
as.list(c(currentstate,parameters)),
{
ae <- aquaenv(S=S, t=t, SumCO2=SumCO2, TA=TA, SumSiOH4=0,
SumBOH3=0, SumH2SO4=0, SumHF=0)
Rc <- kc * ((ae$CO2_sat) - (ae$CO2))
Rp <- kp * (1-ae$omega_calcite)^n
dSumCO2 <- Rc - Rp
dTA <- -2*Rp
ratesofchanges <- c(dSumCO2, dTA)
processrates <- c(Rc=Rc, Rp=Rp)
return(list(ratesofchanges, list(processrates, ae)))
}
)
}
@
\end{scriptsize}
Solving the model:
\begin{scriptsize}
<<ImplicitPHModSolution, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
with (as.list(parameters),
{
initialstate <<- c(SumCO2=SumCO2_t0, TA=TA_t0)
times <<- seq(0,modeltime,(modeltime/outputsteps))
})
output <<- vode(initialstate,times,boxmodel,parameters, hmax=1)
@
\end{scriptsize}
Visualisation of the results can be done similarly as above.
\subparagraph{The explicit pH modelling approach}$\,$\\
A list of parameters:
\begin{scriptsize}
<<ExplicitPHModParams, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
parameters <- list(
S = 35 , # psu
t = 15 , # degrees C
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
)
@
\end{scriptsize}
The model function:
\begin{scriptsize}
<<ExplicitPHModFunction, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
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)))
}
)
}
@
\end{scriptsize}
Solving the model:
\begin{scriptsize}
<<ExplicitPHModSolution, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
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 <- vode(initialstate,times,boxmodel,parameters, hmax=1)
@
\end{scriptsize}
Visualisation of the results can again be done similarly as above.
\subparagraph{The fractional stoichiometric approach} $\;$\\
A list of parameters:
\begin{scriptsize}
<<FracStoichModParams, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
parameters <- list(
S = 35 , # psu
t = 15 , # degrees C
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
)
@
\end{scriptsize}
The model function:
\begin{scriptsize}
<<FracStoichModFunction, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
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
rhoc <- ae$c2 + 2*ae$c3
rhop <- 2*ae$c1 + ae$c2
dHRc <- rhoc*Rc/(-ae$dTAdH)
dHRp <- rhop*Rp/(-ae$dTAdH)
dH <- dHRc + dHRp
ratesofchanges <- c(dSumCO2, dH)
processrates <- c(Rc=Rc, Rp=Rp)
outputvars <- c(dHRc=dHRc, dHRp=dHRp, rhoc=rhoc, rhop=rhop)
return(list(ratesofchanges, list(processrates, outputvars, ae)))
}
)
}
@
\end{scriptsize}
Solving the model:
\begin{scriptsize}
<<FracStoichModSolution, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
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))
})
@
\end{scriptsize}
Visualisation of the results can again be done similarly as above.
\subsection[Titration simulation: the function titration]{Titration simulation: the function \code{titration}}
With the function \code{titration} \aq$\,$ provides a powerful tool to simulate titrations.
A two dimensional object of class \code{aquaenv} will be created where the second dimension is the amount of titrant added.
For this purpose, three extra elements are added to the \code{aquaenv} object that will be created:\\
\noindent
\begin{tabular}{llp{.6\textwidth}} \hline
\textbf{element} & \textbf{unit} & \textbf{explanation}\\
delta\_conc\_titrant & mol/kg-solution & the offset in concentration of the titrant that is caused by adding the titrant to the sample\\
delta\_mass\_titrant & kg & the amount of mass of titrant solution added\\
delta\_moles\_titrant & mol & the amount of moles of titrant added\\
\end{tabular}\\
\noindent
Each one of these elements is a suitable \code{xval} for plotting an \code{aquaenv} object generated by \code{titration}.
\subsubsection{Titration with HCl}
The standard titration type is titration with hydrochloric acid ($\rm HCl$). A simple example will illustrate this.\\
\noindent
An object of class \code{aquaenv} needs to be created to define the initial conditions of the titration. That is
temperature, salinity, depth, the concentrations of all summed quantities and the initial pH (or $\rm [TA]$).
<<HClTit1, fig=FALSE, echo=TRUE>>=
ae_init <- aquaenv(S = 35, t = 15, SumCO2 = 0.0035, SumNH4 = 0.00002, pH = 11.3)
@
Then \code{titration} can be run to create the object describing the simulated titration.
In this example the titrant is $\rm HCl$ of the relatively low concentration of 0.01 mol/kg-solution.
The sample solution amounts to 10 g. To sweep a considerable pH range quite a lot of titrant needs to be
added: 20 g. This means the salinity of the solution in the titration vessel will change due to dilution with the titrant solution.
For this reason, the salinity of the titrant solution needs to be given via the parameter \code{S\_titrant}.
However, we assume the titrant does not contain borate, sulfate or fluoride, that is why we do not set the
flag \code{seawater\_titrant} to \code{TRUE}.
<<HClTit2, fig=FALSE, echo=TRUE>>=
ae <- titration(ae_init, mass_sample = 0.01, mass_titrant = 0.02,
conc_titrant = 0.01, S_titrant = 0.5, steps = 100)
@
To get a quick overview, all elements of the obtained \code{aquaenv} object can be plotted (not shown)
but also a selection of elements can be plotted as a function of the added titrant mass.
<<HClTit3, fig=TRUE, echo=TRUE, width=10, height=7>>=
what <- c("TA", "pH", "CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH",
"NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F", "fCO2")
plot(ae, xval = ae$delta_mass_titrant, xlab = "HCl solution added [kg]",
what = what, size = c(12,8), mfrow = c(4,4))
@
It can also be plotted against titrant concentration, the moles of added titrant (not shown), or against other variables such as pH (figure)
<<HClTit4, fig=FALSE, echo=TRUE, eval=FALSE>>=
plot(ae, xval = ae$delta_conc_titrant, what = what,
xlab = "[HCl] offset added [mol/kg-soln]",
size = c(14,10), mfrow = c(4,4))
plot(ae, xval=ae$delta_moles_titrant, xlab = "HCl added [mol]",
what = what, size = c(14,10), mfrow = c(4,4))
@
<<HClTit5, fig=TRUE, echo=TRUE, width=10, height=7>>=
plot(ae, xval = ae$pH, xlab = "free scale pH", what = what,
size = c(12,8), mfrow = c(4,4))
@
The function \code{plot.aquaenv} also offers the possibility of creating Bjerrum plots from objects obtained with \code{titration}.
The simplest way to do that is:
<<HClTit6, fig=TRUE, echo=TRUE, width=10, height=7>>=
plot(ae, bjerrum = TRUE)
@
Or just select a few elements and customize the plotting style
<<HClTit7, fig=TRUE, echo=TRUE, width=10, height=7>>=
what <- c("CO2", "HCO3", "CO3")
plot(ae, what = what, bjerrum = TRUE, lwd = 4,
palette = c("cyan", "magenta", "yellow"),
bg = "gray", legendinset = 0.1, legendposition = "topleft")
@
Generally Bjerrum plots are done on the log scale. This can be accomplished using the flag \code{log}
<<HClTit9, fig=TRUE, echo=TRUE, width=10, height=7>>=
what <- c("CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH",
"NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F")
plot(ae, what = what, bjerrum = TRUE, log = TRUE)
@
Furthermore, we can zoom in to the region of most interest to marine scientists
<<HClTit10, fig=TRUE, echo=TRUE, width=10, height=7>>=
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)]))
@
\subsubsection{Titration with NaOH}
Similar to the titration with $\rm HCl$, also a titration with $\rm NaOH$ can be simulated
(Please note that due to runtime constraints of the vignette for this package, this titration simulation is not run during construction of the vignette
and, therefore, no resulting plots of the following code chunks is included. The user can generate the plots by extracting and executing the relevant code-chunks.)
<<NaOHTit1, fig=FALSE, echo=TRUE, eval=FALSE>>=
ae <- titration(aquaenv(S = 35, t = 15, SumCO2 = 0.0035, SumNH4 = 0.00002, pH = 2),
mass_sample = 0.01, mass_titrant = 0.02, conc_titrant = 0.01,
S_titrant = 0.5, steps = 50, type = "NaOH")
@
Plotting everything
<<NaOHTit2, fig=FALSE, echo=TRUE, eval=FALSE>>=
plot(ae, xval = ae$delta_mass_titrant, xlab = "NaOH solution added [kg]",
mfrow = c(10,10))
@
Plotting selectively
<<NaOHTit3, fig=FALSE, echo=TRUE, eval=FALSE>>=
what <- c("TA", "pH", "CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH",
"NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F", "fCO2")
plot(ae, xval = ae$delta_mass_titrant, xlab = "NaOH solution added [kg]",
what = what, size = c(12,8), mfrow = c(4,4))
plot(ae, xval = ae$pH, xlab = "free scale pH", what = what,
size = c(12,8), mfrow = c(4,4))
@
Bjerrum plots
<<NaOHTit4, echo=TRUE, eval=FALSE>>=
what <- c("CO2", "HCO3", "CO3")
plot(ae, what=what, bjerrum=TRUE)
@
<<NaOHTit5, echo=TRUE, eval=FALSE>>=
what <- c("CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH",
"NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F")
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)]))
@
\subsubsection{Titration with a titrant with high concentrations and a large sample volume - classical Bjerrum plots}
The Bjerrum plots created in the previous two sections do not really look like the classical textbook ones.
This is because we simulated a titration with a small sample volume and a titrant with low concentrations that is most representative to e.g. pore water titrations.
As a result the total concentrations like, e.g., total carbonate decreased due to dilution.
In simulating a titration with a rather large volume and a titrant with high concentrations the volume and salinity corrections do not matter
any more and graphs known from textbooks \citep[e.g.][]{Zeebe2001} are produced.
<<TextbookTit1, fig=FALSE, echo=TRUE>>=
ae <- titration(aquaenv(S = 35, t = 15, SumCO2 = 0.0035, SumNH4 = 0.00002, pH = 11.3),
mass_sample = 100, mass_titrant = 0.5, conc_titrant = 3,
S_titrant = 0.5, steps = 100)
@
Plotting everything (not shown)
<<TextbookTit2, fig=FALSE, echo=TRUE, eval=FALSE>>=
plot(ae, xval = ae$delta_mass_titrant,
xlab = "HCl solution added [kg]", mfrow = c(10, 10))
@
Plotting selectively and with different elements for \code{xval} (not shown)
<<TextbookTit3, fig=FALSE, echo=TRUE, eval=FALSE>>=
what <- c("TA", "pH", "CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH",
"NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F", "fCO2")
plot(ae, xval = ae$delta_mass_titrant, xlab = "HCl solution added [kg]",
what = what, size = c(12, 8), mfrow = c(4,4))
plot(ae, xval = ae$pH, xlab = "free scale pH", what = what,
size = c(12,8), mfrow = c(4,4))
plot(ae, xval = ae$delta_conc_titrant, xlab = "[HCl] offset added [mol/kg-soln]",
what = what, size = c(12,8), mfrow = c(4,4))
plot(ae, xval = ae$delta_moles_titrant, xlab = "HCl added [mol]",
what = what, size = c(12,8), mfrow = c(4,4))
@
Creating different kinds of Bjerrum plots (not shown)
<<TextbookTit4, fig=FALSE, echo=TRUE, eval=FALSE>>=
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)
@
and the classical textbook one
<<TextbookTit5, fig=TRUE, echo=TRUE, width=10, height=7>>=
what <- c("CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH",
"NH4", "NH3", "H2SO4", "HSO4", "SO4", "HF", "F")
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)]))
@
\subsection[Calculating information from titration curves: the function TAfit]{Calculating information from titration curves: the function \code{TAfit}}
\subsubsection{A little theory}
While titrating a sample of natural seawater with $\rm HCl$ one sees two clear equivalence points \citep{Dickson1981}.
The second equivalence point is the equivalence point of total alkalinity and the difference beween the second and the first
equivalence point signifies the total amount of $\sum \rm CO_2$ of the sample \cite{Hansson1973}.\\
\noindent
This can be illustrated with \aq$\,$. The respective titration curve can be plotted, together with its first and second derivative.
Furthermore, the equivalence points can be marked with vertical lines (Please note that for a titrant concentration of 0.01 mol/kg-solution
and 0.01 kg of sample, the value of the concentration (in mol/kg-solution) of total alkalinity and total carbonate equals the value of the total
amount (in mol)).
<<TAfit1, fig=TRUE, echo=TRUE, width=10, height=7>>=
ae_init <- aquaenv(S = 35, t = 15, SumCO2 = 0.0035, SumNH4 = 0.00002, pH = 11.3)
ae <- titration(ae_init, mass_sample = 0.01, mass_titrant = 0.02,
conc_titrant = 0.01, S_titrant = 0.5, steps = 100)
plot(ae, xval = ae$delta_mass_titrant, xlab = "HCl solution added [kg]",
what = "pH", xlim = c(0,0.015))
par(new = TRUE)
plot(ae$delta_mass_titrant[1:100], diff(ae$pH), type = "l",
col = "red", xlim = c(0,0.015), ylab = "", xlab = "", yaxt = "n")
par(new = TRUE)
plot(ae$delta_mass_titrant[2:100], diff(diff(ae$pH)), type = "l",
col = "blue", xlim = c(0,0.015), ylab = "", xlab = "", yaxt = "n")
abline(h =0, col = "blue")
abline(v = ae$TA[[1]])
abline(v = ae$TA[[1]] - ae$SumCO2[[1]])
@
Following classical chemical textbooks \citep[e.g.][]{Skoog1982}, one can determine [TA] and $\rm [\sum CO_2]$ of a sample
by graphically determining those equivalence points. However, there is no mechanistic understanding of the contents of the solution involved in doing so.\\
\noindent
Other methods, called ``Gran evaluations'' \citep{Gran1952, Hansson1973, Dickson1981, Haraldsson1997, Anderson1999},
try to linearize the mechanistic model of what is going on in the solution during titration.
They define the so called linear ``Gran functions'' and try to find their roots to determine the equivalence points.
We will illustrate that by plotting the Gran functions F0 (blue) and F2 (-F1, green) and again mark the equivalence points with vertical lines.
The y=zero line for the Gran functions is indicated by a horizontal line.
<<TAfit2, fig=TRUE, echo=TRUE, width=10, height=7>>=
plot(ae, xval = ae$delta_mass_titrant, xlab = "HCl solution added [kg]",
what = "pH", xlim = c(0,0.015))
prot1 <- c()
for (i in 1:length(ae$pH))
prot1 <- c(prot1, (10^-(ae$pH[[i]])+ae$HSO4[[i]]+ae$HF[[i]]+
ae$CO2[[i]]-ae$CO3[[i]]-ae$BOH4[[i]]-ae$OH[[i]]))
par(new = TRUE)
plot(ae$delta_mass_titrant, prot1, type = "l", col = "blue", xlim = c(0,0.015),
ylab = "", xlab = "", yaxt = "n", ylim = c(-0.015,0.015))
prot2 <- c()
for (i in 1:length(ae$pH))
prot2 <- c(prot2, (10^-(ae$pH[[i]])+ae$HSO4[[i]]+ae$HF[[i]]-
ae$HCO3[[i]]-2*ae$CO3[[i]]-ae$BOH4[[i]]-ae$OH[[i]]))
par(new = TRUE)
plot(ae$delta_mass_titrant, prot2, type = "l", col = "green", xlim = c(0,0.015),
ylab = "", xlab = "", yaxt = "n", ylim = c(-0.015,0.015))
abline(v = ae$TA[[1]])
abline(v = ae$TA[[1]] - ae$SumCO2[[1]])
abline(h = 0)
@
One can see that the Gran functions are not linear. This is due to volume and salinity change effects during the titration.
This can be overcome by either employing ``modified Gran functions'' \citep[see][]{Haraldsson1997} that correct for the volume changes or by using
a titration with a titrant with high concentrations and a large sample volume (Please note that here the value of the concentration of total
alkalinity and total carbonate does not equal their total amount and need to be converted with the factor 100/3)
<<TAfit3, fig=TRUE, echo=TRUE, width=10, height=7>>=
ae <- titration(aquaenv(S = 35, t = 15, SumCO2 = 0.0035, SumNH4 = 0.00002,
pH = 11.3),
mass_sample = 100, mass_titrant = 0.5, conc_titrant = 3,
S_titrant = 0.5, steps = 100)
plot(ae, xval = ae$delta_mass_titrant, xlab = "HCl solution added [kg]",
what = "pH", xlim = c(0, 0.5))
prot1 <- c()
for (i in 1:length(ae$pH))
prot1 <- c(prot1, (10^-(ae$pH[[i]])+ae$HSO4[[i]]+ae$HF[[i]]+
ae$CO2[[i]]-ae$CO3[[i]]-ae$BOH4[[i]]-ae$OH[[i]]))
par(new = TRUE)
plot(ae$delta_mass_titrant, prot1, type = "l", col = "blue",
xlim = c(0,0.5), ylab = "", xlab = "", yaxt = "n",
ylim = c(-0.015, 0.015))
prot2 <- c()
for (i in 1:length(ae$pH))
prot2 <- c(prot2, (10^-(ae$pH[[i]])+ae$HSO4[[i]]+ae$HF[[i]]-
ae$HCO3[[i]]-2*ae$CO3[[i]]-ae$BOH4[[i]]-ae$OH[[i]]))
par(new = TRUE)
plot(ae$delta_mass_titrant, prot2, type = "l", col = "green",
xlim = c(0,0.5), ylab = "", xlab = "", yaxt = "n",
ylim = c(-0.015,0.015))
abline(v = (ae$TA[[1]]*100/3))
abline(v = ((ae$TA[[1]] - ae$SumCO2[[1]])*100/3))
abline(h = 0)
@
Another proposed method of determining [TA] and $\rm [\sum CO_2]$ is to not only determine the two equivalence points, but to
fit the whole titration curve with a theoretical titration curve based on a mechanistic model
of what is going on in the solution during the titration \citep{Dickson1981, DOE1994, Anderson1999}.
The function \code{titration} of \aq$\,$ provides exactly such a theoretical titration curve and the function \code{TAfit} makes use of this fact to
determine [TA] and $\rm [\sum CO_2]$ of a sample by non linear curve fitting.
\subsubsection[Determining TA and SumCO2 by non linear curve fitting]{Determining TA and $\rm [\sum CO_2]$ by non linear curve fitting}
\paragraph{Proof of concept}$\,$\\
First, a proof of concept will show that the function \code{TAfit} is implemented consistently.
Some "data" can be generated with the \code{titration} function.
<<TAfit4, fig=FALSE, echo=TRUE>>=
initial_ae <- aquaenv(S = 35, t = 15, SumCO2 = 0.002, TA = 0.0022)
ae <- titration(initial_ae, mass_sample = 0.01, mass_titrant = 0.003,
conc_titrant = 0.01, S_titrant = 0.5, steps = 20)
@
Now, the input data for the \code{TAfit} routine can be generated: a table with the added mass of the titrant and the resulting free scale pH:
<<TAfit5, fig=FALSE, echo=TRUE>>=
titcurve <- cbind(ae$delta_mass_titrant, ae$pH)
@
Note that For the \code{TAfit} all total quantities except $\rm \sum CO_2$ ($\rm \sum NH_4^+$, $\rm \sum H_2S$, $\rm \sum H_3PO4$, $\rm \sum Si(OH)_4$, $\rm \sum HNO3$, $\rm \sum HNO2$, $\rm \sum B(OH)_3$, $\rm \sum H_2SO_4$, $\rm \sum HF$)
need to be known. However, the latter three can be calculated from salinity as it is done in this example.
<<TAfit6, fig=FALSE, echo=TRUE>>=
fit1 <- TAfit(initial_ae, titcurve, conc_titrant = 0.01,
mass_sample = 0.01, S_titrant = 0.5)
fit1
@
Thus, we see that \code{TAfit} calculates the correct $\rm \sum CO_2$ and $\rm TA$ values.\\
\noindent
Trying the \cite{Roy1993b} (K\_CO2 and K\_HCO3) and \cite{Perez1987} (K\_HF) values.
(Please note that due to runtime constraints of the vignette for this package, this calculation is not run during construction of the vignette.
The user can perform the calculation by extracting and executing the relevant code-chunks.)
<<TAfit6a, fig=FALSE, echo=TRUE, keep.source=TRUE, eval=FALSE>>=
initial_ae_ <- aquaenv(S = 35, t = 15, SumCO2 = 0.002, TA = 0.0022,
k1k2 = "lueker", khf = "perez")
ae_ <- titration(initial_ae_, mass_sample = 0.01, mass_titrant = 0.003,
conc_titrant = 0.01,
S_titrant = 0.5, steps = 20, k1k2 = "roy",
khf = "perez")
titcurve_ <- cbind(ae_$delta_mass_titrant, ae_$pH)
fit1_ <- TAfit(initial_ae_, titcurve_, conc_titrant = 0.01,
mass_sample = 0.01, S_titrant = 0.5, k1k2 = "roy",
khf = "perez", verbose = TRUE)
fit1_
@
\noindent
\code{TAfit} can also take E (V) values as input variables, so we generate E values using E0 = 0.4 V and the Nernst equation.
(But before that, we calculate a titration with a titrant with the same salinity as seawater such that S does not change during the titration.
Otherwise we would need to calculate the S profile for the titration extra to use it to convert to the total scale in the following step.)
However, to do so we first need to convert our pH curve to the seawater pH scale. According to \citet[p.7, ch.4, sop.3][]{DOE1994} and \cite{Dickson2007}, the Nernst equation relates E to
the total proton concentration.
<<TAfit7, fig=FALSE, echo=TRUE, eval=FALSE>>=
ae <- titration(initial_ae, mass_sample = 0.01, mass_titrant = 0.003,
conc_titrant = 0.01, steps = 20, seawater_titrant = TRUE)
titcurve <- cbind(ae$delta_mass_titrant, ae$pH)
tottitcurve <- cbind(ae$pH, convert(ae$pH, "pHscale", "free2tot", S = 35, t = 15))
Etitcurve <- cbind(ae$delta_mass_titrant,
(0.4 - ((PhysChemConst$R/10)*initial_ae$T/PhysChemConst$F)*
log(10^-tottitcurve[,2])))
@
Again, \code{TAfit} can be executed, this time also calculating E0. Note that the flag \code{verbose = TRUE} causes
\code{TAfit} to show the progress of the fitting procedure in a plot window.
(Please note that due to runtime constraints of the vignette for this package, this calculation is not run during construction of the vignette.
The user can perform the calculation by extracting and executing the relevant code-chunks.)
<<TAfit8, fig=FALSE, echo=TRUE, eval=FALSE>>=
fit2 <- TAfit(initial_ae, Etitcurve, conc_titrant = 0.01, mass_sample = 0.01,
Evals = TRUE, verbose = TRUE, seawater_titrant = TRUE)
fit2
@
Furthermore, \code{TAfit} can fit K\_CO2 as well, however, one single value for the whole titration curve is fitted, i.e. there is no correction for K\_CO2 changes due to changing S due to mixing with the titrant
<<TAfit9, fig=FALSE, echo=TRUE>>=
fit3 <- TAfit(initial_ae, titcurve, conc_titrant = 0.01, mass_sample = 0.01,
S_titrant = 0.5, K_CO2fit = TRUE)
fit3
initial_ae$K_CO2
@
One can see that the fitted value for K\_CO2 is not the same as the value in the initial \code{aquaenv} object, which is the "correct" value. That is, because during data creation K\_CO2 changed along the course of the titration due to changes in salinity. Assuming that the titrant has the same salinity as the sample (and is made up of natural seawater, i.e. containing SumBOH4, SumH2SO4 and SumHF as functions of S),
then the "correct" K\_CO2 should be fitted. This can be accomplished in \code{TAfit} by not giving the argument \code{S\_titrant} (i.e. assuming the titrant has the same salinity as the sample) and
setting the flag \code{seawater\_titrant} to \code{TRUE}
<<TAfit10, fig=FALSE, echo=TRUE>>=
ae <- titration(initial_ae, mass_sample = 0.01, mass_titrant = 0.003,
conc_titrant = 0.01, steps = 20, seawater_titrant = TRUE)
titcurve <- cbind(ae$delta_mass_titrant, ae$pH)
fit4 <- TAfit(initial_ae, titcurve, conc_titrant = 0.01, mass_sample = 0.01,
K_CO2fit = TRUE, seawater_titrant = TRUE)
fit4
@
Furthermore, TA, SumCO2, K\_CO2 and E0 can be fitted at the same time.
(Please note that due to runtime constraints of the vignette for this package, this calculation is not run during construction of the vignette.
The user can perform the calculation by extracting and executing the relevant code-chunks.)
<<TAfit11, fig=FALSE, echo=TRUE, eval=FALSE>>=
Etitcurve <- cbind(titcurve[,1], (0.4 - ((PhysChemConst$R/10)*initial_ae$T/
PhysChemConst$F)*log(10^-titcurve[,2])))
fit5 <- TAfit(initial_ae, Etitcurve, conc_titrant = 0.01, mass_sample = 0.01,
K_CO2fit = TRUE, seawater_titrant = TRUE, Evals = TRUE)
fit5
@
Sometimes, the obtained titration curve is not equally spaced on the x axis. \code{TAfit} can deal with such curves if the flag
\code{equalspaced} is set to \code{FALSE}. (Please note that due to runtime constraints of the vignette for this package, this calculation is not run during construction of the vignette.
The user can perform the calculation by extracting and executing the relevant code-chunks.)
<<TAfit12, fig=FALSE, echo=TRUE, eval=FALSE>>=
neqsptitcurve <- rbind(titcurve[1:9,], titcurve[11:20,])
fit6 <- TAfit(initial_ae, neqsptitcurve, conc_titrant = 0.01,
mass_sample = 0.01, seawater_titrant = TRUE, equalspaced = FALSE,
verbose = TRUE, debug = TRUE)
fit6
@
Finally, some "noise" is added to the generated data
<<TAfit13, fig=FALSE, echo=TRUE>>=
noisetitcurve <- titcurve * rnorm(length(titcurve), mean = 1, sd = 0.01) #one percent error possible
fit7 <- TAfit(initial_ae, noisetitcurve, conc_titrant = 0.01, mass_sample = 0.01,
seawater_titrant = TRUE, verbose = FALSE, debug = TRUE)
fit7
@
The flag \code{verbose=TRUE} (default is \code{FALSE}) prompts to show the traject of the fitting procedure in a plot window. However, each new fit is plotted over the first one
and Sweave includes only the first plot in each code chunk in the resulting \LaTeX file. Therefore, we use the flag \code{debug=TRUE} to
visualize the final fit
<<TAfit14, fig=TRUE, echo=TRUE, width=10, height=7>>=
ylim <- range(noisetitcurve[,2], calc)
xlim <- range(tit$delta_mass_titrant, noisetitcurve[,1])
plot(noisetitcurve[,1], noisetitcurve[,2], xlim = xlim, ylim = ylim,
type = "l", xlab = "delta mass titrant", ylab = "pH (free scale)")
par(new = TRUE)
plot(tit$delta_mass_titrant, calc, xlim = xlim, ylim = ylim, type = "l",
col = "red", xlab = "", ylab = "")
@
\paragraph{Test with generated data from literature}$\,$\\
\cite{Dickson1981} provided a synthetic dataset to test total alkalinity fitting programs. This dataset is included in \aq$\,$ as
\code{sample\_dickson1981}.
Following quantities are given
<<TAfit15, fig=FALSE, echo=TRUE>>=
conc_titrant <- 0.3 # mol/kg-soln
mass_sample <- 0.2 # kg
S_titrant <- 14.835 # is aequivalent to the ionic strength of 0.3 mol/kg-soln
SumBOH3 <- 0.00042 # mol/kg-soln
SumH2SO4 <- 0.02824 # mol/kg-soln
SumHF <- 0.00007 # mol/kg-soln
@
Note that all concentrations are in mol/kg-solution and the mass of the sample is in kg. Note further that the salinity of the titrant has been
calculated from its ionic strenght of 0.3 mol/kg-soln.\\
\noindent
In the original dataset as represented in \code{sample\_dickson1981}, the mass of titrant is given in g which needs to be converted to kg
<<TAfit16, fig=FALSE, echo=TRUE>>=
sam <- cbind(sample_dickson1981[,1]/1000, sample_dickson1981[,2]) # convert mass of titrant from g to kg
@
Then an attempt to recalculate the [TA] and $\rm [\sum CO_2]$ values given in \cite{Dickson1981} ([TA]=0.00245 mol/kg-soln and $\rm [\sum CO_2]$ 0.00220 mol/kg-soln)
can be done
<<TAfit17, fig=FALSE, echo=TRUE>>=
dicksonfit <- TAfit(aquaenv(S = 35, t = 25, SumBOH3 = SumBOH3,
SumH2SO4 = SumH2SO4, SumHF=SumHF), sam, conc_titrant,
mass_sample, S_titrant = S_titrant, debug = TRUE)
dicksonfit
@
This shows the fit is not accurate. Why is that so?
\subparagraph[Does the salinity correction (Stitrant) matter?]{Does the salinity correction (\code{S\_titrant}) matter?}$\;$\\
\noindent
Let us calculate a theoretical titration without salinity correction
<<TAfit18, fig=FALSE, echo=TRUE>>=
dicksontitration1 <- titration(aquaenv(S = 35, t = 25, SumCO2 = 0.00220,
SumBOH3 = SumBOH3, SumH2SO4 = SumH2SO4, SumHF = SumHF, TA = 0.00245),
mass_sample = mass_sample, mass_titrant = 0.0025,
conc_titrant = conc_titrant, steps = 50, type = "HCl")
@
and one with salinity correction
<<TAfit19, fig=FALSE, echo=TRUE>>=
dicksontitration2 <- titration(aquaenv(S = 35, t = 25, SumCO2 = 0.00220,
SumBOH3 = SumBOH3, SumH2SO4 = SumH2SO4, SumHF = SumHF, TA = 0.00245),
mass_sample = mass_sample, mass_titrant = 0.0025,
conc_titrant = conc_titrant, S_titrant = S_titrant, steps = 50, type = "HCl")
@
Now the difference between both curves (in red and blue) and the ``Dickson'' curve (in black) can be visualized
<<TAfit20, fig=TRUE, echo=TRUE, width=10, height=7>>=
plot(dicksontitration1, xval = dicksontitration1$delta_mass_titrant,
what = "pH", xlim = c(0,0.0025), ylim = c(3,8.2), col = "red",
xlab = "delta mass titrant")
par(new = TRUE)
plot(dicksontitration2, xval = dicksontitration2$delta_mass_titrant,
what = "pH", xlim = c(0,0.0025), ylim = c(3,8.2), col = "blue", xlab = "")
par(new = TRUE)
plot(sam[,1], sam[,2], type = "l", xlim = c(0,0.0025),
ylim = c(3,8.2), xlab = "", ylab = "")
@
That means, the salinity correction makes no significant difference (the red and the blue curve cannot be discerned), because the relation between the total amount of sample and the added amount of titrant is very large:
salinity only drops from 35 to 34.75105.\\
\noindent
But there is an offset between the "Dickson" curve and our curve
<<TAfit21, fig=TRUE, echo=TRUE, width=10, height=7>>=
plot(dicksontitration2$pH - sam[,2])
@
\subparagraph[Does fitting KCO2 as well improve the fit?]{Does fitting \code{K\_CO2} as well improve the fit?}$\;$\\
<<TAfit22, fig=FALSE, echo=TRUE>>=
dicksonfit2 <- TAfit(aquaenv(S = 35, t = 25, SumBOH3 = SumBOH3,
SumH2SO4 = SumH2SO4, SumHF = SumHF), sam, conc_titrant,
mass_sample, S_titrant = S_titrant, debug = TRUE, K_CO2fit = TRUE)
dicksonfit2
@
Yes it does, but it is not optimal yet.\\
<<>>=
@
\noindent
There still remains one major difference between the calculations carried out in \cite{Dickson1981} and the calculations in \aq$\,$:
\cite{Dickson1981} uses fixed values for the equilibrium constants and does not calculate them as functions of temperature and salinity.
Furthermore, the values that are used in \cite{Dickson1981} are not exactly the same as are obtained in \aq$\,$ for the same salinity and
temperature.\\
Let us calculate a theoretical titration curve employing exactly the same equilibrium constant values as used in \cite{Dickson1981} and plot the result together with the ``Dickson'' curve
<<TAfit23, fig=TRUE, echo=TRUE, width=10, height=7>>=
dicksontitration3 <- titration(aquaenv(S = 35, t = 25, SumCO2 = 0.00220,
SumBOH3 = SumBOH3, SumH2SO4 = SumH2SO4, SumHF = SumHF, TA = 0.00245,
k_w = 4.32e-14, k_co2 = 1e-6, k_hco3 = 8.20e-10, k_boh3 = 1.78e-9,
k_hso4 = (1/1.23e1), k_hf = (1/4.08e2)),
mass_sample = mass_sample, mass_titrant = 0.0025, conc_titrant = conc_titrant,
steps = 50, type = "HCl", S_titrant = S_titrant,
k_w = 4.32e-14, k_co2 = 1e-6, k_hco3 = 8.20e-10, k_boh3 = 1.78e-9,
k_hso4 = (1/1.23e1), k_hf = (1/4.08e2))
plot(dicksontitration3, xval = dicksontitration3$delta_mass_titrant,
what = "pH", xlim = c(0,0.0025), ylim = c(3,8.2), col = "blue",
xlab = "delta mass titrant")
par(new = TRUE)
plot(sam[,1], sam[,2], type = "l", xlim = c(0,0.0025), ylim = c(3,8.2),
xlab = "", ylab = "")
@
Plotting the differences between both curves reveals that they are the same down to 1 umol/kg-soln.
<<TAfit24, fig=TRUE, echo=TRUE, width=10, height=7>>=
plot(dicksontitration3$pH - sam[,2])
@
Calculating [TA] and $\rm [\sum CO_2]$ using \code{TAfit} and exactly the same equilibrium constant values as used in \cite{Dickson1981}
<<TAfit25, fig=FALSE, echo=TRUE>>=
dicksonfit3 <- TAfit(aquaenv(S = 35, t = 25, SumBOH3 = SumBOH3, SumH2SO4 = SumH2SO4,
SumHF = SumHF, k_w = 4.32e-14, k_co2 = 1e-6, k_hco3 = 8.20e-10,
k_boh3 = 1.78e-9, k_hso4 = (1/1.23e1), k_hf = (1/4.08e2)),
sam, conc_titrant, mass_sample, S_titrant = S_titrant, debug = TRUE,
k_w = 4.32e-14, k_co2 = 1e-6, k_hco3 = 8.20e-10, k_boh3 = 1.78e-9,
k_hso4 = (1/1.23e1), k_hf = (1/4.08e2))
dicksonfit3
@
reveals that now exactly the same values are calculated as are given in \cite{Dickson1981}.
\section[Extending AquaEnv]{Extending \textbf{\textsf{AquaEnv}}}
It is simple for the user to create own functions that use \aq$\,$ and extend its functionality.
We will demonstrate that by creating simple analogons for the \aq$\,$ functions \code{titration} and \code{TAfit}.\\
\noindent
The function \code{simpletitration} will take the following arguments\\
\begin{tabular}{lp{.7\textwidth}}
\code{aquaenv} & an object of class \code{aquaenv}: minimal definition, contains all information about the system: S, t, p, total concentrations of nutrients etc \\
\code{volume} & the volume of the (theoretical) titration vessel in l \\
\code{amount} & the amount of titrant added in mol\\
\code{steps} & the amount of steps the amount of titrant is added in \\
\code{type} & the type of titrant: either "HCl" or "NaOH"\\
\end{tabular}
The function is defined as
\begin{scriptsize}
<<extend1, fig=FALSE, echo=TRUE, keep.source=TRUE>>=
simpletitration <- function(aquaenv, # an object of class aquaenv: minimal definition,
# contains all information about the system:
# T, S, d, total concentrations of nutrients etc
volume, # the volume of the (theoretical) titration vessel in l
amount, # the amount of titrant added in mol
steps, # the amount of steps the amount of titrant is added in
type) # the type of titrant: either "HCl" or "NaOH"
{
directionTAchange <- switch(type, HCl = -1, NaOH = +1)
TAconcchangeperstep <- convert(((amount/steps)/volume), "conc", "molar2molin", aquaenv$t, aquaenv$S)
aquaenvtemp <- aquaenv
for (i in 1:steps)
{
TA <- aquaenvtemp$TA + (directionTAchange * TAconcchangeperstep)
aquaenvtemp <- aquaenv(ae=aquaenvtemp, TA=TA)
aquaenv <- merge(aquaenv, aquaenvtemp)
}
aquaenv[["DeltaCTitrant"]] <- convert((amount/volume)/steps*(1:(steps+1)),
"conc", "molar2molin", aquaenv$t, aquaenv$S)
return(aquaenv) # object of class aquaenv which contains a titration simulation
}
@
\end{scriptsize}
and can be used to create a bjerrum plot
<<extend2, fig=TRUE, echo=TRUE, width=10, height=7>>=
ae <- simpletitration(aquaenv(S = 35, t = 15, SumCO2 = 0.003500,
SumNH4 = 0.000020, pH = 11.3),
volume =100, amount = 1.5, steps = 100, type = "HCl")
what <- c("CO2", "HCO3", "CO3", "BOH3", "BOH4", "OH", "NH4", "NH3",
"H2SO4", "HSO4", "SO4", "HF", "F")
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)]))
@
The function \code{simpletitration} in turn can be used to create a simple analogon to \code{TAfit} with the arguments\\
\begin{tabular}{lp{.7\textwidth}}
\code{ae} & an object of class \code{aquaenv}: minimal definition, contains all information about the system: S, t, p, total concentrations of nutrients etc \\
\code{pHmeasurements} & a table containing the titration curve: basically a series of pH values (pH on free proton scale)\\
\code{volume} & the volume of the titration vessel\\
\code{amount} & the total amount of the titrant added\\
\code{TAguess}=0.0025 & a first guess for [TA] and [SumCO2] to be used as initial values for the optimization procedure\\
\code{type}="HCl"& the type of titrant: either "HCl" or "NaOH"\\
\end{tabular}\\
\noindent
defined as
\begin{scriptsize}
<<extend3, fig=FALSE, echo=TRUE, keep.source=TRUE>>=
simpleTAfit <- function(ae, # an object of class aquaenv: minimal definition,
# contains all information about the system:
# T, S, d, total concentrations of nutrients etc
pHmeasurements, # a table containing the titration curve:
# basically a series of pH values (pH on free proton scale)
volume, # the volume of the titration vessel
amount, # the total amount of the titrant added
TAguess=0.0025, # a first guess for [TA] and [SumCO2] to be used as
# initial values for the optimization procedure
type="HCl") # the type of titrant: either "HCl" or "NaOH"
{
ae$Na <- NULL # make sure ae gets cloned as "skeleton": cloneaquaenv determines "skeleton"
# TRUE or FALSE from the presence of a value for Na
residuals <- function(state)
{
ae$SumCO2 <- state[[1]]
pHcalc <- simpletitration(aquaenv(ae=ae, TA=state[[2]]), volume=volume,
amount=amount, steps=(length(pHmeasurements)-1), type=type)$pH
residuals <- pHmeasurements-pHcalc
return(residuals)
}
require(minpack.lm)
out <- nls.lm(fn=residuals, par=c(TAguess, TAguess)) #guess for TA is also used as guess for SumCO2
result <- list(out$par[[2]], out$par[[1]], out$deviance)
attr(result[[1]], "unit") <- "mol/kg-soln"
attr(result[[2]], "unit") <- "mol/kg-soln"
names(result) <- c("TA", "SumCO2", "sumofsquares")
return(result) # a list of three values
# ([TA] in mol/kg-solution, [SumCO2] in mol/kg-solution, sum of the squared residuals)
}
@
\end{scriptsize}
The function \code{simpleTAfit} can be used to calculate TA and SumCO2
<<extend4, fig=FALSE, echo=TRUE>>=
pHmeasurements <- ae$pH
fit <- simpleTAfit(aquaenv(S = 35, t = 15, SumNH4 = 0.00002),
pHmeasurements, volume = 100, amount = 1.5)
fit
@
<<CleaningUp, fig=FALSE, echo=FALSE>>=
graphics.off()
@
\appendix
\section{Abbreviations for references used throughout the code and in the helpfiles}
\begin{tabular}{ll}
Atkins1996 & \cite{Atkins1996} \\
Boudreau1996 & \cite{Boudreau1996}\\
DOE1994 & \cite{DOE1994}\\
Dickson1979a & \cite{Dickson1979a}\\
Dickson1981 & \cite{Dickson1981}\\
Dickson1984 & \cite{Dickson1984}\\
Dickson1987 & \cite{Dickson1987}\\
Dickson1990 & \cite{Dickson1990}\\
Dickson2007 & \cite{Dickson2007}\\
Emerson2008 & \cite{Emerson2008}\\
Feistel2008 & \cite{Feistel2008}\\
Fofonoff1983 & \cite{Fofonoff1983}\\
Follows2006 & \cite{Follows2006}\\
Hofmann2008 & \cite{Hofmann2008}\\
Khoo1977 & \cite{Khoo1977}\\
Lewis1998 & \cite{Lewis1998}\\
Lueker2000 & \cite{Lueker2000}\\
Millero1981 & \cite{Millero1981}\\
Millero1988 & \cite{Millero1988}\\
Millero1995 & \cite{Millero1995}\\
Millero1995a & \cite{Millero1995a}\\
Millero2006 & \cite{Millero2006}\\
Mucci1983 & \cite{Mucci1983}\\
Perez1987a & \cite{Perez1987a}\\
Riordan2005 & \cite{Riordan2005}\\
Roy1993b & \cite{Roy1993b}\\
Sundquist1979 & \cite{Sundquist1979}\\
Weiss1970 & \cite{Weiss1970}\\
Weiss1974 & \cite{Weiss1974}\\
Wischmeyer2003 & \cite{Wischmeyer2003}\\
Zeebe2001 & \cite{Zeebe2001}\\
\end{tabular}
\clearpage
\section[References for the elements of an object of class aquaenv]{References for the elements of an object of class \code{aquaenv}} \label{app: references}
\begin{footnotesize}
\begin{longtable}{p{.15\textwidth}|p{.9\textwidth}}
\textbf{element}& \textbf{references} \\ \hline
\code{p}, \code{P}, \code{Pa}, \code{p} & The relation between pressure and depth given in \cite{Fofonoff1983} is used. The standard value for
atmospheric pressure \code{Pa} at sea level as well as the definition of total pressure and gauge pressure is taken from \cite{Feistel2008}.\\
\code{Cl} & \citet[chapter 5, p. 11]{DOE1994}, and \citet[p. 100, footnote 3]{Zeebe2001}\\
\code{I} & \citet[chapter 5, p. 13, 15]{DOE1994}, \citet[p.12]{Zeebe2001}, and \citet[ p.257]{Roy1993b}.
Note that the approximation I/(mol/kg-solution) $\approx$ 0.0199201 \; \code{S} is given in \citet[p. 428.]{Millero1982}.
This relationship converted into mol/kg-$\rm H_2O$ and the last digits adjusted (from 0.0199201 to 0.019924) results in the formula used here.\\
density & \citet{Millero1981} and \citet[chapter 5, p. 6f]{DOE1994}.\\
\code{Br}, \code{ClConc}, \code{Na}, \code{Mg}, \code{Ca}, \code{K}, \code{Sr} & \citet[ chapter 5, p.11]{DOE1994}\\
\code{molal2molin} & \citet[ p.257]{Roy1993b}, and \citet[chapter 5, p. 15]{DOE1994}\\
\code{free2tot}, \code{tot2free} & \citet[p.2302]{Dickson1984}, \citet[chapter 4, p.16]{DOE1994}, \citet[p.57, 261]{Zeebe2001}\\
\code{free2sws}, \code{tot2sws}, \code{sws2free}, \code{sws2tot} & \citet[p.2303]{Dickson1984}, \citet[p.57]{Zeebe2001}\\
\code{K0\_CO2} & \citet{Weiss1974}, \citet[chapter 5, p. 13]{DOE1994} (here it is stated that the unit is mol/(kg-solution*atm)), \citet[ p.663]{Millero1995}, \citet[ p.257]{Zeebe2001}\\
\code{K0\_O2} & derived from a formula for the oxygen saturation concentration in ml-$\rm O_2$/kg-solution by \citet{Weiss1970} using the first virial coefficient
of oxygen \citep[][p. 41, 1029]{Atkins1996} and the atmospheric oxygen fugacity \citep{Williams2004}\\
\code{K\_W} & \citet[p.670]{Millero1995} (\textbf{original reference}, but slightly different formula for seawater pH), \citet[chapter 5, p. 18]{DOE1994} (NOT the original reference! \citet{DOE1994} cites in an update from 1997 \citet{Millero1995}!
However the version of the formula used here is the one converted to total pH scale given in \citet{DOE1994}), and \citet[p. 258]{Zeebe2001}. Constant type (stoichiometric),
pH scale (total, converted to free here) , and
concentration unit (mol/kg-solution squared): \citet[chapter 5, p. 12,18]{DOE1994}, pH scale also in \citet[p. 258]{Zeebe2001}.\\
\code{K\_HSO4} & \citet[chapter 5 page 13]{DOE1994}, \citet[p. 260]{Zeebe2001}, \citet{Dickson1990a} (original reference). Constant type (stoichiometric), pH scale (free) , and
concentration unit (mol/kg-$\rm H_2O$ converted to mol/kg-solution here): \citet[chapter 5, p. 13]{DOE1994}. Note that it is also possible to use the constant according to \cite{Khoo1977}, as cited in, e.g., \cite{Roy1993b}, \cite{Millero1995}, and \cite{Lewis1998}. In \cite{Lewis1998} it is stated that the constant resulting from this equation is in mol/kg-H$_2$O and on the free pH scale.\\
\code{K\_HF} & \citet[p. 91]{Dickson1979} (original reference), \citet[c. 5, p. 15]{DOE1994}, \citet[p. 257]{Roy1993b}, \citet[p. 1783]{Dickson1987}, \citet[p. 664]{Millero1995}, \citet[p. 260]{Zeebe2001}
(converted to molinty and total scale). Constant type (stoichiometric), pH scale (free) , and
concentration unit (mol/kg-$\rm H_2O$ converted to mol/kg-solution here): \citet[chapter 5, p. 15, 16]{DOE1994}. In \textsf{AquaEnv}, it is also possible to use the constant according to \cite{Perez1987a}. \\
\code{K\_CO2}, \code{K\_HCO3} & \citet[p. 254]{Roy1993b} (original reference), \citet[chapter 5, p.14]{DOE1994} (in a version converted to mol/kg-$\rm H_2O$), \citet[p. 664]{Millero1995}, \citet[p. 255]{Zeebe2001}.
Constant type (stoichiometric) and concentration unit (mol/kg-$\rm H_2O$ converted to mol/kg-solution here): \citet[chapter 5, p. 14, 15]{DOE1994}, pH scale (total, converted to free here):
In \citet[chapter 5, p. 12]{DOE1994} the total scale is stated for the formula for high salinities and thus can be inferred for the formula for low salinities.
The scale is also indirectly stated for both formulations in the original reference \citet{Roy1993b}. Note that in \citet{Roy1993b} a function for fresh water (based on \cite{Millero1979} which in turn is on a temperature relationship from \cite{Harned1943} and \cite{Harned1941} respectively) and a function for seawater is derived. In \citet{Millero1995} it is stated that for S$<$5 the fresh water formula (based on \cite{Millero1979}) should be used and for S$>=$5 the seawater formula derived in \citet{Roy1993b}. However, both formulations do not always
intersect at \code{S}=5. The true intersection with respect to salinity \code{S} is a function of temperature \code{t}. Here, we first calculate this intersection by
numerical root finding and then decide which formulation to use. This practise results in a continuous function with respect to \code{S}. (Note that there is a typesetting error in \citet{Roy1993b}: one of the numerical values for the function for $\rm K^*_{CO_2}$ is given as 310.48919, but correct is 2310.48919. However, in \citet{Millero1995} this value is stated correctly.) In \textsf{AquaEnv}, it is also possible to use the constants according to \cite{Lueker2000} and \cite{Millero2006}.\\
\code{K\_BOH3} & \citet[p. 763]{Dickson1990} (original, but mol/kg-$\rm H_2O$ version), \citet[ch. 5, p. 14]{DOE1994}, \citet[p. 262]{Zeebe2001}, \citet[p.669]{Millero1995} (mol/kg-$\rm H_2O$ version) , agrees with data in \citet{Roy1993a}.
Constant type (stoichiometric) and concentration unit (mol/kg-solution): \citet[chapter 5, p. 14]{DOE1994}, pH scale (total): \citet[chapter 5, p. 12]{DOE1994} and \citet[p.263]{Zeebe2001}. \\
\code{K\_NH4} & \citet{Millero1995a} (original reference), \citet[p.671]{Millero1995}. Constant type (stoichiometric) and concentration unit (mol/kg-solution): \citet[p.671]{Millero1995}, pH scale (seawater, converted to free here):
\citet{Lewis1998} (in corrections of \citet{Millero1995}).\\
\code{K\_H2S} & \citet{Millero1988} (original reference), \citet[p.671]{Millero1995}. Constant type (stoichiometric) and concentration unit (mol/kg-solution): \citet[p.671]{Millero1995}, pH scale (seawater, converted to free here):
\citet{Lewis1998} (in corrections of \citet{Millero1995}).\\
\code{K\_H3PO4}, \code{K\_H2PO4}, \code{K\_HPO4} & \citet[p.670]{Millero1995} (original reference, but formula for seawater scale pH), \citet[ch. 5, p 16,17]{DOE1994}, agrees with data in \citet{Dickson1979a}.
Constant type (stoichiometric), concentration unit (mol/kg-solution), and pH scale (total, converted to free here): \citet[chapter 5, p. 12, 16, 17]{DOE1994}.\\
\code{K\_SiOH4} & \citet{Millero1988} (original reference), \citet[chapter 5, p 17]{DOE1994}, \citet[p.671]{Millero1995} (formula for seawater scale pH)
Constant type (stoichiometric), concentration unit (mol/kg-$\rm H_2O$ converted to mol/kg-solution here by omitting the conversion summand ln(1-0.001005 S)), and pH scale (total, converted to free here): \citet[chapter 5, p. 12, 17]{DOE1994}.\\
\code{K\_SiOOH3} & \citet{Wischmeyer2003} (original reference), corrected due to personal communication with Dieter Wolf-Gladrow (one of the authors). The corrected version can be obtained from either Dieter Wolf-Gladrow or Andreas F Hofmann
(a.hofmann@nioo.knaw.nl). Constant type (stoichiometric), concentration unit (mol/kg-solution), and pH scale (total, converted to free here): \citet{Wischmeyer2003}.\\
\code{K\_HNO2} & Constant value, not a function of temperature and salinity! Obtained as a hybrid pk value (featuring the activity of the proton but the concentration of other species (see \citet{Zeebe2001} for a treatment of different types of
equilibrium constants) in molar concentration (mol/l) on the NBS pH scale \citep{Durst1975} from \citet{Riordan2005}. Used as an approximation for the stoichiometric $\rm K^*_{\rm HNO_2}$ in mol/kg-solution on the free proton pH scale here.\\
\code{K\_H2SO4} & Constant value, not a function of temperature and salinity! Obtained as a standard pK value from \citet[p. 1045]{Atkins1996}. Used as an approximation for the stoichiometric $\rm K^*_{\rm H_2SO_4}$ in mol/kg-solution on the free proton pH scale here.\\
\code{K\_HS} & Constant value, not a function of temperature and salinity! Obtained as a standard pK value from \citet[p. 1045]{Atkins1996}. Used as an approximation for the stoichiometric $\rm K^*_{\rm Hs^-}$ in mol/kg-solution on the free proton pH scale here.\\
\code{Ksp\_calcite}, \code{Ksp\_aragonite} & \citet{Mucci1983} (original reference), \citet{Boudreau1996}. Note that in there are errors in \citet{Boudreau1996}: $b_0$ for calcite is not 0.7712 but 0.77712 and $b_1$ for
aragonite is not 0.001727 but 0.0017276.\\
\code{pH} & As given in \citet{Dickson1984}, p. 2303 (use of "m") and \citet{Dickson1979a}, p. 91f all concentrations appearing in the definition
of the total and the seawater pH scale are \textbf{molal} (mol/kg-$\rm H_2O$) concentrations. But in \citet{Roy1993b}, p. 257 and in \citet{DOE1994},
chapter 4, SOP 6, p. 1 it is stated, that concentrations for the seawater and total pH scale are in mol/kg-solution.
To be consistent with \citet{DOE1994} \textbf{molin} concentrations (mol/kg-solution) are chosen for calculating the pH.\\
\code{revelle} & N/A (function redundant in current version of \aq)\\
\code{dTAdH}, \code{dTAdSumCO2}, \code{dTAdSumBOH3}, \code{dTAdSumH2SO4}, \code{dTAdSumHF}, \code{dTAdSumH3PO4}, \code{dTAdSumSumSiOH4}, \code{dTAdSumH2S}, \code{dTAdSumNH4}, \code{dTAdSumHNO3}, \code{dTAdSumHNO2} & \citet{Hofmann2008}\\
\code{c1}, \code{c2}, \code{c3}, \code{b1}, \code{b2}, \code{so1}, \code{so2}, \code{so3}, \code{f1}, \code{f2}, \code{p1}, \code{p2}, \code{p3}, \code{p4}
\code{si1}, \code{si2}, \code{si3}, \code{s1}, \code{s2}, \code{s3}, \code{n1}, \code{n2}, \code{na1}, \code{na2}, \code{ni1}, \code{ni2}
& \citet{Skoog1982}, \citet{Stumm1996}, \citet{Hofmann2010a}\\
\code{dTAdKdKdS}, \code{dTAdKdKdT}, \code{dTAdKdKdp}, \code{dTAdKdKdSumH2SO4}, \code{dTAdKdKdSumHF}
& \citet{Hofmann2009}\\ \hline
\end{longtable}
\end{footnotesize}
\noindent
The values for \code{K\_W},\code{K\_HSO4}, \code{K\_HF}, \code{K\_CO2}, \code{K\_HCO3}, \code{K\_BOH3}, \code{K\_NH4}, \code{K\_H2S}, \code{K\_H3PO4}, \code{K\_H2PO4}, \code{K\_HPO4}, \code{K\_SiOH4}, \code{K\_SiOOH3}, \code{Ksp\_calcite}, \code{Ksp\_aragonite} obtained as functions of
salinity \code{S} and temperature \code{t} from the above references are pressure corrected using the gauge pressure \code{p} according to \citet{Millero1995}
with corrections by \citet{Lewis1998}.\\
\noindent
In general it is to be said that all corrections from \cite{Lewis1998} have been applied.
\bibliography{AquaEnv}
\end{document}