Comparing Reaktoro with PHREEQC results for Quartz dissolution rate

[leilahashemi]

Description

Hi,
I am quite new to Reaktoro and tried to start with a simple example from PHREEQC for quartz dissolution rate and rewrite in Reaktoro and compare [Si] or pH vs time. In the PHREEQC example, the initial condition involves pure water coming into contact with quartz. However, I am encountering difficulties when trying to replicate this in Reaktoro. Specifically, I haven't been able to remove quartz from the initial chemical system in Reaktoro, to then introduce the reaction at a later stage, as is done in PHREEQC. I would appreciate any help you can give.

Best regards,
Leila

Example

from reaktoro import *

# The surface area of a cube per volume (in m2/m3)
V = 0.16  # Solution volume [L]
A0 = 23.13 # Initial surface of Quartz [m2]
m0 = 158.8 # initial mass of quartz [mole]

def ratefn(props: ChemicalProps):
    aprops = AqueousProps(props)  # we need an AqueousProps object to compute the saturation ratio of the mineral
    k0 = pow(10.0, -13.7)  # the reaction rate constant at 25 oC from Palandri and Kharaka (2004)
    Omega = aprops.saturationRatio("Quartz")  # the current saturation ratio of the mineral omega = IAP/K
    return (A0 / V)  * pow((props.speciesAmount("Quartz")/m0), 0.67) * k0 * (1 - Omega)  # (1 - omega) is the current measure of disequilibrium! omega < 1 implies mineral is undersaturated

db = PhreeqcDatabase("phreeqc.dat") 
system = ChemicalSystem(db,
    AqueousPhase().set(ActivityModelDebyeHuckelPHREEQC()),#("H2O H+ OH- H4SiO4").set(ActivityModelHKF()),
     MineralPhase("Quartz"),
    GeneralReaction("Quartz").setRateModel(ratefn)  # define the reaction and set its reaction rate model!: SiO2 + 2 H2O = H4SiO4
    ) 
 
state = ChemicalState(system)

state.temperature(25.0, "C")
state.pressure(1.0, "bar") 
state.set("H2O", 1.0, "kg") #change to mole later, 0.16 L * 0.99704 g/cm3 =  of solution  
state.add("Quartz", 158.8, "mol") # scalePhaseVolume("Quartz", 1.0, "cm3")  # start with 1 cm3 cube of Quartz crystal 
print("::::::::::::state:::::::::::")  
print(state) 

print("::::::::::::a_props:::::::::::") 
aprops = AqueousProps(state)
print(aprops)  

solver = KineticsSolver(system)  # the chemical kinetics solver

table = Table()  # used to create table of data for later output and plotting

dt = 31536000 / 3  # time step (in seconds): 1year = 3600 * 24 * 365 = 31536000

# Initiate the time stepping for the kinetics modeling
for i in range(16):
    result = solver.solve(state, dt)  # compute the chemical state of the system after it reacted for given time length

    assert result.succeeded(), f"Calculation did not succeed at time step #{i}."

    props = state.props()  # get the current thermodynamic and chemical properties of the system 
    
    aprops = AqueousProps(state)

    table.column("No")   << i  # from seconds to 4 month 
    table.column("Time")   << i*dt / (60.0 * 60 * 24 * 365)  # from seconds to 4 month 
    table.column("pH") << aprops.pH()   # in ?
    table.column("QuartzVolume") << props.phaseProps("Quartz").volume()   # in ?
    table.column("QuartzAmount") << props.speciesAmount("Quartz") * 1e3  # in mmole
    table.column("H4SiO4")    << props.speciesAmount("H4SiO4") * 1e3 # in mmol 

table.save("data_Quartz.txt")

print(table)

print("::::::::::::Phases:::::::::::")   
for phase in system.phases():
    print(phase.name())
    for species in phase.species():
        print(":: " + species.name())

print("::::::::::::a_props:::::::::::")  
print(aprops) mentMolality("Si")         , "molal") 


from reaktplot import *
import pandas as pd

data = pd.read_csv('QuartzDissolution.txt', delim_whitespace=True)  

fig1 = Figure()
fig1.title("H4SiO4 SPECIES AMOUNTS OVER TIME")
fig1.xaxisTitle("Time [year]")
fig1.yaxisTitle("Amount of H4SiO4 [mmol/L]")  
fig1.drawLine(table["Time"], table["H4SiO4"], "H4SiO4") 
fig1.drawLine(data['time'].divide(3.1536e7), data['Si'].multiply(1e3), "Si_PHREEQC") 

fig1.show()
fig1.save("Quartz_H4SiO4.pdf")   

fig3 = Figure()
fig3.title("pH OVER TIME")
fig3.xaxisTitle("Time [year]")
fig3.yaxisTitle("pH") 
fig3.drawLine(table["Time"], table["pH"], "pH")  
fig3.drawLine(data['time'].divide(3.1536e7), data['pH'], "pH_PHREEQC") 
fig3.show()
fig3.save("Quartz_pH.pdf")


