kinetics with more than one element

[uliw]

I am trying to model pyrite dissolution and set up a rate function that does increase the total Fe, however, the S-concentration does not change. Using the odeint approach, it is evident that I only return the rate for Fe, but I am unclear on how to map a second return value so that it affects the S-concentration. I use the following code

import phreeqpython
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import pprint

def rate_pyrite(pyrite_dissolved, time, sol, A0, m0):
    # acid solution parameters
    a1 = 2.82e02  # A
    e1 = 56900
    n1 = -0.500
    n3 = 0.500
    # neutral solution parameters
    a2 = 2.64e05
    e2 = 56900
    n2 = 0.5
    r = 8.31451

    temp = sol.copy()  # make a temporary copy of the solution
    temp.add("FeS2", pyrite_dissolved[0], "mol")  # add Pyrite to solution
    m = m0 - pyrite_dissolved[0]

    # get solution properties
    pyrite_sr = temp.sr("pyrite")
    hplus_ac = temp.activity("H+", units="mol")
    o2_ac = temp.activity("O2", units="mol")
    fe3plus_ac = temp.activity("Fe+3",  units="mol")
    tk = temp.temperature + 273.15  

    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 * fe3plus_ac**n3
        neutral_rate = a2 * math.e ** (-e2 / (r * tk)) * o2_ac * n2  #
        rate = (acid_rate + neutral_rate) * (1 - pyrite_sr) * sa

    temp.forget()  # cleanup the no longer needed temporary solution
    return rate

# ----------------- code starts below  ------------------------#
pp = phreeqpython.PhreeqPython()

solution = pp.add_solution(
    {
        "units": "mol/kgw",
        "temp": 30,
        "O(0)": "1 O2(g) -0.7",
        "pH": 1,
        "Fe": 0.1,
        "S" : 0,
    }
)
solution.equalize(
    ["O2(g)", "CO2(g)"],
    to_si=[-0.7, -3.5],
)

hours = 60 * 60
time = np.linspace(0, 30 * hours, 100)
surface_area = 50
mass = 0.05  # mol

# solve differential equation
yy = odeint(
    rate_pyrite,
    0,
    time,
    args=(
        solution,
        surface_area,
        mass,
    ),
)

pprint.pp(solution.species)
plt.plot(time / hours, yy * 1e3)
plt.xlabel("Hours")
plt.ylabel("mmol/l")
plt.title("Pyrite Dissolution")
plt.grid()
plt.show()

[uliw]

never mind, I had a look at the source, and now understand that the above only provides the rates, but that the actual chemisrty is handled afterwards by

y = np.insert(np.diff(yy[:,0]), 0, 0)
for i in range(len(time)):
    t = time[i]
    self.add(element, y[i])
    yield(t, self)