kinetic rate function yields different results than phreeqc when ph changes with time

[uliw]

when I write a rate function where the addition of an element also changes the pH of the solution, the result will differ from what I get with phreeqc. A good example is the oxidation of pyrite which depends on the activity of H+ and Fe+3, and will also change both during the dissolution process.

Looking at the code, it appears that the activities will always be taken from the initial solution, and as such do not reflect the fact that the oxidation rate changes in concert with both activities.

Is there an easy way to calculate such a system with phreeqpython? I can provide code examples if need be.

[AbelHeinsbroek]

Hi @uliw

Please provide some code examples and expected results so I can take a more detailed look!

[uliw]

Ok, attached 3 files 1) phreeqc input file 2) plot of Tot Fe and pH based on phreeqc input file 3) plot of Tot Fe (and S) and pH based on phreeqpython file The python code is below. Note, the code assumes that the kinetics method accepts the "unit" parameter (see my pull request). I have been unable to create a version that operates with the default "mmol" unit.

You can see from the plots that the overall concentrations for Fe and S are similar, but that the resulting pH is very different.

Both files describe the oxidation of 1mol pyrite via Fe+3 and O2 under acidic conditions (pH = 1). The rate function and the available Iron are the same in both files.

pyrite_dissolution.txt test_figure.pdf pyrite_dissolution.pdf

""" Calculate pyrite dissolution
Note: All units must be in mmol!
"""
from phreeqpython import PhreeqPython
import numpy as np
import math
import matplotlib.pyplot as plt

pp = PhreeqPython()

mass = 1  # pyrite mass in mol!
goethite = 2 * mass  # mass in mol"
area = 50  # surface area in m^2/kgw
temperature = 25  # C
pH = 1
time_unit = "hours"  # years, hours, minutes, seconds
run_time = 140  # in time_unit
steps = 100
display_unit = "mol"  # "mol" or "mmol"
plot_fn = "pyrite_dissolution.pdf"

solution = pp.add_solution(
    {
        "units": "mol/kgw",
        "temp": temperature,
        "O(0)": "1 O2(g) -0.7",
        "pH": pH,
        "Fe": goethite,
    }
)

# we work in euqilibrium with atmosphere
solution.equalize(
    ["O2(g)", "CO2(g)"],
    to_si=[-0.7, -3.5],
)
# get time units
time_convert = {
    "seconds": 1,
    "minutes": 60,
    "hours": 60 * 60,
    "days": 60 * 60 * 24,
    "months": 60 * 60 * 24 * 30.5,
    "years": 60 * 60 * 24 * 365,
}
time_unit = time_convert[time_unit]
run_time = run_time * time_unit

# get mol units
mol_convert = {"mol": 1, "mmol": 1000}

# setup result arrays
t = np.array([])
fe = []  # total Fe concentration
ph = []  # pH
s = []  # sulfur

def rate_pyrite(
    sol,  # solution object
    fes2_dissolved,  # amount dissolved pyrite
    m0,  # initial mass in mol/l
    a0,  # surface area m^2/kgw
    v,  # volume, currently not used,
):
    """calculate the pyrite dissolution per time step. Note
    the equation parameter requie that all units must be in mol/kgw
    irrespective of the unit definition in the solution.
    """
    # acid sol parameters
    a1 = 2.82e02
    e1 = 56900
    n1 = -0.500
    n3 = 0.500
    # neutral sol parameters
    a2 = 2.64e05
    e2 = 56900
    n2 = 0.5
    r = 8.31451

    # get solution properties
    pyrite_sr = sol.sr("pyrite")
    hplus_ac = sol.activity("H+", units="mol")
    o2_ac = sol.activity("O2", units="mol")
    fe3_ac = sol.activity("Fe+3", units="mol")
    tk = sol.temperature + 273.15  # T in Kelvin
    m = m0 - fes2_dissolved

    if m0 <= 0:
        sa = a0
    else:
        sa = a0 * ((m / m0) ** 0.67)

    if m < 0:
        rate = 0
    elif m == 0 and pyrite_sr < 1:
        rate = 0
    else:
        acid_rate = a1 * math.e ** (-e1 / (r * tk)) * hplus_ac**n1 * fe3_ac**n3
        neutral_rate = a2 * math.e ** (-e2 / (r * tk)) * o2_ac * n2  #
        rate = (acid_rate + neutral_rate) * (1 - pyrite_sr) * sa

    return rate  # mol

# ----------------- code starts below  ------------------------#

for time, sol in solution.kinetics(
    "Pyrite",
    rate_function=rate_pyrite,
    time=np.linspace(0, run_time, steps),
    m0=mass,
    args=(area, 0),
    unit="mol",
):
    t = np.append(t, time)
    fe.append(sol.total_element("Fe", units="mol"))
    s.append(sol.total_element("S", units="mol"))
    ph.append(-math.log10(sol.activity("H+", units="mol")))

