Carbon Chemistry
Our carbon chemistry model follows the best practice described by Dickson et al. (2007), as implemented by e.g. cbsyst (Branson, 2023) and CO2SYS (Humphreys et al., 2022) in Python.
The carbon chemistry model is primarily used for diagnosing the partial pressure of carbon dioxide in the surface water for gas exchange with the air, but is capable of diagnosing other species such as carbonate concentration, useful in the calculation of calcite dissolution. The model works by computing the pH from the total dissolved inorganic carbon and total alkalinity (and optionally silicate, phosphate, boron, sulfate, and fluoride concentration), along with temperature and salinity, which can then be used to find the concentration of different species. We will first describe how to use the model, followed by the underlying chemistry and parameterisations.
Using the model
To use the carbon chemistry model it is constructed by simply writing:
using OceanBioME
carbon_chemistry = CarbonChemistry()`CarbonChemistry` model which solves for pCO₂ and pHComputing $fCO_2$ and $pH$
To compute the fugacity of carbon dioxide ($fCO_2$) you call the model with the DIC, Alkalinity, Temperature, and Salinity:
DIC = 2145
Alk = 2448
T = 25.4
S = 36.45
carbon_chemistry(; DIC, Alk, T, S)460.8490818837613This is sufficient when computing $fCO_2$ at the surface, but if we wanted to know $fCO_2$ at depth where there is higher pressure we can specify the pressure in bar like:
carbon_chemistry(; DIC, Alk, T, S, P = 100)457.1758936301019We may also be interested in the pH so we can request that be outputted:
carbon_chemistry(; DIC, Alk, T, S, return_pH = true)8.001440897031806These function calls assume a constant boron, sulfate, and fluoride ratio relative to the salinity (as described below), but can be specified instead:
carbon_chemistry(; DIC, Alk, T, S, boron = 0)305.0230153445105And the silicate and phosphate concentrations are assumed to be zero but can similarly be specified:
carbon_chemistry(; DIC, Alk, T, S, silicate = 2, phosphate = 1)462.961060028668The same code can also be used to compute $fCO_2$ when the pH is already known by passing it in the same way:
carbon_chemistry(; DIC, pH = 8.1, T, S)358.4572891858634Finally, for slightly improved accuracy you may wish to specify that the seawater density calculation, is an intermediary step in the calculations above, is computed using the full TEOS-10 Feistel (2008) standard by setting this as the density_function when you construct the carbon chemistry model:
using OceanBioME.Models: teos10_density
carbon_chemistry = CarbonChemistry(; density_function = teos10_density)`CarbonChemistry` model which solves for pCO₂ and pHBut this comes at the cost of significantly increased computational expense and GPU incompatability.
During the conversion from practical to absolute salinity units the location can then also be entered for (very slightly, in this example ~$10^{-4}\mu$atm) improved accuracy:
carbon_chemistry(; DIC, Alk, T, S, lon = -31.52, lat = 33.75)460.812912879999The default uses the polynomial approximation described in Roquet et al. (2015) as provided by SeawaterPolynomials.jl.
Computing the carbonate concentration
So that this model can be used in calcite dissolution models it can also return the carbonate saturation by calling the function carbonate_saturation
using OceanBioME.Models.CarbonChemistryModel: carbonate_saturation
carbonate_saturation(carbon_chemistry; DIC, Alk, T, S)8.188275727436483e-39This function takes all of the same arguments (e.g. boron) as carbon_chemistry above.
Chemistry
pH computation
When carbon dioxide is dissolved in seawater it dissociates into carbonate and bicarbonate species according to the equilibria:
\[\ce{CO_2(g)} \ce{<=> CO_2(aq)},\]
\[\ce{CO_2(aq)} + \ce{H_2O} \ce{<=> H_2CO_3(aq)},\]
\[\ce{CO_2(aq) + H_2O}\ce{<=> H^+ HCO^-_3},\]
\[\ce{HCO^-_3} \ce{<=> H^+ + CO^{2-}_3},\]
from which we define the total dissolved inorganic carbon (DIC) to be $DIC = [\ce{CO_2(aq)}] + [\ce{HCO^-_3}] + [\ce{CO^{2-}_3}]$. The equilibrium constants for these equations are defined as:
\[K_0=\frac{[\ce{CO_2(aq)}]}{\ce{pCO_2}},\]
\[K_1=\frac{\ce{[HCO^-_3][H^+]}}{\ce{[CO_2(aq)]}},\]
and
\[K_2=\frac{\ce{[CO^{2-}_3][H^+]}}{\ce{[HCO^-_3]}}.\]
These equilibria depend on the total hydrogen ion concentration $[H^+]$, which depends on the concentration of acids and buffering of bases.
Alkalinity: acid-base buffering
Alkalinity is a measure of the buffering capacity of seawater and we can approximate it as:
\[Alk = [\ce{HCO_3^-}] + 2[\ce{CO_3^{2-}}] + [\ce{B(OH)_4^-}] + [\ce{OH^-}] + [\ce{HPO_4^{2-}}] + 2[\ce{PO_4^{3-}}] + [\ce{SiO(OH)_3^-}] + [\ce{NH_3}] + [\ce{HS^-}] - [\ce{H^+}] - [\ce{HSO_4^-}] - [\ce{HF}] - [\ce{H_3PO_4}] + \text{ minor acids and bases},\]
where the "minor" species are either in sufficiently low concentrations as to be neglected (sometimes the photphate, silicate and sulfate species can also be neglected), and we shall neglect $[\ce{NH_3}]$.
Each of these acid and base species are also in an equilibrated dissociation state given below.
Hydrogen sulfate
Sulfate in the form of hydrogen sulfate dissociates to form sulfate ions in the equilibrium
\[\ce{HSO_4^-}\ce{<=> H^+ + SO_4^{2-}},\]
with an equilibrium constant given by
\[K_S=\frac{\ce{[SO_4^{2-}][H^+]}}{\ce{[HSO_4^-]}}.\]
Boric acid
Boron in the form of boric acid equilibrates with water in the reaction
\[\ce{B(OH)_3+H_2O}\ce{<=> H^+ + B(OH)_4^-},\]
with equilibrium constant
\[K_B=\frac{\ce{[B(OH)_4^-][H^+]}}{\ce{[B(OH)_3]}}.\]
Hydrogen fluoride
Hydrogen fluoride dissociated in the equilibrium
\[\ce{HF}\ce{<=> H^+ + F^-},\]
with equilibrium constant
\[K_F=\frac{\ce{[F^-][H^+]}}{\ce{[HF]}}.\]
Phosphoric acid
Phosphoric acid undergoes a three stage dissociation given by the equilibrium
\[\ce{H_3PO_4}\ce{<=> H^+ + H_2PO_4^-},\]
\[\ce{H_2PO_4^-}\ce{<=> H^+ + HPO_4^{2-}},\]
\[\ce{HPO_4^{2-}}\ce{<=> H^+ + PO_4^{3-}},\]
with equilibrium constants
\[K_{P1} = \frac{\ce{[H^+][H_2PO_4^-]}}{\ce{[H_3PO_4]}},\]
\[K_{P2} = \frac{\ce{[H^+][HPO_4^{2-}]}}{\ce{[H_2PO_4^-]}},\]
\[K_{P3} = \frac{\ce{[H^+][PO_4^{3-}]}}{\ce{[HPO_4^{2-}]}}.\]
Silicic acid
Silicic acid dissociates in the equilibrium
\[\ce{Si(OH)_4}\ce{<=> H^+ + SiO(OH)_3^-},\]
with equilibrium constant
\[K_{Si} = \frac{\ce{[H^+][SiO(OH)_3^-]}}{\ce{[Si(OH)_4]}}.\]
Water
Finally, water dissociates in the equilibrium
\[\ce{H_2O}\ce{<=> H^+ + OH^-},\]
with equilibrium constant
\[K_w = \ce{[H^+][OH^-]}.\]
Alkalinity equilibration
From these rate constants we can rewrite the total alkalinity (given above) in terms of only the rate constants, total hydrogen ion concentration ($\ce{H^+}$), the total dissolved inorganic carbon ($[DIC]$), boron ($[\text{B}]$), phosphate ($[\text{P}]$), silicate ($[\text{Si}]$), Sulfate ($[\text{Sulf}]$), and fluorine ($[\text{F}]$) content of the water, by rearranging the equations above. This results in a form of the total alkalinity:
\[\begin{align} Alk &\approx \frac{[DIC][\ce{H^+}]K_1}{[\ce{H^+}]^2 + K_1[\ce{H^+}] + K_1K_2}\\ &+ \frac{2[DIC]K_1K_2}{[\ce{H^+}]^2 + K_1[\ce{H^+}] + K_1K_2}\\ &+ \frac{[\text{B}]}{1+[\ce{H^+}]/K_B}\\ &+ \frac{K_w}{[\ce{H^+}]}\\ &+ \frac{[\text{P}][\ce{H^+}]K_{P1}K_{P2}}{[\ce{H^+}]^3+K_{P1}[\ce{H^+}]^2+K_{P1}K_{P2}[\ce{H^+}] + K_{P1}K_{P2}K_{P3}}\\ &+ \frac{2[\text{P}]K_{P1}K_{P2}K_{P3}}{[\ce{H^+}]^3+K_{P1}[\ce{H^+}]^2+K_{P1}K_{P2}[\ce{H^+}] + K_{P1}K_{P2}K_{P3}}\\ &+ \frac{[\text{Si}]}{1+[\ce{H^+}]/K_{Si}}\\ &- \frac{[\ce{H^+}]}{1+\text{S}/K_S}\\ &- \frac{[\text{Sulf}]}{1+K_S/[\ce{H^+}]/(1+S/K_S)}\\ &- \frac{[\text{F}]}{1+K_F/[\ce{H^+}]}\\ &- \frac{[\ce{H^+}]^2}{[\ce{H^+}]^3+K_{P1}[\ce{H^+}]^2+K_{P1}K_{P2}[\ce{H^+}] + K_{P1}K_{P2}K_{P3}}. \end{align}\]
This gives us a second equation from total alkalinity (which we must already know) with one unknown, $[\ce{H^+}]$. Our model solves for $[\ce{H^+}]$ using a bisection method to an accuracy of $10^{-10}$, having approximated the equilibrium constants from parametrisations described below. We can then determine the pH as, $pH = -\log_{10}([\ce{H^+}])$.
Carbon dioxide
Now that we know $[DIC]$, $Alk$, and $[\ce{H^+}]$ we can return to the equation for total dissolved inorganic carbon to find the concentration of aqueous carbon dioxide
\[[CO_2(aq)] = \frac{[DIC][\ce{H^+}]^2}{[\ce{H^+}]^2 + K_1[\ce{H^+}] + K_1K_2},\]
from which we determine gas phase carbon dioxide concentration as
\[fCO_2 = \frac{[CO_2(aq)]}{K_0},\]
in atmospheres.
Carbonate concentration and calcite saturation
Similarly we can also diagnose the calcite concentration
\[[CO_3^{2-}] = \frac{[DIC]K_1K_2}{[\ce{H^+}]([\ce{H^+}] + K_1)+K_1K_2}.\]
This concentration is important in the dissolution of calcium carbonate which reacts according to the equilibrium
\[\ce{CaCO_3(aq)}\ce{<=> Ca^+ + CO_3^{2-}},\]
which has an equilibrium constant
\[K_{SP} = [\ce{Ca^{2+}}][\ce{CO_3^{2-}}].\]
The calcite saturation can then be defined as $\Omega=\frac{[\ce{CO_3^{2-}}]}{[\ce{CO_{3, saturation}^{2-}}]}$ which can be diagnosed as:
\[\Omega = \frac{[\ce{Ca^{2+}}][\ce{CO_3^{2-}}]}{K_{SP}}.\]
Missing pieces
In most cases the chemistry described above requires information about more elements that is usually available. This means that we must parameterise their concentrations, usually this results in the boron, sulfate, and fluoride concentrations being set as constant ratios of the salinity. Usually these ratios are:
\[\begin{align} [\text{B}] &= \frac{0.000232}{10.811}\frac{S}{1.80655},\\ [\text{S}] &= \frac{0.14}{96.06} \frac{S}{1.80655},\\ [\text{F}] &= \frac{0.000067}{18.9984} \frac{S}{1.80655}.\\ \end{align}\]
We use these ratios by default in the model (but they can be changed when calling the models functions described below). Additionally, silicate and phosphate concentrations are often unavailable (both in observations and models), so are by default set to zero and their alkalinity contribution assumed to be small.
The "ionic strength" must also be parameterised for some equilibrium constants and is usually assumed to be:
\[Is = \frac{19.924}{1000 - 1.005S},\]
but these parameters are also variable in the model.
Model parameterisation
The chemical system described above has a large number of equilibrium constants, the constants are typically assumed to depend on the temperature, salinity, and pressure (hydrostatic pressure from depth underwater). [//]: # (This sentence is not currently actually true, but will be by the time of PR merge) By default, this model parameterises them based on the "best practice" guidelines of Dickson et al. (2007). These parameterisations are:
using OceanBioME.Models.CarbonChemistryModel: K0, K1, K2, KB, KW, KS, KF, KP1, KP2, KP3, KSi, KSP_calcite, KSP_aragonite
K0() # Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)Solubility constant
ln(k₀/k°) = -60.2409 + 9345.17 / T + 23.3585 (log(T) - log(100)) + 0.0 T² + (0.023517 + -0.00023656 T + 4.7036e-7 T²)SK1() # Millero (1995, Geochim. Cosmochim. Acta, 59, 664)First carbon dioxide dissociation constant
log₁₀(k₁/k°) = 61.2172 + -3633.86 / T + -9.6777 log(T) + 0.011555 S + -0.0001152 S²K2() # Millero (1995, Geochim. Cosmochim. Acta, 59, 664)Second carbon dioxide dissociation constant
log₁₀(k₂/k°) = -25.929 + -471.78 / T + 3.16967 log(T) + 0.01781 S + -0.0001122 S²KB() # Dickson (1990, Deep-Sea Res., 37, 755–766)Boric acid dissociation constant
ln(kᵇ/k°) = 148.0248 + (-8966.9 + -2890.53 √S + -77.942 S + 1.728 √S³ + -0.0996 S²) / T
+ 137.1942 * √S
+ 1.62142 * S
+ (-24.4344 + -25.085 √S + -0.2474 S ) * log(T)
+ 0.053105 * √S * TKW() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677)Water dissociation constant
ln(kʷ/k°) = 148.9652
+ -13847.26 / T
+ -23.6521 log(T)
+ (-5.977 + 118.67 / T + 1.0495 log(T)) √S
+ -0.01615 * SKS() # Dickson (1990, Chem. Thermodyn., 22, 113–127)Bisulfate dissociation constant
ln(kˢ/k°) = 141.328
+ -4276.1 / T
+ -23.093 log(T)
+ (324.57 + -13856.0 / T + -47.986 log(T)) √Is
+ (-771.54 + 35474.0 / T + 114.723 log(T)) Is
+ -2698.0 √Is³ / T
+ 1776.0 Is² / T
+ log(1 + -0.001005 S)KF() # Dickson and Riley (1979, Mar. Chem., 7, 89–99)Hydrogen fluoride dissociation constant
ln(kᶠ/k°) = -9.68
+ 874.0 / T
+ 0.111 √S
+ log(1 + 0.0 S)
+ log(1 + 0.0 S / Kˢ)KP1() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677)Phosphate dissociation constant
ln(kᵖⁿ/k°) = 115.525
+ -4576.752 / T
+ -18.453 log(T)
+ (0.69171 + -106.736 / T) √S
+ (-0.01844 + -0.65643 / T) SKP2() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677)Phosphate dissociation constant
ln(kᵖⁿ/k°) = 172.0883
+ -8814.715 / T
+ -27.927 log(T)
+ (1.3566 + -160.34 / T) √S
+ (-0.05778 + 0.37335 / T) SKP3() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677)Phosphate dissociation constant
ln(kᵖⁿ/k°) = -18.141
+ -3070.75 / T
+ 0.0 log(T)
+ (2.81197 + 17.27039 / T) √S
+ (-0.09984 + -44.99486 / T) SKSi() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677)Silicic acid constant
ln(kˢⁱ/k°) = 117.385
+ -8904.2 / T
+ -19.334 log(T)
+ (3.5913 + -458.79 / T) √Is
+ (-1.5998 + 188.74 / T) Is
+ (0.07871 + -12.1652 / T) Is²
+ log(1 + -0.001005 S)KSP_calcite() # Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341)Calcite solubility
ln(kₛₚ) = -171.9065 + -0.077993 T + 2839.319 / T + 71.595 log(T)
ln(kₛₚˢ) = ln(kₛₚ) + (-0.77712 + 0.0028426 T + 178.34 / T) √S
+ -0.07711 S + 0.0041249 √S³KSP_aragonite() # Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341)Calcite solubility
ln(kₛₚ) = -171.945 + -0.077993 T + 2903.293 / T + 71.595 log(T)
ln(kₛₚˢ) = ln(kₛₚ) + (-0.068393 + 0.0017276 T + 88.135 / T) √S
+ -0.10018 S + 0.0059415 √S³The pressure corrections from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341) are:
K1().pressure_correctionEquilibrium constant pressure correction
ΔV = -25.5 + 0.1271 Tc + 0.0 * Tc²
Δκ = -0.00308 + 8.77e-5 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTK2().pressure_correctionEquilibrium constant pressure correction
ΔV = -15.82 + -0.0219 Tc + 0.0 * Tc²
Δκ = 0.00113 + -0.0001475 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKB().pressure_correctionEquilibrium constant pressure correction
ΔV = -29.48 + 0.1622 Tc + -0.002608 * Tc²
Δκ = -0.00284 + 0.0 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKW().pressure_correctionEquilibrium constant pressure correction
ΔV = -20.02 + 0.1119 Tc + -0.001409 * Tc²
Δκ = -0.00513 + 7.94e-5 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKS().pressure_correctionEquilibrium constant pressure correction
ΔV = -18.03 + 0.0466 Tc + 0.000316 * Tc²
Δκ = -0.00453 + 9.0e-5 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKF().pressure_correctionEquilibrium constant pressure correction
ΔV = -9.78 + -0.009 Tc + -0.000942 * Tc²
Δκ = -0.00391 + 5.4e-5 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKP1().pressure_correctionEquilibrium constant pressure correction
ΔV = -14.51 + 0.1211 Tc + -0.000321 * Tc²
Δκ = -0.00267 + 4.27e-5 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKP2().pressure_correctionEquilibrium constant pressure correction
ΔV = -23.12 + 0.1758 Tc + -0.002647 * Tc²
Δκ = -0.00515 + 9.0e-5 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKP3().pressure_correctionEquilibrium constant pressure correction
ΔV = -26.57 + 0.202 Tc + -0.003042 * Tc²
Δκ = -0.00408 + 7.14e-5 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKSP_calcite().pressure_correctionEquilibrium constant pressure correction
ΔV = -48.76 + 0.5304 Tc + -0.0 * Tc²
Δκ = -0.01176 + 0.0003692 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RTKSP_aragonite().pressure_correctionEquilibrium constant pressure correction
ΔV = -45.96 + 0.5304 Tc + -0.0 * Tc²
Δκ = -0.01176 + 0.0003692 T
ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RT