Calcite solubility in water, rainwater, and seawater#

Written by Svetlana Kyas (ETH Zurich) on Mar 30th, 2022

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!

Studying the solubility of carbon dioxide in rainwater is important for understanding climate change. Anthropogenic CO2 absorbed by the oceans increases the concentration of hydrogen ions (H+), and thus its acidification. In a review article by Figuerola et al., A Review and Meta-Analysis of Potential Impacts of Ocean Acidification on Marine Calcifiers From the Southern Ocean, Front. Mar. Sci., 2021, ocean acidification (OA) (accompanied by CO reduction, see figure below) is identified as a critical problem for the shells and skeletons of marine calcifiers such as foraminifera, corals, echinoderms, mollusks, and bryozoans. The OA process reduces the carbonate saturation horizon (the depth below which calcium carbonate dissolves), which likely increases the vulnerability of many resident marine calcifiers to dissolution. In addition, ocean warming could further exacerbate the effects of the OA process on these particular species. For this reason, understanding the dependence of calcite solubility on various factors is an important issue in modeling.

In this tutorial, we examine the solubility of calcite in water, rainwater, and seawater. We also look at how it is affected by changes in temperature, pressure and the amount of carbon dioxide dissolved.

Infographic of the ocean acidification process

Infographic of the ocean acidification process, Source: frontiersin.org

First, we initialize the chemical system with aqueous, gaseous, and calcite phases.

from reaktoro import *

# Create the database
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(ActivityModelPitzer())

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

# Create a mineral phase
mineral = MineralPhase("Calcite")

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

Next, the equilibrium specifications, equilibrium conditions, and equilibrium solver for the equilibrium calculations are initialized. To constrain the charge of the chemical state, we need to make it open to the Cl-. Finally, we create aqueous properties to evaluate pH in the forthcoming calculations.

# Define equilibrium specs
specs = EquilibriumSpecs (system)
specs.temperature()
specs.pressure()
specs.charge()
specs.openTo("Cl-")

# Define conditions to be satisfied at the chemical equilibrium state
conditions = EquilibriumConditions(specs)
conditions.charge(0.0, "mol") # to make sure the mixture is charge neutral

# Define the equilibrium solver
solver = EquilibriumSolver(specs)

# Define aqueous properties
aprops = AqueousProps(system)

The functions below define the chemical states corresponding to pure water, rainwater saturated with carbon dioxide, and seawater, respectively:

def water():

    state = ChemicalState(system)
    state.add("H2O(aq)", 1.0, "kg")

    return state

def rainwater():

    state = ChemicalState(system)
    # Rainwater composition
    state.set("H2O(aq)", 1.0, "kg")
    state.set("Na+"    , 2.05, "mg") # Sodium, 2.05 ppm = 2.05 mg/L ~ 2.05 mg/kgw
    state.set("K+"     , 0.35, "mg") # Potassium
    state.set("Ca+2"   , 1.42, "mg") # Calcium
    state.set("Mg+2"   , 0.39, "mg") # Magnesium
    state.set("Cl-"    , 3.47, "mg") # Chloride
    state.set("SO4-2"  , 2.19, "mg")
    state.set("NO3-"   , 0.27, "mg")
    state.set("NH4+"   , 0.41, "mg")
    state.set("CO2(aq)", 0.36, "mol")  # rainwater is saturated with CO2

    return state

def seawater():

    state = ChemicalState(system)
    # Seawater composition
    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")

    return state

Solubility of calcite for different temperatures and pressures#

To calculate the solubility of calcite in a given brine, we define the function solubility_of_calcite(), which calculates the pH of the aqueous solution along with the amount of calcite dissolved in it. As the names suggest, the instances of pandas.DataFrame, i.e., df_water and df_rainwater, correspond to the results of chemical equilibrium calculations for water and rainwater, respectively.

import pandas as pd
df_water = pd.DataFrame(columns=["T", "P", "pH", "deltaCalcite"])
df_rainwater = pd.DataFrame(columns=["T", "P", "pH", "deltaCalcite"])

# Initial amount of calcite
n0Calcite = 10.0

def solubility_of_calcite(state, T, P, tag):

    conditions.temperature(T, "celsius")
    conditions.pressure(P, "bar")

    # Equilibrate the solution with given initial chemical state and desired conditions at equilibrium
    res = solver.solve(state, conditions)

    # Check calculation succeeded
    assert res.succeeded()

    # Update aqueous properties
    aprops.update(state)

    # Fetch the amount of final calcite in the equilibrium state
    nCalcite = float(state.speciesAmount("Calcite"))

    if tag == "water":
        df_water.loc[len(df_water)] = [T, P, float(aprops.pH()), n0Calcite - nCalcite]
    elif tag == "rainwater":
        df_rainwater.loc[len(df_rainwater)] = [T, P, float(aprops.pH()), n0Calcite - nCalcite]

Next, we initialize the array of temperatures from 20 °C till 90 °C and pressures P = 1, 10, 100 bars and perform solubility calculations for water and rainwater.

import numpy as np
temperatures = np.arange(20.0, 91.0, 5.0)   # in celsius
pressures = np.array([1, 10, 100])          # in bar

state = water()
state.set("Calcite", n0Calcite, "mol")
[solubility_of_calcite(state, T, P, "water") for P in pressures for T in temperatures];

state = rainwater()
state.set("Calcite", n0Calcite, "mol")
[solubility_of_calcite(state, T, P, "rainwater") for P in pressures for T in temperatures];

Let us check that the calcite solubility determined for 25 °C and 1 bar actually corresponds to the values from Wikipedia, i.e., the solubility in water is 0.013 g/L (25 °C). In the following, we use 100.0869 g/mol as the molar mass of calcite:

# Fetch Calcite solubility for T = 25 celsius and P = 1 bar
deltaCalcite = df_water[(df_water['P']== 1.0) & (df_water['T'] == 25.0)]['deltaCalcite'].iloc[0]
print(f"Solubility of calcite in water equals to {1e3*deltaCalcite:.3f} "
      f"mmol/kgw = ... = {deltaCalcite * 0.1000869 * 1e3:.3f} g/L")

Solubility of calcite in water equals to 0.116 mmol/kgw = ... = 0.012 g/L

In the case of pure water, therefore, about 0.116 mmol of calcite will dissolve. This dissolved quantity is independent of the initial quantity value used for calcite (provided it exceeds the solubility limit).

We plot below the solubility of calcite in water and in CO2-saturated rainwater at 1 bar:

from reaktplot import *

df_water_P1 = df_water[df_water['P'] == 1.0]
df_rainwater_P1 = df_rainwater[df_rainwater['P'] == 1.0]

fig = Figure()
fig.title("CALCITE SOLUBILITY IN WATER AND<br>CO2-SATURATED RAINWATER AT 1 BAR")
fig.xaxisTitle('Temperature [°C]')
fig.yaxisTitle('Calcite Solubility [mol/kgw]')
fig.drawLineWithMarkers(df_water_P1["T"], df_water_P1["deltaCalcite"], name="water")
fig.drawLineWithMarkers(df_rainwater_P1["T"], df_rainwater_P1["deltaCalcite"], name="rainwater")

fig.show()