# convert to numpy arrays
fe = np.array(fe)
s = np.array(s)

# convert to display units
c = mol_convert[display_unit]
fe = fe * 1000
s = s * 1000

fig, ax = plt.subplots()
ax2 = ax.twinx()

(l1,) = ax.plot(t / time_unit, fe, label="$\Sigma$ Fe", color="C0")
(l2,) = ax.plot(t / time_unit, s, label="$\Sigma$ S", color="C1")
(l3,) = ax2.plot(t / time_unit, ph, color="C2", label="pH")

# plt.figure(figsize=[10, 5])
ax.legend(handles=[l1, l2, l3])
ax.set_xlabel("Hours")
ax.set_ylabel(f"Concentration [{display_unit}/l]")
plt.title(
    f"FeS$_2$={mass}, Goethite={goethite}[mol/kgw], pH={pH:.1f}, T={solution.temperature:.1f}"
)
ax.grid()
ax2.set_ylabel("pH")
ax2.set_ylim([0, 7])
fig.tight_layout()
fig.savefig(plot_fn)
plt.show()

[AbelHeinsbroek]

Hi @uliw,

I suspect part of the problem stems from the way PhreeqPython currently handles solution operations. At this moment after every command, e.q. add, equalize or saturate the solution is saved:

For example,

s = pp.add_solution({
    'temp': 25,
    'pH': 1
})

s.equalize(['CO2(g)', 'Gibbsite', 'Goethite', 'O2(g)'], to_si=['-3.5', '0', '0', '-0.699'])
s.add('NaOH', 1)

Results in the following steps:

1. Create phreeqc solution
2. Equalize
3. Add 1 mmol of NaOH

and generates the following PHREEQC code:

SOLUTION 0
  temp 25
  pH 1
SAVE SOLUTION 0
END 

USE SOLUTION 0
EQUILIBRIUM PHASES 1 
CO2(g) -3.5 10
Gibbsite 0 10
Goethite 0 10
O2(g) -0.699 10
SAVE SOLUTION 0
END
USE SOLUTION 0
REACTION 1 
NaOH 0.001
1 mol 
SAVE SOLUTION 0
END

Resulting in a pH of .

However, in PHREEQC often the following input is used:

SOLUTION 0
  temp 25
  pH 1

EQUILIBRIUM PHASES 1 
CO2(g) -3.5 10
Gibbsite 0 10
Goethite 0 10
O2(g) -0.699 10

REACTION 1 
NaOH 0.001
1 mol 

SAVE SOLUTION 0
END 

Which has a subtle difference as now the following steps are taken:

1. Create PHREEQC solution
2. Add 1 mmol NaOH whilst being in equilibrium with the given phases.

And results in a different pH of 3.46

I'm trying to come up with an intuitive way of making both different inputs possible. For instance with something like:

s = pp.add_solution({
    'temp': 25,
    'pH': 1
})
s.chain()
s.equalize(['CO2(g)', 'Gibbsite', 'Goethite', 'O2(g)'], to_si=['-3.5', '0', '0', '-0.699'])
s.add('NaOH', 1)
s.end()

[uliw]

I see, so the problem is not caused by the kinetics code, but in the step where the pyrite is added to the solution. Happy test once you have a fix.

[uliw]

I looked at this again, but I think my main problem is that pyrite oxidation depends on pH and that pH itself changes because of pyrite oxidation. However, the kinetics code does the following

sol = solution.copy()
# more code
sol.forget()

so the pH used inside the rate function always uses the initial pH rather than the pH at a given time step.

[AbelHeinsbroek]

Correct. The problem is that the kinetics method doesn't take the 'INCREMENTAL REACTIONS' into account.

However, if we use incremental reactions we don't need to integrate numerically, but can simply use forward Euler to compute the rate. Keep in mind that a PHREEQC kinetics function returns values per second,

""" Calculate pyrite dissolution
Note: All units must be in mmol!
"""
from phreeqpython import PhreeqPython
import numpy as np
import math
import matplotlib.pyplot as plt
import pandas as pd
from scipy.integrate import odeint
pp = PhreeqPython('./phreeqc.dat')
pp.ip.debug = False
pp.ip.run_string('''PHASES
# from https://github.com/HydrogeoIU/PHREEQC-Kinetic-Library
Pyrite
Fe1S2 + 1.0000 H2O = 0.5000 O2    + 2.0000 HS- + 1.0000 Fe+2
    -analytic   7.990927e+003   2.426724e+000   -3.481962e+005  -3.112050e+003  1.492998e+007   -8.237808e-004
    -Vm 23.9400''')
mass = 1  # pyrite mass in mol!
goethite = 2 * mass  # mass in mol"
area = 50  # surface area in m^2/kgw
temperature = 25  # C
pH = 1
time_unit = "hours"  # years, hours, minutes, seconds
run_time = 140  # in time_unit
steps = 100
display_unit = "mol"  # "mol" or "mmol"
plot_fn = "pyrite_dissolution.pdf"

