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 pH

Computing $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.8490818837613

This 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.1758936301019

We may also be interested in the pH so we can request that be outputted:

carbon_chemistry(; DIC, Alk, T, S, return_pH = true)

8.001440897031806

These 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.0230153445105

And 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.961060028668

The 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.4572891858634

Finally, 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 pH

But 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.812912879999

The 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 calcite_saturation

using OceanBioME.Models.CarbonChemistryModel: calcite_saturation

calcite_saturation(carbon_chemistry; DIC, Alk, T, S)

0.001301242857483965

This 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²)S

K1() # 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 * T

KW() # 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 * S

KS() # 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) S

KP2() # 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) S

KP3() # 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) S

KSi() # 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_correction

Equilibrium 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 / RT

K2().pressure_correction

Equilibrium constant pressure correction
    ΔV = -15.82 + -0.0219 Tc + 0.0 * Tc²
    Δκ = 0.00113 + -0.0001475 T
    ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RT

KB().pressure_correction

Equilibrium constant pressure correction
    ΔV = -29.48 + 0.1622 Tc + -0.002608 * Tc²
    Δκ = -0.00284 + 0.0 T
    ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RT

KW().pressure_correction

Equilibrium 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 / RT

KS().pressure_correction

Equilibrium 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 / RT

KF().pressure_correction

Equilibrium 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 / RT

KP1().pressure_correction

Equilibrium 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 / RT

KP2().pressure_correction

Equilibrium 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 / RT

KP3().pressure_correction

Equilibrium 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 / RT

KSP_calcite().pressure_correction

Equilibrium constant pressure correction
    ΔV = -48.76 + 0.5304 Tc + -0.0 * Tc²
    Δκ = -0.01176 + 0.0003692 T
    ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RT

KSP_aragonite().pressure_correction

Equilibrium 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