### What is your operating system?

Windows

### Additional Information

_No response_

[allanleal]

Hi @leilahashemi , thanks for your question. I'm attaching here a Jupyter Notebook that does this benchmark (I prepared this a while a go, but haven't yet posted on the website).

12 - Chemical Kinetics - Quartz dissolution in water.zip

It produces the following plot:

This is the corresponding PHREEQC script:

# Kinetic dissolution of quartz
#
RATES
 Quartz
 # d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
 # Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
 #  k = 10^-13.7 mol/m2/s (25 C)
 #  A0, initial surface of quartz (m2) recalc's to mol/s
 #  V, solution volume in contact with A0
 -start
  1 A0 = parm(1)
  2 V = parm(2)
  10 rate = (A0 / V) * (m/m0)^0.67 * 10^-13.7 * (1 -  SR("Quartz"))
  20 save rate * time
 -end
KINETICS
Quartz
 -formula SiO2
 -m0  158.8                        # initial moles of quartz
 -parms 23.13  0.16
 -time_step 5 year in 15           # 15 time steps
INCREMENTAL_REACTIONS true

SOLUTION 1
USER_GRAPH
 -chart_title "Quartz dissolution"
 -axis_titles Years "mmol / L"
 -axis_scale x_axis 0 5
 -initial_solutions  true
 -start
 10 graph_x total_time / 3.1536e7  # time in years on x-axis
 20 graph_y tot("Si") * 1e3        # parameter on y-axis
 -end
END

Please let me know if this helps. I can go through your script later if this does not help.

[leilahashemi]

Hi @allanleal
Thank you for replying and sharing the notebook. It works.
Best regards,
Leila

[allanleal]

Excellent! Thanks for reporting back.

[leilahashemi]

Just a follow up update, might be useful to share!

I found my code would work the same as the notebook by moving the reaction to the end of the kinetic reaction loop (so skip reaction in the first step) and also decreasing time-step (or increasing number of steps) helped to get closer match. However, there was still difference in the initial pH in both (e.g., 16 instead of 7 comparing to PHREEQC result)). That could be fixed as follows, where I first did the reaction with zero quartz and then added it in the second step where the loop of kinetic reaction starts :

from reaktoro import *

# The surface area of a cube per volume (in m2/m3)
V = 0.16  # Solution volume [L]
A0 = 23.13 # Initial surface of Quartz [m2]
m0 = 158.8 # initial mass of quartz [mole]

def ratefn(props: ChemicalProps):
    aprops = AqueousProps(props)  # we need an AqueousProps object to compute the saturation ratio of the mineral
    k0 = pow(10.0, -13.7)  # the reaction rate constant at 25 oC from Palandri and Kharaka (2004)
    # q = props.phaseProps("Quartz").volume()  # the amount of a species in the system (in mol)
    Omega = aprops.saturationRatio("Quartz")  # the current saturation ratio of the mineral omega = IAP/K
    rate = (A0 / V)  * pow((props.speciesAmount("Quartz")/m0), 0.67) * k0 * (1 - Omega)  # (1 - omega) is the current measure of disequilibrium! omega < 1 implies mineral is undersaturated
    # print(rate)
    return rate

db = PhreeqcDatabase("phreeqc.dat") 

system = ChemicalSystem(db,
    AqueousPhase().set(ActivityModelDebyeHuckelPHREEQC()), 
    # AqueousPhase("H2O H+ OH- H4SiO4").set(ActivityModelDebyeHuckelPHREEQC()), 
    MineralPhase("Quartz"),
    GeneralReaction("Quartz").setRateModel(ratefn)  # define the reaction and set its reaction rate model!: SiO2 + 2 H2O = H4SiO4
    )  

state = ChemicalState(system)

state.temperature(25.0, "C")
state.pressure(1.0, "bar") 
state.set("H2O", 1.0, "kg") #change to mole later, 0.16 L * 0.99704 g/cm3 =  of solution 
state.add("Quartz", 0, "mol") # scalePhaseVolume("Quartz", 1.0, "cm3")  # start with 1 cm3 cube of Quartz crystal 
print("::::::::::::state:::::::::::")  
print(state)
  
print("::::::::::::a_props:::::::::::") 
aprops = AqueousProps(state)
print(aprops)  

solver = KineticsSolver(system)  # the chemical kinetics solver

table = Table()  # used to create table of data for later output and plotting
timeSteps = 20
dt = 31536000 / timeSteps  # time step (in seconds): 1year = 3600 * 24 * 365 = 31536000
# Initiate the time stepping for the kinetics modeling

i = 0
result = solver.solve(state, dt)  # compute the chemical state of the system after it reacted for given time length
assert result.succeeded(), f"Calculation did not succeed at time step #{i}."
props = state.props()  # get the current thermodynamic and chemical properties of the system 
   