solution = pp.add_solution(
    {
        "units": "mol/kgw",
        "temp": temperature,
        "O(0)": "1 O2(g) -0.7",
        "pH": pH,
        "Fe": goethite,
    }
)

# we work in euqilibrium with atmosphere
solution.equalize(
    ["O2(g)", "CO2(g)"],
    to_si=[-0.7, -3.5],
)
# get time units
time_convert = {
    "seconds": 1,
    "minutes": 60,
    "hours": 60 * 60,
    "days": 60 * 60 * 24,
    "months": 60 * 60 * 24 * 30.5,
    "years": 60 * 60 * 24 * 365,
}
time_unit = time_convert[time_unit]
run_time = run_time * time_unit

# get mol units
mol_convert = {"mol": 1, "mmol": 1000}

# setup result arrays
t = np.array([])
fe = []  # total Fe concentration
ph = []  # pH
s = []  # sulfur

def rate_pyrite(
    sol,  # solution object
    fes2_dissolved,  # amount dissolved pyrite
    m0,  # initial mass in mol/l
    a0,  # surface area m^2/kgw
    v,  # volume, currently not used,
):
    """calculate the pyrite dissolution per time step. Note
    the equation parameter requie that all units must be in mol/kgw
    irrespective of the unit definition in the solution.
    """
    # acid sol parameters
    a1 = 2.82e02
    e1 = 56900
    n1 = -0.500
    n3 = 0.500
    # neutral sol parameters
    a2 = 2.64e05
    e2 = 56900
    n2 = 0.5
    r = 8.31451

    # get solution properties
    pyrite_sr = sol.sr("pyrite")
    hplus_ac = sol.activity("H+", units="mol")
    o2_ac = sol.activity("O2", units="mol")
    fe3_ac = sol.activity("Fe+3", units="mol")
    tk = sol.temperature + 273.15  # T in Kelvin
    m = m0 - fes2_dissolved

    if m0 <= 0:
        sa = a0
    else:
        sa = a0 * ((m / m0) ** 0.67)

    if m < 0:
        rate = 0
    elif m == 0 and pyrite_sr < 1:
        rate = 0
    else:
        acid_rate = a1 * math.e ** (-e1 / r / tk) * hplus_ac**n1 * fe3_ac**n3
        neutral_rate = a2 * math.e ** (-e2 / r / tk) * o2_ac * n2  #
        rate = (acid_rate + neutral_rate) * (1 - pyrite_sr) * sa

    return rate  # mol

# ----------------- code starts below  ------------------------#

sol = solution.copy()
steps = 100
dissolved = 0

for time in np.linspace(0, 5e5, steps):

    rate = rate_pyrite(sol, dissolved, mass, area, 0) * 5e5/steps

    dissolved+=rate

    sol.add('Pyrite', rate, 'mol')
    sol.equalize(
    ["O2(g)", "CO2(g)"],
    to_si=[-0.7, -3.5],
    )

    t = np.append(t, time)
    fe.append(sol.total_element("Fe", units="mol"))
    s.append(sol.total_element("S", units="mol"))
    ph.append(-math.log10(sol.activity("H+", units="mol")))

# convert to numpy arrays
fe = np.array(fe)
s = np.array(s)

# convert to display units
c = mol_convert[display_unit]
fe = fe * 1000
s = s * 1000

fig, ax = plt.subplots()
ax2 = ax.twinx()

(l1,) = ax.plot(t/3600, fe, label="$\Sigma$ Fe", color="C0")
(l2,) = ax.plot(t/3600, s, label="$\Sigma$ S", color="C1")
(l3,) = ax2.plot(t/3600, ph, color="C2", label="pH")

# plt.figure(figsize=[10, 5])
ax.legend(handles=[l1, l2, l3])
ax.set_xlabel("Hours")
ax.set_ylabel(f"Concentration [{display_unit}/l]")
plt.title(
    f"FeS$_2$={mass}, Goethite={goethite}[mol/kgw], pH={pH:.1f}, T={solution.temperature:.1f}"
)
ax.grid()
ax2.set_ylabel("pH")
ax2.set_ylim([2.5, 3.5])
fig.tight_layout()
fig.savefig(plot_fn)
#plt.show()

df = pd.read_csv('pyrite_data.txt', sep='\t')
df = df.iloc[::5,:]
df = df.set_index(' Time(hours)')
df['          pH'].plot(style='x', color='tab:green')
(df['      Tot_Fe']*1000).plot(style='x', ax=ax)
(df['       Tot_S']*1000).plot(style='x', ax=ax)

Which results in the same curve as PHREEQC does:

[uliw]

argh, I did not see the forest for all the trees! Thanks!

[AbelHeinsbroek]

no problem! I still need to write proper documentation for this module, if I even find the time...

Can I close this issue?

[uliw]

sure, go ahead. I'll do a notebook for my students. Let me know if you want a copy for your examples