Carbon dioxide (CO₂) dissolution in seawater#

Written by Svetlana Kyas (ETH Zurich) on Mar 31th, 2022.
Last revised on March 31st, 2023 by Allan Leal.

Attention

Always make sure you are using the latest version of Reaktoro. Otherwise, some new features documented on this website will not work on your machine and you may receive unintuitive errors. Follow these update instructions to get the latest version of Reaktoro!

Ocean acidification remains a significant problem resulting from human activities, primarily the burning of fossil fuels. As atmospheric carbon dioxide levels increase, the ocean absorbs more carbon dioxide, triggering a series of chemical reactions in seawater that have a negative impact on the ocean and its inhabitants. These reactions can lead to a decrease in the production of calcium carbonate shells in bivalves and other aquatic organisms, as well as a range of other physiological challenges for marine organisms. This reduction in shell production has implications for the entire ocean food chain, as many species rely on these shells for protection and as a food source. Therefore, addressing ocean acidification is critical for the health and sustainability of ocean ecosystems.

This tutorial demonstrates the use of Reaktoro for modeling ocean acidification, by computing the pH of a seawater solution for increasing concentrations of dissolved CO₂.

The CO2 cycle between the atmosphere and the ocean

The CO2 cycle between the atmosphere and the ocean, Source: wikipedia.org

We start by defining the chemical system containing aqueous and gaseous phases.

from reaktoro import *
import pandas as pd

db = SupcrtDatabase("supcrtbl")

# Create an aqueous phase automatically selecting all species with provided elements
aqueousphase = AqueousPhase(speciate("H O C Ca Mg K Cl Na S N"))
aqueousphase.set(chain(
    ActivityModelHKF(),
    ActivityModelDrummond("CO2"),
))

# Create a gaseous phase
gaseousphase = GaseousPhase("CO2(g)")
gaseousphase.set(ActivityModelPengRobinson())

# Create the chemical system
system = ChemicalSystem(db, aqueousphase, gaseousphase)

# Create the equilibrium solver
solver = EquilibriumSolver(system)

Next, we initialize the chemical state corresponding to the seawater content, equilibrate it, and evaluate the pH level obtained after equilibration:

state = ChemicalState(system)
state.setTemperature(25, "celsius")
state.setPressure(1.0, "bar")
state.add("H2O(aq)",      1.0, "kg")
state.add("Ca+2" ,   412.3, "mg")
state.add("Mg+2" ,  1290.0, "mg")
state.add("Na+"  , 10768.0, "mg")
state.add("K+"   ,   399.1, "mg")
state.add("Cl-"  , 19353.0, "mg")
state.add("HCO3-",   141.7, "mg")
state.add("SO4-2",  2712.0, "mg")

solver.solve(state)

aprops = AqueousProps(state)
print("pH of seawater = ", float(aprops.pH()))

pH of seawater =  7.699061790440688

Finally, we define the auxiliary variables with the initial values of the CO2 amount and its increment. We run the loop with nsteps steps adding CO2 into seawater and re-equilibrate it again. We collect the amount of added carbon dioxide and the corresponding pH value in pandas.DataFrame.

co2_0 = 0.0
co2_delta = 0.1
nsteps = 50

df = pd.DataFrame(columns=["amountCO2", "pH"])
df.loc[len(df)] = [co2_0, float(aprops.pH())]

for i in range(nsteps):

    # Add more CO2 to the problem
    state.add("CO2(g)", co2_delta, "mmol")

    # Equilibrate state with updated problem
    solver.solve(state)

    # Update aqueous properties
    aprops.update(state)

    # Update CO2 amount
    co2_0 += co2_delta

    # Append new calculated value to the dataframe
    df.loc[len(df)] = [co2_0, float(aprops.pH())]

***WARNING***
The chemical equilibrium/kinetics calculation did not converge.
This can occur due to various factors, such as:
  1. Infeasible Problem Formulation: The problem you have defined lacks a feasible solution within the given constraints and assumptions.
  2. Out-of-Range Thermodynamic Conditions: The conditions you have specified for the system may fall outside the valid range supported by the thermodynamic models used.
  3. Numerical Instabilities: Convergence issues may arise from numerical problems during the execution of the chemical/kinetic equilibrium algorithm. Consider reporting the issue with a minimal reproducible example if you believe the algorithm is responsible for this issue.
Disable this warning message with Warnings.disable(906) in Python and Warnings::disable(906) in C++.
***WARNING***
The chemical equilibrium/kinetics calculation did not converge.
This can occur due to various factors, such as:
  1. Infeasible Problem Formulation: The problem you have defined lacks a feasible solution within the given constraints and assumptions.
  2. Out-of-Range Thermodynamic Conditions: The conditions you have specified for the system may fall outside the valid range supported by the thermodynamic models used.
  3. Numerical Instabilities: Convergence issues may arise from numerical problems during the execution of the chemical/kinetic equilibrium algorithm. Consider reporting the issue with a minimal reproducible example if you believe the algorithm is responsible for this issue.
Disable this warning message with Warnings.disable(906) in Python and Warnings::disable(906) in C++.
***WARNING***
The chemical equilibrium/kinetics calculation did not converge.
This can occur due to various factors, such as:
  1. Infeasible Problem Formulation: The problem you have defined lacks a feasible solution within the given constraints and assumptions.
  2. Out-of-Range Thermodynamic Conditions: The conditions you have specified for the system may fall outside the valid range supported by the thermodynamic models used.
  3. Numerical Instabilities: Convergence issues may arise from numerical problems during the execution of the chemical/kinetic equilibrium algorithm. Consider reporting the issue with a minimal reproducible example if you believe the algorithm is responsible for this issue.
Disable this warning message with Warnings.disable(906) in Python and Warnings::disable(906) in C++.
***WARNING***
The chemical equilibrium/kinetics calculation did not converge.
This can occur due to various factors, such as:
  1. Infeasible Problem Formulation: The problem you have defined lacks a feasible solution within the given constraints and assumptions.
  2. Out-of-Range Thermodynamic Conditions: The conditions you have specified for the system may fall outside the valid range supported by the thermodynamic models used.
  3. Numerical Instabilities: Convergence issues may arise from numerical problems during the execution of the chemical/kinetic equilibrium algorithm. Consider reporting the issue with a minimal reproducible example if you believe the algorithm is responsible for this issue.
Disable this warning message with Warnings.disable(906) in Python and Warnings::disable(906) in C++.

We plot below pH as a function of the added amount of CO2 into the seawater.

from reaktplot import *

fig = Figure()
fig.title("PH DEPENDENCE ON AMOUNT OF ADDED CO2 TO THE SEAWATER")
fig.xaxisTitle('CO2 [mmol]')
fig.yaxisTitle('pH')
fig.drawLineWithMarkers(df["amountCO2"], df["pH"], name="")

fig.show()

As expected, the pH decreases as we continue to add CO2 into the seawater, the solution gradually becomes more acidic.