#!/usr/bin/env python
# coding: utf-8

# # Phosphate accumulation in carbonate-rich brines
# 
# <p class="acknowledgement">Written by Svetlana Kyas (ETH Zurich) on Mar 31th, 2022</p>
# 
# ```{attention}
# Always make sure you are using the [latest version of Reaktoro](https://anaconda.org/conda-forge/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](updating_reaktoro_via_conda) 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 CO<sub>2</sub> 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**](geobiology-phreeqc-fixed-fugacity.ipynb)). 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):

# In[ ]:


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:

# In[ ]:


# 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<sup>-13</sup>.

# In[ ]:


# 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 log<sub>10</sub>(pCO<sub>2</sub>) = −3.5 to 0 bar. We note that up to 2.1 moles phosphate occurs in equilibrium with Na<sub>2</sub>(HPO<sub>4</sub>)·7H<sub>2</sub>O salts.
# 
# The following block defines the array of CO<sub>2</sub> 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:

# In[ ]:


# 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"))]


# The modeled pH of saturated phosphate brines depends on the temperature and the partial CO<sub>2</sub> pressures. At present-day pCO<sub>2</sub> levels (log<sub>10</sub>(pCO<sub>2</sub>) = −3.5), solutions are highly alkaline (pH approximate to 10), consistent with high pHs measured in modern soda lakes. However, in CO<sub>2</sub>-rich atmospheres on the early Earth (log<sub>10</sub>(pCO<sub>2</sub>) = −2 to 0), brines range from moderately alkaline (with pH = 9) to slightly acidic (pH = 6.5) because of acidification by CO<sub>2</sub> (see the plot below). Below, we plot pH levels with the [bokeh](https://bokeh.org/) 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 CO<sub>2</sub> is more soluble in solutions at lower temperatures, making pH slightly higher.

# In[ ]:


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 CO<sub>2</sub>  of 1 bar and becomes alkaline (pH ~ 9) at a partial pressure of CO<sub>2</sub> of 0.01 bar (ppCO2 = 2). The relatively high pCO<sub>2</sub> values acidify the solution, which suggests that increased phosphate concentrations could have occurred in CO<sub>2</sub>-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 CO<sub>3</sub><sup>−2</sup> vs. HCO<sub>3</sub><sup>−</sup> ions in saturated Na–HCO<sub>3</sub>–CO<sub>3</sub> brines and the pH of the solution is controlled by the pCO<sub>2</sub> according to our model. Below, we plot pH and pCO<sub>2</sub> as a function of the equivalent percentage of CO<sub>3</sub><sup>−2</sup> 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}]}}$.

# In[ ]:


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 CO<sub>2</sub> partial pressure.