aprops = AqueousProps(state)
    
table.column("No")   << i  # from seconds to 4 month 
table.column("Time")   << i*dt / (60.0 * 60 * 24 * 365)  # from seconds to 4 month 
table.column("pH") << aprops.pH()   # in ?
table.column("QuartzVolume") << props.phaseProps("Quartz").volume()   # in ?
table.column("QuartzAmount") << props.speciesAmount("Quartz") * 1e3  # in mmole
table.column("H4SiO4")    << props.speciesAmount("H4SiO4") * 1e3 # in mmol 
 
# Adding Quartz
state.add("Quartz", 158.8, "mol") # scalePhaseVolume("Quartz", 1.0, "cm3")  # start with 1 cm3 cube of Quartz crystal 
 
for i in range(timeSteps*5):
    result = solver.solve(state, dt)  # compute the chemical state of the system after it reacted for given time length

    assert result.succeeded(), f"Calculation did not succeed at time step #{i+1}."

    props = state.props()  # get the current thermodynamic and chemical properties of the system 
    
    aprops = AqueousProps(state)    

    table.column("No")   << i+1  # from seconds to 4 month 
    table.column("Time")   << (i+1)*dt / (60.0 * 60 * 24 * 365)  # from seconds to 4 month 
    table.column("pH") << aprops.pH()   # in ?
    table.column("QuartzVolume") << props.phaseProps("Quartz").volume()   # in ?
    table.column("QuartzAmount") << props.speciesAmount("Quartz") * 1e3  # in mmole
    table.column("H4SiO4")    << props.speciesAmount("H4SiO4") * 1e3 # in mmol 
 

table.save("data_Quartz.txt")

print(table)

print("::::::::::::Phases:::::::::::")   
for phase in system.phases():
    print(phase.name())
    for species in phase.species():
        print(":: " + species.name())

print("::::::::::::a_props:::::::::::")  
print(aprops) 

from reaktplot import *
import pandas as pd

data = pd.read_csv('QuartzDissolution.txt', delim_whitespace=True) 
# data = pd.read_csv('QuartzDissolution.txt', delim_whitespace=True, skiprows=[1])
# phreeqcdata = {
#     "mSi": [0.02634, 0.04606, 0.06082, 0.07186, 0.08013, 0.08631, 0.09094, 0.09441, 0.09700, 0.09894, 0.10040, 0.10150, 0.10230, 0.10290, 0.10340],
#     "Time": [0.3333, 0.6667, 1.0000, 1.3333, 1.6667, 2.0000, 2.3333, 2.6667, 3.0000, 3.3333, 3.6667, 4.0000, 4.3333, 4.6667, 5.0000],
#     "pH":[6.92915,	6.88802, 6.86163, 6.84378, 6.83132, 6.82244, 6.81602, 6.81134, 0.68079, 6.80536, 6.80348, 6.80208, 6.80105, 6.80027, 6.79969]
# }

fig1 = Figure()
fig1.title("H4SiO4 SPECIES AMOUNTS OVER TIME")
fig1.xaxisTitle("Time [year]")
fig1.yaxisTitle("Amount of H4SiO4 [mmol/L]")  
fig1.drawLine(table["Time"], table["H4SiO4"], "H4SiO4_Reaktoro") 
fig1.drawLine(data['time'].divide(3.1536e7), data['Si'].multiply(1e3), "H4SiO4_PHREEQC") 

fig1.show()
fig1.save("Quartz_H4SiO4.pdf")  
 
fig3 = Figure()
fig3.title("pH OVER TIME")
fig3.xaxisTitle("Time [year]")
fig3.yaxisTitle("pH") 
fig3.drawLine(table["Time"], table["pH"], "pH_Reaktoro")  
fig3.drawLine(data['time'].divide(3.1536e7), data['pH'], "pH_PHREEQC") 
fig3.show()
fig3.save("Quartz_pH.pdf") 

[allanleal]

However, there was still difference in the initial pH in both (e.g., 16 instead of 7 comparing to PHREEQC result)).

Hi @leilahashemi , thanks for sharing your solution. Alternatively, you could also have executed the line below before starting the loop:

solver.precondition(state)

This would have reacted the system with a zero time step, causing the aqueous species to equilibrate, but the mineral to remain inert in the process. solver.precondition(state) is basically the same as solver.solve(state, dt) with dt = 0.0. This would have resolved the strange pH = 16, which is a value corresponding to the solution in disequilibrium.

Here is the documentation of KineticSolver.precondition: https://reaktoro.org/api/classReaktoro_1_1KineticsSolver.html#ae3d484bb706fb8b3e10709c53d276097

[leilahashemi]

Thank you very much @allanleal for the information!

[allanleal]

Thanks, Leila. I'll close this discussion but feel free to repost if you need to.