Phosphate accumulation in carbonate-rich brines
Phosphate accumulation in carbonate-rich brines#
Written by Svetlana Kyas (ETH Zurich) on Mar 31th, 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!
The simulations performed in this tutorial are based on the article “A carbonate-rich lake solution to the phosphate problem of the origin of life” by Toner and Catling (2020) and propose a possible approach to the phosphate problem in the origin of life. The phosphate problem arises due to the low phosphate solubility as it tends to precipitate in a form of apatite in the presence of calcium. At the same time, it is one of the most important compounds in the formation of the first cellular compounds such as phospholipids or nucleic acids.
However, Toner2020 show that early Earth conditions differed in many ways from today, including significantly higher carbon dioxide levels in the atmosphere. Elevated CO2 partial pressures may have led to carbonation of brines and then precipitation of calcium carbonate in the brines on the early Earth. Lower calcium concentrations in the brines would then have allowed higher phosphate concentrations, creating suitable conditions for the formation of the first cells on Earth.
Note
This tutorial is one of two tutorials that follow the paper Toner2020 and attempts to replicate the geobiological simulations performed in it (see also the second tutorial Phosphate accumulation in carbonate-rich brines). This work was done in collaboration with Cara Magnabosco and Laura Murzakhmetov, ETH-Zurich.
The thermodynamic database is loaded from the database file provided by the publication of Toner and Catling (2020):
from reaktoro import *
import numpy as np
import pandas as pd
db = PhreeqcDatabase.fromFile('phreeqc-toner-catling.dat') # if running from tutorials folder
We define a chemical system based on the database and provided aqueous and mineral phases. Moreover, to evaluate pH and phosphate amount in the aqueous phase, we will need aqueous and chemical properties:
# Define aqueous phase
solution = AqueousPhase(speciate("H O C Na Cl P"))
solution.set(chain(
ActivityModelHKF(),
ActivityModelDrummond("CO2")
))
# Define mineral phases
minerals = MineralPhases("Natron Nahcolite Trona Na2CO3:H2O Na2CO3:7H2O Halite Na2(HPO4):7H2O")
# Define chemical system
system = ChemicalSystem(db, solution, minerals)
# Define aqueous and chemical properties
props = ChemicalProps(system)
aprops = AqueousProps(system)
To tell the solver that fugacity in this chemical system will be constrained, we have to define equilibrium specifications and corresponding conditions. The first specifies what will be a constraint and the second by what value (specified below for the range of fugacities). We also reset the equilibrium option’s field epsilon to reset the default lower bound of species to 10-13.
# Define equilibrium specifications
specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
specs.fugacity("CO2")
# Define conditions to be satisfied at the chemical equilibrium state
conditions = EquilibriumConditions(specs)
# Define equilbirium conditions
opts = EquilibriumOptions()
opts.epsilon = 1e-13
To determine the maximum possible phosphate concentrations in such brine, we model solutions saturated with sodium phosphate, carbonate, and chloride salts at temperatures between 0 and 50 °C and gas pressure of log10(pCO2) = −3.5 to 0 bar. We note that up to 2.1 moles phosphate occurs in equilibrium with Na2(HPO4)·7H2O salts.
The following block defines the array of CO2 partial pressures and the data blocks that will store the results for different temperatures. We perform equilibrium calculations for different pressures in the for loop:
# Auxiliary arrays
num_log10pCO2s = 71
co2pressures = np.flip(np.linspace(-5.0, 2.0, num=num_log10pCO2s))
temperatures = np.array([0, 25, 50])
# Output dataframe
data = pd.DataFrame(columns=["T", "ppCO2", "pH", "mCO3", "mHCO3", "x", "amount_P"])
for T in temperatures:
for log10pCO2 in co2pressures:
conditions.temperature(T, "celsius")
conditions.pressure(1.0, "atm")
conditions.fugacity("CO2", 10 ** log10pCO2, "bar")
state = ChemicalState(system)
state.set("H2O" , 1.0, "kg")
state.set("Nahcolite" , 10.0, "mol")
state.set("Halite" , 10.0, "mol")
state.set("Na2(HPO4):7H2O", 10.0, "mol")
state.set("CO2" , 100.0, "mol")
solver = EquilibriumSolver(specs)
solver.setOptions(opts)
# Equilibrate the solution with the given initial chemical state and desired conditions at the equilibrium
res = solver.solve(state, conditions)
# Stop if the equilibration did not converge or failed
if not res.succeeded(): continue
props.update(state)
aprops.update(state)
mCO3 = float(state.speciesAmount("CO3-2"))
mHCO3 = float(state.speciesAmount("HCO3-"))
x = 100 * 2 * mCO3 / (mHCO3 + 2 * mCO3)
data.loc[len(data)] = [T, log10pCO2, float(aprops.pH()),
mCO3, mHCO3, x,
float(props.elementAmountInPhase("P", "AqueousPhase"))]
***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++.
The modeled pH of saturated phosphate brines depends on the temperature and the partial CO2 pressures. At present-day pCO2 levels (log10(pCO2) = −3.5), solutions are highly alkaline (pH approximate to 10), consistent with high pHs measured in modern soda lakes. However, in CO2-rich atmospheres on the early Earth (log10(pCO2) = −2 to 0), brines range from moderately alkaline (with pH = 9) to slightly acidic (pH = 6.5) because of acidification by CO2 (see the plot below). Below, we plot pH levels with the bokeh plotting library. We note that those are the maximum values for the corresponding solution, as it is saturated with respect to carbonate minerals. For undersaturated solutions, the pH is always lower. We see that temperature also affects the pH levels, because CO2 is more soluble in solutions at lower temperatures, making pH slightly higher.
from bokeh.plotting import figure, show, gridplot
from bokeh.models import HoverTool
from bokeh.io import output_notebook
from bokeh.models import ColumnDataSource
output_notebook()
# ----------------------------------- #
# Plot P amount
# ----------------------------------- #
hovertool1 = HoverTool()
hovertool1.tooltips = [("pH", "@pH"),
("T", "@T °C"),
("ppCO2", "@ppCO2")]
p1 = figure(
title="DEPENDENCE PH ON TEMPERATURE",
x_axis_label=r'PARTIAL PRESSURE CO2 [-]',
y_axis_label='PH [-]',
sizing_mode="scale_width",
height=300)
p1.add_tools(hovertool1)
colors = ['teal', 'darkred', 'indigo']
for T, color in zip(temperatures, colors):
df_T = ColumnDataSource(data[data['T'] == T])
p1.line("ppCO2", "pH", legend_label=f'{T} °C', line_width=2, line_cap="round", line_color=color, source=df_T)
p1.legend.location = 'top_right'
# ----------------------------------- #
# Plot Ca amount
# ----------------------------------- #
hovertool2 = HoverTool()
hovertool2.tooltips = [("amount(P)", "@amount_P mol"),
("T", "@T °C"),
("ppCO2", "@ppCO2")]
p2 = figure(
title="DEPENDENCE PHOSPHATE AMOUNT ON TEMPERATURE",
x_axis_label=r'PARTIAL PRESSURE CO2 [-]',
y_axis_label='AMOUNT OF P IN SOLUTION [MOL]',
sizing_mode="scale_width",
height=300)
p2.add_tools(hovertool2)
colors = ['teal', 'darkred', 'indigo']
for T, color in zip(temperatures, colors):
df_T = ColumnDataSource(data[data['T'] == T])
p2.line("ppCO2", "amount_P", legend_label=f'{T} °C', line_width=2, line_cap="round", line_color=color, source=df_T)
p2.legend.location = 'top_left'
grid = gridplot([[p1], [p2]])
show(grid)
We see that the pH is neutral around a partial pressure of CO2 of 1 bar and becomes alkaline (pH ~ 9) at a partial pressure of CO2 of 0.01 bar (ppCO2 = 2). The relatively high pCO2 values acidify the solution, which suggests that increased phosphate concentrations could have occurred in CO2-rich atmospheres on the early Earth. We also plot the amount of phosphorus in the brine supporting this hypothesis. The second plot also indicates that the solubility of phosphates increases with growing temperature and pressure.
The relative proportion of CO3−2 vs. HCO3− ions in saturated Na–HCO3–CO3 brines and the pH of the solution is controlled by the pCO2 according to our model. Below, we plot pH and pCO2 as a function of the equivalent percentage of CO3−2 ions relative to the total carbonate alkalinity defined as \(\mathsf{x = \frac{2\, [\rm{CO}_3^{-2}]}{[\rm{HCO}_3^-] + 2\, [\rm{CO}_3^{-2}]}}\).
from bokeh.models import LinearAxis, Range1d, Legend
df_T = data[data['T'] == 25.0]
df_T_pH = df_T["pH"]
hovertool = HoverTool()
hovertool.tooltips = [("pH", "@pH"),
("ppCO2", "@ppCO2"),
("x", "@x")]
p = figure(y_range=(df_T_pH.iloc[0]-0.1, df_T_pH.iloc[-1]+0.1),
title="PH AND PCO2 DEPENDENCE ON X",
x_axis_label=r'x',
y_axis_label='PH [-]',
sizing_mode="scale_width",
height=300)
p.add_tools(hovertool)
r11 = p.line("x", "pH", line_width=3, line_cap="round", line_color="midnightblue", source=df_T)
r12 = p.circle("x", "pH", size=6, fill_color=None, line_color="midnightblue", source=df_T)
p.extra_y_ranges = {"foo": Range1d(start=co2pressures[-1]-0.1, end=1.01*co2pressures[0]+0.1)}
r21 = p.line("x", "ppCO2", y_range_name="foo", line_width=2, line_cap="round", line_color="deeppink", source=df_T)
r22 = p.square("x", "ppCO2", y_range_name="foo", size=6, fill_color=None, line_color="deeppink", source=df_T)
legend = Legend(items=[
("pH" , [r11, r12]),
("ppCO2", [r21, r22])
], location="center")
p.add_layout(legend, 'right')
p.add_layout(LinearAxis(y_range_name="foo",axis_label="ppCO2 [-]"), 'right')
show(p)
We see the increase of the brine alkalinity with respect to increate of x. At the same time, higher x corresponds to the lower CO2 partial pressure.