Question

[Khaled]

Hi everyone,

I have some questions related to the code below. I wanted to find the rate of precipitation, but the graphs are not logical it seems there's a mistake? can anyone guide what's wrong, the amount precipitated should be asymptotic with time? Also, I tried taking the values manually from the results page and added the amounts and it gave me a linear relation between amount and time for aragonite, and also, I want to add to the code some part to give me the thickness of the precipitation. I want an explanation. Thanks

RATES
Aragonite
-start
 1 rem PARM(1) = log(specific area), #log(m^2 per mole calcite)
 2 rem PARM(2) =  (M/M0)
 40 log_k2= log (1.9e-2) + ((71.2/0.008314)*((1/(25+273.15))-(1/(45+273.15))))   
 50 lograte = log_k2 + log(PARM(1)) + log((M/M0)^PARM(2))
 60 moles = exp(lograte) * TIME * (1 - SR("Aragonite"))                       
 200 SAVE moles
210 PRINT moles, exp(lograte), TIME, log_k2, PARM(1)
-end
Gypsum
-start
 1 rem PARM(1) = log(specific area), #log(m^2 per mole siderite)
 2 rem PARM(2) = (M/M0)
 40 log_k2= log(1e-9) + ((35/0.008314)*((1/(25+273.15))-(1/(45+273.15))))   
 50 lograte = log_k2 + log(PARM(1)) + log((M/M0)^PARM(2))                         
 60 moles = exp(lograte) * TIME *(1-SR("Gypsum"))
 200 SAVE moles
-end

KINETICS 1
Aragonite
       -formula CaCO3     
       -m0    1e-3              #initial mole of the mineral
       -parms 300.1  0           #surface area (m2/mol)
Gypsum
       -formula CaSO4.2H2O       
       -m0     1e-3              #initial mole of the mineral     
       -parms 166.37  0          #surface area (m2/mol)
-time  86400 in 24
-cvode
INCREMENTAL_REACTIONS true

SOLUTION 1
    temp     65
    pH        7
    pe        4
    redox     pe
    units     ppm
    density   1
   
    Br        97.5
    C(4)      195
    Ca        690
    Cl        39832.5
    Co        0.1
    Fe        0.0375
    Mg        2792.5
    Na        16500
    S(6)      5187.5
    Si        1.125
    N(-3)     0.5625
    -water    1 # kg

USER_GRAPH 1 Aragonite mole number as a function of time
   -headings time Aragonite SI(Aragonite)
   -axis_titles "time (Hours)" "Aragonite (mmol)" "SI"
   #-axis_scale y_axis 0 5 1
   #-axis_scale x_axis 0 30 1
   -axis_scale sy_axis -1 1
   -start
   10 graph_x total_time/3600
   20 graph_y KIN("Aragonite")*1000
   30 graph_sy SI("Aragonite")
-end
USER_GRAPH 2 Gypsum mole number as a function of time
   -headings time Gypsum SI(gyp)
   -axis_titles "time (Hours)" "Gypsum (mmol)" "SI"
   #-axis_scale y_axis 0 5 1
   #-axis_scale x_axis 0 30 1
   -start
   10 graph_x total_time/3600
   20 graph_y KIN("Gypsum")*1000
   30 graph_sy SI("Gypsum")
-end
END

[dlparkhurst]

You have a high-ionic strength water, where the activity of water is lower. In this case, anhydrite is the stable phase relative to gypsum.

Your kinetic reaction with anhydrite is fast, so anhydrite remains near equilibrium while gypsum dissolves. It takes about an hour for all of the gypsum to dissolve, at which point the reactions stop. All of the gypsum has dissolved and reprecipitated as anhydrite.

[Khaled]

but am asking about Aragonite not anhydrite? so what is the problem with the code I sent, because the rate should be asymptotic not linear? How can I edit this in the code. Thanks 

[dlparkhurst]

Sorry, I was not reading carefully, but the basics were right. Aragonite is in equilibrium and precipitates throughout. Gypsum is undersaturated and dissolves throughout. Gypsum dissolves at a constant rate until gypsum is totally removed in about an hour.

If you add 25 mmol of gypsum in your KINETIC definition and run for 100 hours, then there is sufficient gypsum to reach equilibrium. The SI(gyp) asymptotically reaches equilibrium.

[Khaled]

What if I want to make this code for flow system, I mean how to make the same input but flowing can it be done using phreeqc or better to couple it with another software?

Thanks

[dlparkhurst]

IPhreeqc and IPhreeqcCOM are packages that allow you to use PHREEQC in scripting languages and programming languages C++, Fortrans, and C.

If you want to use PHREEQC with a flow and transport code PhreeqcRM is a package that is designed to have cells corresponding to your model cells. The basic idea is to alternate between PhreeqcRM calculations and transporting with a flow code. The updated concentrations would be passed back and forth at each time step. PhreeqcRM has parallelization with either OpenMPI or MPI.

[Khaled]

Thanks for your reply

I was searching and found a software called SYSCAD, do you have an idea about it if is it better than PHREEQRM? Or what do you think the best option to link phreeqc with for flow systems?

Thanks

[dlparkhurst]

Your choice.